[contrasting] add spatio temporal test for return level confidence interval for polynomial models

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

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

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):

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
                                )

test/test_extreme_fit/test_estimator/test_spatio_temporal_estimator/test_gev_spatio_temporal_polynomial_extremes_mle.py

+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()