diff --git a/extreme_fit/distribution/gev/gev_params.py b/extreme_fit/distribution/gev/gev_params.py index f1e0864bb3a2a67fadb42663555c355245f9cb9d..f2ae17420cf8649414ae4e8656b63b6190b46f3b 100644 --- a/extreme_fit/distribution/gev/gev_params.py +++ b/extreme_fit/distribution/gev/gev_params.py @@ -135,6 +135,15 @@ class GevParams(AbstractExtremeParams): cls.SHAPE: 'zeta', }[param_name] + @classmethod + def greek_letter_from_param_name_confidence_interval(cls, param_name): + assert param_name in cls.PARAM_NAMES + return { + cls.LOC: 'mu', + cls.SCALE: 'sigma', + cls.SHAPE: 'xi', + }[param_name] + @classmethod def full_name_from_param_name(cls, param_name): assert param_name in cls.PARAM_NAMES diff --git a/extreme_fit/distribution/gev/main_fevd_mle.R b/extreme_fit/distribution/gev/main_fevd_mle.R index e39cd52b203972b114c074549b144f344efc6269..a43fda6bdc5810d7ee6213faa6b6f3c9d41782d6 100644 --- a/extreme_fit/distribution/gev/main_fevd_mle.R +++ b/extreme_fit/distribution/gev/main_fevd_mle.R @@ -23,7 +23,7 @@ colnames(coord) = c("T") coord = data.frame(coord, stringsAsFactors = TRUE) # res = fevd_fixed(x_gev, data=coord, method='MLE', verbose=TRUE, use.phi=FALSE) # res = fevd_fixed(x_gev, data=coord, location.fun= ~T, scale.fun= ~T, method='MLE', type="GEV", verbose=FALSE, use.phi=FALSE) -res = fevd_fixed(x_gev, data=coord, location.fun= ~poly(T, 2), method='MLE', type="GEV", verbose=FALSE, use.phi=FALSE) +res = fevd_fixed(x_gev, data=coord, shape.fun= ~poly(T, 2), method='MLE', type="GEV", verbose=FALSE, use.phi=FALSE) print(res) # Some display for the results diff --git a/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_extract_eurocode_return_level.py b/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_extract_eurocode_return_level.py index 38a173b18cc4861de006f69f8bf19d6d75947cbd..6ddb75c8f1bc10d693e6c697f510d46a974ac6ba 100644 --- a/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_extract_eurocode_return_level.py +++ b/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_extract_eurocode_return_level.py @@ -43,7 +43,12 @@ class ExtractEurocodeReturnLevelFromCiMethod(AbstractExtractEurocodeReturnLevel) @property def transformed_temporal_covariate(self): - return self.estimator.dataset.coordinates.transformation.transform_float(self.temporal_covariate) + if isinstance(self.temporal_covariate, (float, int)): + return self.estimator.dataset.coordinates.transformation.transform_float(self.temporal_covariate) + elif isinstance(self.temporal_covariate, np.ndarray): + return self.estimator.dataset.coordinates.transformation.transform_array(self.temporal_covariate) + else: + raise NotImplementedError @cached_property def confidence_interval_method(self): diff --git a/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_result_from_extremes.py b/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_result_from_extremes.py index 4870914970d9502bd0266d3d6967617ea5653f02..00b50eb60de2656c06ae58aefbf41956c7a6e4cf 100644 --- a/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_result_from_extremes.py +++ b/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_result_from_extremes.py @@ -60,9 +60,19 @@ class AbstractResultFromExtremes(AbstractResultFromModelFit): 'type': r.c("return.level") } if self.param_name_to_dim: - d = {GevParams.greek_letter_from_param_name(param_name) + '1': r.c(transformed_temporal_covariate) for - param_name in self.param_name_to_dim.keys()} - print(d) + if isinstance(transformed_temporal_covariate, (int, float, np.int, np.float, np.int64)): + d = {GevParams.greek_letter_from_param_name_confidence_interval(param_name) + '1': r.c(transformed_temporal_covariate) for + param_name in self.param_name_to_dim.keys()} + elif isinstance(transformed_temporal_covariate, np.ndarray): + d = {} + for param_name in self.param_name_to_dim: + for coordinate_idx, _ in self.param_name_to_dim[param_name]: + idx_str = str(coordinate_idx + 1) + d2 = {GevParams.greek_letter_from_param_name_confidence_interval(param_name) + idx_str: r.c(transformed_temporal_covariate)} + d.update(d2) + else: + raise NotImplementedError + kwargs = { "vals": r.list(**d ) diff --git a/test/test_extreme_fit/test_estimator/test_spatio_temporal_estimator/__init__.py b/test/test_extreme_fit/test_estimator/test_spatio_temporal_estimator/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/test/test_extreme_fit/test_estimator/test_spatio_temporal_estimator/test_gev_spatio_temporal_polynomial_extremes_mle.py b/test/test_extreme_fit/test_estimator/test_spatio_temporal_estimator/test_gev_spatio_temporal_polynomial_extremes_mle.py new file mode 100644 index 0000000000000000000000000000000000000000..3db108e7c76cfd5ac00c636a79b07ea53ce28cdb --- /dev/null +++ b/test/test_extreme_fit/test_estimator/test_spatio_temporal_estimator/test_gev_spatio_temporal_polynomial_extremes_mle.py @@ -0,0 +1,107 @@ +import unittest + +import numpy as np +import pandas as pd + +from extreme_fit.distribution.gev.gev_params import GevParams +from extreme_fit.estimator.margin_estimator.utils import fitted_linear_margin_estimator_short +from extreme_fit.model.margin_model.polynomial_margin_model.gev_altitudinal_models import StationaryAltitudinal +from extreme_fit.model.margin_model.polynomial_margin_model.models_based_on_pariwise_analysis.gev_with_constant_shape_wrt_altitude import \ + AltitudinalShapeConstantTimeLocationLinear +from extreme_fit.model.margin_model.polynomial_margin_model.models_based_on_pariwise_analysis.gev_with_linear_shape_wrt_altitude import \ + AltitudinalShapeLinearTimeLocationLinear, AltitudinalShapeLinearTimeLocScaleLinear, \ + AltitudinalShapeLinearTimeStationary +from extreme_fit.model.margin_model.polynomial_margin_model.polynomial_margin_model import \ + NonStationaryQuadraticLocationModel, \ + NonStationaryQuadraticScaleModel, NonStationaryQuadraticLocationGumbelModel, NonStationaryQuadraticScaleGumbelModel +from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import \ + ConfidenceIntervalMethodFromExtremes +from extreme_fit.model.result_from_model_fit.result_from_extremes.eurocode_return_level_uncertainties import \ + EurocodeConfidenceIntervalFromExtremes +from extreme_trend.abstract_gev_trend_test import fitted_linear_margin_estimator +from extreme_fit.model.margin_model.utils import \ + MarginFitMethod +from extreme_fit.model.utils import r, set_seed_r +from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates +from spatio_temporal_dataset.coordinates.spatial_coordinates.abstract_spatial_coordinates import \ + AbstractSpatialCoordinates +from spatio_temporal_dataset.coordinates.spatio_temporal_coordinates.abstract_spatio_temporal_coordinates import \ + AbstractSpatioTemporalCoordinates +from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_coordinates import \ + AbstractTemporalCoordinates +from spatio_temporal_dataset.coordinates.temporal_coordinates.generated_temporal_coordinates import \ + ConsecutiveTemporalCoordinates +from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset +from spatio_temporal_dataset.spatio_temporal_observations.abstract_spatio_temporal_observations import \ + AbstractSpatioTemporalObservations + + +class TestGevTemporalQuadraticExtremesMle(unittest.TestCase): + + def setUp(self) -> None: + nb_data = 100 + set_seed_r() + r(""" + N <- {} + loc = 0; scale = 1; shape <- 0.1 + x_gev <- rgev(N, loc = loc, scale = scale, shape = shape) + start_loc = 0; start_scale = 1; start_shape = 1 + """.format(nb_data)) + + # Compute coordinates + altitudes = [300, 600] + temporal_coordinates = ConsecutiveTemporalCoordinates.from_nb_temporal_steps(nb_data) + spatial_coordinates = AbstractSpatialCoordinates.from_list_x_coordinates(altitudes) + self.coordinates = AbstractSpatioTemporalCoordinates.from_spatial_coordinates_and_temporal_coordinates( + spatial_coordinates, + temporal_coordinates) + + # Compute observations + a = np.array(r['x_gev']) + data = np.concatenate([a, a], axis=0) + df2 = pd.DataFrame(data=data, index=self.coordinates.index) + observations = AbstractSpatioTemporalObservations(df_maxima_gev=df2) + + self.dataset = AbstractDataset(observations=observations, coordinates=self.coordinates) + self.fit_method = MarginFitMethod.extremes_fevd_mle + + def function_test_gev_spatio_temporal_margin_fit_non_stationary_quadratic(self, model_class): + # Create estimator + estimator = fitted_linear_margin_estimator_short(model_class=model_class, + dataset=self.dataset, + fit_method=self.fit_method) + + # Assert that indicators are correctly computed + self.assertAlmostEqual(estimator.result_from_model_fit.nllh, estimator.nllh()) + self.assertAlmostEqual(estimator.result_from_model_fit.aic, estimator.aic()) + self.assertAlmostEqual(estimator.result_from_model_fit.bic, estimator.bic()) + # Assert we can compute the return level + covariate_for_return_level = np.array([400, 25])[::1] + confidence_interval = EurocodeConfidenceIntervalFromExtremes.from_estimator_extremes(estimator, + ci_method=ConfidenceIntervalMethodFromExtremes.ci_mle, + temporal_covariate=covariate_for_return_level) + return_level = estimator.function_from_fit.get_params(covariate_for_return_level).return_level(50) + self.assertAlmostEqual(return_level, confidence_interval.mean_estimate) + self.assertFalse(np.isnan(confidence_interval.confidence_interval[0])) + self.assertFalse(np.isnan(confidence_interval.confidence_interval[1])) + + def test_gev_spatio_temporal_margin_fit_1(self): + self.function_test_gev_spatio_temporal_margin_fit_non_stationary_quadratic(StationaryAltitudinal) + + def test_gev_spatio_temporal_margin_fit_1_bis(self): + self.function_test_gev_spatio_temporal_margin_fit_non_stationary_quadratic(AltitudinalShapeLinearTimeStationary) + + # def test_gev_spatio_temporal_margin_fit_2(self): + # self.function_test_gev_spatio_temporal_margin_fit_non_stationary_quadratic( + # AltitudinalShapeConstantTimeLocationLinear) + # + # def test_gev_spatio_temporal_margin_fit_3(self): + # self.function_test_gev_spatio_temporal_margin_fit_non_stationary_quadratic( + # AltitudinalShapeLinearTimeLocationLinear) + # + # def test_gev_spatio_temporal_margin_fit_4(self): + # self.function_test_gev_spatio_temporal_margin_fit_non_stationary_quadratic( + # AltitudinalShapeLinearTimeLocScaleLinear) + +if __name__ == '__main__': + unittest.main()