[contrasting] add Gumbel altitudinal models. But so far only the first three...

[contrasting] add Gumbel altitudinal models. But so far only the first three models are passing tests (models with quadratic & cross terms are crashing)

extreme_fit/model/margin_model/polynomial_margin_model/altitudinal_models.py → extreme_fit/model/margin_model/polynomial_margin_model/gev_altitudinal_models.py

@@ -5,6 +5,7 @@ from extreme_fit.model.margin_model.polynomial_margin_model.spatio_temporal_poly
 
 
 class AbstractAltitudinalModel(AbstractSpatioTemporalPolynomialModel):
+    DISTRIBUTION_STR = 'Gev'
 
     def load_margin_function(self, param_name_to_dims=None):
         return super().load_margin_function(self.param_name_to_list_dim_and_degree_for_margin_function)
@@ -22,7 +23,7 @@ class AbstractAltitudinalModel(AbstractSpatioTemporalPolynomialModel):
 
     def dim_to_str_number(self, param_name, dim):
         list_dim_and_degree = self.param_name_to_list_dim_and_degree[param_name]
-        dims  = [d for d, _ in list_dim_and_degree]
+        dims = [d for d, _ in list_dim_and_degree]
         if dim not in dims:
             return '0'
         else:
@@ -31,7 +32,7 @@ class AbstractAltitudinalModel(AbstractSpatioTemporalPolynomialModel):
 
     @property
     def name_str(self):
-        name = 'Gev'
+        name = self.DISTRIBUTION_STR
         name += self.dim_to_str_number(GevParams.LOC, self.coordinates.idx_temporal_coordinates)
         name += self.dim_to_str_number(GevParams.SCALE, self.coordinates.idx_temporal_coordinates)
         if isinstance(self, AbstractAddCrossTermForLocation):
@@ -49,6 +50,16 @@ class StationaryAltitudinal(AbstractAltitudinalModel):
         }
 
 
+class NonStationaryAltitudinalScaleLinear(AbstractAltitudinalModel):
+
+    @property
+    def param_name_to_list_dim_and_degree(self):
+        return {
+            GevParams.LOC: [(self.coordinates.idx_x_coordinates, 1)],
+            GevParams.SCALE: [(self.coordinates.idx_x_coordinates, 1), (self.coordinates.idx_temporal_coordinates, 1)]
+        }
+
+
 class NonStationaryAltitudinalLocationLinear(AbstractAltitudinalModel):
 
     @property
@@ -124,3 +135,8 @@ class NonStationaryAltitudinalLocationLinearScaleLinearCrossTermForLocation(Abst
 class NonStationaryAltitudinalLocationQuadraticScaleLinearCrossTermForLocation(AbstractAddCrossTermForLocation,
                                                                                NonStationaryAltitudinalLocationQuadraticScaleLinear):
     pass
+
+
+class NonStationaryAltitudinalScaleLinearCrossTermForLocation(AbstractAddCrossTermForLocation,
+                                                              NonStationaryAltitudinalScaleLinear):
+    pass

extreme_fit/model/margin_model/polynomial_margin_model/gumbel_altitudinal_models.py

+from extreme_fit.distribution.gev.gev_params import GevParams
+from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import GumbelTemporalModel
+from extreme_fit.model.margin_model.polynomial_margin_model.gev_altitudinal_models import AbstractAltitudinalModel, \
+    StationaryAltitudinal, NonStationaryAltitudinalScaleLinear, NonStationaryAltitudinalLocationLinear, \
+    NonStationaryAltitudinalLocationQuadratic, NonStationaryAltitudinalLocationLinearScaleLinear, \
+    NonStationaryAltitudinalLocationQuadraticScaleLinear, NonStationaryAltitudinalScaleLinearCrossTermForLocation, \
+    NonStationaryCrossTermForLocation, NonStationaryAltitudinalLocationLinearCrossTermForLocation, \
+    NonStationaryAltitudinalLocationQuadraticCrossTermForLocation, \
+    NonStationaryAltitudinalLocationLinearScaleLinearCrossTermForLocation, \
+    NonStationaryAltitudinalLocationQuadraticScaleLinearCrossTermForLocation, AbstractAddCrossTermForLocation
+from extreme_fit.model.margin_model.polynomial_margin_model.polynomial_margin_model import PolynomialMarginModel
+from extreme_fit.model.margin_model.polynomial_margin_model.spatio_temporal_polynomial_model import \
+    AbstractSpatioTemporalPolynomialModel
+from extreme_fit.model.margin_model.utils import MarginFitMethod
+from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
+
+
+class AbstractGumbelAltitudinalModel(AbstractAltitudinalModel):
+    DISTRIBUTION_STR = 'Gum'
+
+    def __init__(self, coordinates: AbstractCoordinates, params_user=None, starting_point=None,
+                 fit_method=MarginFitMethod.extremes_fevd_mle, nb_iterations_for_bayesian_fit=5000,
+                 params_initial_fit_bayesian=None, params_class=GevParams, max_degree=2):
+        super().__init__(coordinates, params_user, starting_point, fit_method, nb_iterations_for_bayesian_fit,
+                         params_initial_fit_bayesian, "Gumbel", params_class, max_degree)
+
+    @property
+    def nb_params(self):
+        return super().nb_params - 1
+
+
+class StationaryGumbelAltitudinal(AbstractGumbelAltitudinalModel, StationaryAltitudinal):
+    pass
+
+
+class NonStationaryGumbelAltitudinalScaleLinear(AbstractGumbelAltitudinalModel, NonStationaryAltitudinalScaleLinear):
+    pass
+
+
+class NonStationaryGumbelAltitudinalLocationLinear(AbstractGumbelAltitudinalModel,
+                                                   NonStationaryAltitudinalLocationLinear):
+    pass
+
+
+class NonStationaryGumbelAltitudinalLocationQuadratic(AbstractGumbelAltitudinalModel,
+                                                      NonStationaryAltitudinalLocationQuadratic):
+    pass
+
+
+class NonStationaryGumbelAltitudinalLocationLinearScaleLinear(AbstractGumbelAltitudinalModel,
+                                                              NonStationaryAltitudinalLocationLinearScaleLinear):
+    pass
+
+
+class NonStationaryGumbelAltitudinalLocationQuadraticScaleLinear(AbstractGumbelAltitudinalModel,
+                                                                 NonStationaryAltitudinalLocationQuadraticScaleLinear):
+    pass
+
+
+# Add cross terms
+
+
+class NonStationaryGumbelCrossTermForLocation(AbstractAddCrossTermForLocation, AbstractGumbelAltitudinalModel, StationaryAltitudinal):
+    pass
+
+
+class NonStationaryGumbelAltitudinalLocationLinearCrossTermForLocation(AbstractAddCrossTermForLocation, AbstractGumbelAltitudinalModel,
+                                                                       NonStationaryAltitudinalLocationLinear):
+    pass
+
+
+class NonStationaryGumbelAltitudinalLocationQuadraticCrossTermForLocation(AbstractGumbelAltitudinalModel,
+                                                                          NonStationaryAltitudinalLocationQuadraticCrossTermForLocation):
+    pass
+
+
+class NonStationaryGumbelAltitudinalLocationLinearScaleLinearCrossTermForLocation(AbstractGumbelAltitudinalModel,
+                                                                                  NonStationaryAltitudinalLocationLinearScaleLinearCrossTermForLocation):
+    pass
+
+
+class NonStationaryGumbelAltitudinalLocationQuadraticScaleLinearCrossTermForLocation(AbstractGumbelAltitudinalModel,
+                                                                                     NonStationaryAltitudinalLocationQuadraticScaleLinearCrossTermForLocation):
+    pass
+
+
+class NonStationaryGumbelAltitudinalScaleLinearCrossTermForLocation(AbstractGumbelAltitudinalModel,
+                                                                    NonStationaryAltitudinalScaleLinearCrossTermForLocation):
+    pass

extreme_fit/model/margin_model/polynomial_margin_model/utils.py

-from extreme_fit.model.margin_model.polynomial_margin_model.altitudinal_models import StationaryAltitudinal, \
+from extreme_fit.model.margin_model.polynomial_margin_model.gev_altitudinal_models import StationaryAltitudinal, \
     NonStationaryAltitudinalLocationLinear, NonStationaryAltitudinalLocationQuadratic, \
     NonStationaryAltitudinalLocationLinearScaleLinear, NonStationaryAltitudinalLocationQuadraticScaleLinear, \
     NonStationaryCrossTermForLocation, NonStationaryAltitudinalLocationLinearCrossTermForLocation, \
     NonStationaryAltitudinalLocationQuadraticCrossTermForLocation, \
     NonStationaryAltitudinalLocationLinearScaleLinearCrossTermForLocation, \
-    NonStationaryAltitudinalLocationQuadraticScaleLinearCrossTermForLocation
+    NonStationaryAltitudinalLocationQuadraticScaleLinearCrossTermForLocation, NonStationaryAltitudinalScaleLinear, \
+    NonStationaryAltitudinalScaleLinearCrossTermForLocation
+from extreme_fit.model.margin_model.polynomial_margin_model.gumbel_altitudinal_models import \
+    StationaryGumbelAltitudinal, NonStationaryGumbelAltitudinalScaleLinear, \
+    NonStationaryGumbelAltitudinalLocationLinear, NonStationaryGumbelAltitudinalLocationQuadratic, \
+    NonStationaryGumbelAltitudinalLocationLinearScaleLinear, NonStationaryGumbelAltitudinalLocationQuadraticScaleLinear, \
+    NonStationaryGumbelAltitudinalLocationLinearCrossTermForLocation, \
+    NonStationaryGumbelAltitudinalLocationQuadraticCrossTermForLocation, \
+    NonStationaryGumbelAltitudinalLocationLinearScaleLinearCrossTermForLocation, \
+    NonStationaryGumbelAltitudinalLocationQuadraticScaleLinearCrossTermForLocation, \
+    NonStationaryGumbelAltitudinalScaleLinearCrossTermForLocation, NonStationaryGumbelCrossTermForLocation
 from extreme_fit.model.margin_model.polynomial_margin_model.spatio_temporal_polynomial_model import \
     NonStationaryLocationSpatioTemporalLinearityModel1, NonStationaryLocationSpatioTemporalLinearityModel2, \
     NonStationaryLocationSpatioTemporalLinearityModel3, NonStationaryLocationSpatioTemporalLinearityModel4, \
     NonStationaryLocationSpatioTemporalLinearityModel5, NonStationaryLocationSpatioTemporalLinearityModelAssertError1, \
     NonStationaryLocationSpatioTemporalLinearityModelAssertError2, \
     NonStationaryLocationSpatioTemporalLinearityModelAssertError3, NonStationaryLocationSpatioTemporalLinearityModel6
-ALTITUDINAL_MODELS = [
-                         StationaryAltitudinal,
-                         NonStationaryAltitudinalLocationLinear,
-                         NonStationaryAltitudinalLocationQuadratic,
-                         NonStationaryAltitudinalLocationLinearScaleLinear,
-                         NonStationaryAltitudinalLocationQuadraticScaleLinear,
 
-                         NonStationaryCrossTermForLocation,
-                         NonStationaryAltitudinalLocationLinearCrossTermForLocation,
-                         NonStationaryAltitudinalLocationQuadraticCrossTermForLocation,
-                         NonStationaryAltitudinalLocationLinearScaleLinearCrossTermForLocation,
-                        NonStationaryAltitudinalLocationQuadraticScaleLinearCrossTermForLocation,
-                     ][:]
+ALTITUDINAL_GEV_MODELS = [
+                             StationaryAltitudinal,
+                             NonStationaryAltitudinalScaleLinear,
+                             NonStationaryAltitudinalLocationLinear,
+                             NonStationaryAltitudinalLocationQuadratic,
+                             NonStationaryAltitudinalLocationLinearScaleLinear,
+                             NonStationaryAltitudinalLocationQuadraticScaleLinear,
 
+                             NonStationaryCrossTermForLocation,
+                             NonStationaryAltitudinalScaleLinearCrossTermForLocation,
+                             NonStationaryAltitudinalLocationLinearCrossTermForLocation,
+                             NonStationaryAltitudinalLocationQuadraticCrossTermForLocation,
+                             NonStationaryAltitudinalLocationLinearScaleLinearCrossTermForLocation,
+                             NonStationaryAltitudinalLocationQuadraticScaleLinearCrossTermForLocation,
+                         ][:]
 
+ALTITUDINAL_GUMBEL_MODELS = [
+                                StationaryGumbelAltitudinal,
+                                NonStationaryGumbelAltitudinalScaleLinear,
+                                NonStationaryGumbelAltitudinalLocationLinear,
+                                NonStationaryGumbelAltitudinalLocationQuadratic,
+                                NonStationaryGumbelAltitudinalLocationLinearScaleLinear,
+                                NonStationaryGumbelAltitudinalLocationQuadraticScaleLinear,
 
+                                NonStationaryGumbelCrossTermForLocation,
+                                NonStationaryGumbelAltitudinalLocationLinearCrossTermForLocation,
+                                NonStationaryGumbelAltitudinalLocationQuadraticCrossTermForLocation,
+                                NonStationaryGumbelAltitudinalLocationLinearScaleLinearCrossTermForLocation,
+                                NonStationaryGumbelAltitudinalLocationQuadraticScaleLinearCrossTermForLocation,
+                                NonStationaryGumbelAltitudinalScaleLinearCrossTermForLocation,
+                            ][:]
 
 MODELS_THAT_SHOULD_RAISE_AN_ASSERTION_ERROR = [NonStationaryLocationSpatioTemporalLinearityModelAssertError1,
                                                NonStationaryLocationSpatioTemporalLinearityModelAssertError2,

projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py

@@ -10,19 +10,22 @@ from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranS
     SafranSnowfall5Days, SafranSnowfall7Days, SafranPrecipitation1Day, SafranPrecipitation3Days, \
     SafranPrecipitation5Days, SafranPrecipitation7Days
 from extreme_data.meteo_france_data.scm_models_data.utils import Season
-from extreme_fit.model.margin_model.polynomial_margin_model.utils import ALTITUDINAL_MODELS
+from extreme_fit.model.margin_model.polynomial_margin_model.utils import ALTITUDINAL_GEV_MODELS, \
+    ALTITUDINAL_GUMBEL_MODELS
 from projects.altitude_spatial_model.altitudes_fit.altitudes_studies import AltitudesStudies
 from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitudes_studies_visualizer_for_non_stationary_models import \
     AltitudesStudiesVisualizerForNonStationaryModels
 
 
 def plot_altitudinal_fit(studies, massif_names=None):
-    model_classes = ALTITUDINAL_MODELS
+    model_classes = ALTITUDINAL_GEV_MODELS + ALTITUDINAL_GUMBEL_MODELS
+    # model_classes = ALTITUDINAL_GUMBEL_MODELS
     visualizer = AltitudesStudiesVisualizerForNonStationaryModels(studies=studies,
                                                                   model_classes=model_classes,
                                                                   massif_names=massif_names,
                                                                   show=False)
     # Plot the results for the model that minimizes the individual aic
+    OneFoldFit.best_estimator_minimizes_mean_aic = False
     visualizer.plot_moments()
     visualizer.plot_shape_map()
 
@@ -31,13 +34,14 @@ def plot_altitudinal_fit(studies, massif_names=None):
     model_class_to_aic_scores = {model_class: [] for model_class in model_classes}
     for one_fold_fit in visualizer.massif_name_to_one_fold_fit.values():
         for model_class, estimator in one_fold_fit.model_class_to_estimator_with_finite_aic.items():
-            print(estimator.n())
             aic_score_normalized = estimator.aic() / estimator.n()
             model_class_to_aic_scores[model_class].append(aic_score_normalized)
     model_class_to_mean_aic_score = {model_class: np.array(aic_scores).mean()
                                      for model_class, aic_scores in model_class_to_aic_scores.items()}
+    print(model_class_to_mean_aic_score)
     sorted_model_class = sorted(model_classes, key=lambda m: model_class_to_mean_aic_score[m])
     best_model_class_for_mean_aic = sorted_model_class[0]
+    print(best_model_class_for_mean_aic)
     for one_fold_fit in visualizer.massif_name_to_one_fold_fit.values():
         one_fold_fit.best_estimator_class_for_mean_aic = best_model_class_for_mean_aic
 
@@ -58,8 +62,8 @@ def plot_moments(studies, massif_names=None):
 
 def main():
     altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][4:7]
-    altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][:]
-    altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900]
+    # altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][:]
+    # altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900]
     study_classes = [SafranSnowfall1Day, SafranSnowfall3Days, SafranSnowfall5Days, SafranSnowfall7Days][:2]
     study_classes = [SafranPrecipitation1Day, SafranPrecipitation3Days, SafranPrecipitation5Days,
                      SafranPrecipitation7Days][:]

projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py

@@ -5,7 +5,9 @@ from cached_property import cached_property
 from scipy.stats import chi2
 
 from extreme_fit.estimator.margin_estimator.utils import fitted_linear_margin_estimator_short
-from extreme_fit.model.margin_model.polynomial_margin_model.altitudinal_models import StationaryAltitudinal
+from extreme_fit.model.margin_model.polynomial_margin_model.gev_altitudinal_models import StationaryAltitudinal
+from extreme_fit.model.margin_model.polynomial_margin_model.gumbel_altitudinal_models import \
+    StationaryGumbelAltitudinal, AbstractGumbelAltitudinalModel
 from extreme_fit.model.margin_model.utils import MarginFitMethod
 from root_utils import classproperty
 
@@ -137,9 +139,10 @@ class OneFoldFit(object):
 
     @property
     def stationary_estimator(self):
-        for estimator in self.sorted_estimators_with_finite_aic:
-            if isinstance(estimator.margin_model, StationaryAltitudinal):
-                return estimator
+        if isinstance(self.best_estimator.margin_model, AbstractGumbelAltitudinalModel):
+            return self.model_class_to_estimator_with_finite_aic[StationaryGumbelAltitudinal]
+        else:
+            return self.model_class_to_estimator_with_finite_aic[StationaryAltitudinal]
 
     @property
     def likelihood_ratio(self):
@@ -151,7 +154,9 @@ class OneFoldFit(object):
 
     @property
     def is_significant(self) -> bool:
-        if isinstance(self.best_estimator.margin_model, StationaryAltitudinal):
+        stationary_model_classes = [StationaryAltitudinal, StationaryGumbelAltitudinal]
+        if any([isinstance(self.best_estimator.margin_model, c) 
+                for c in stationary_model_classes]):
             return False
         else:
             return self.likelihood_ratio > chi2.ppf(q=1 - self.SIGNIFICANCE_LEVEL, df=self.degree_freedom_chi2)

test/test_extreme_fit/test_estimator/test_gev_spatio_temporal_extremes_mle.py

@@ -2,12 +2,12 @@ import unittest
 from random import sample
 
 from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day, SafranPrecipitation1Day
-from extreme_fit.model.margin_model.polynomial_margin_model.altitudinal_models import \
+from extreme_fit.model.margin_model.polynomial_margin_model.gev_altitudinal_models import \
     NonStationaryAltitudinalLocationQuadraticScaleLinearCrossTermForLocation, \
     NonStationaryAltitudinalLocationQuadraticCrossTermForLocation, NonStationaryAltitudinalLocationLinear, \
     NonStationaryAltitudinalLocationLinearCrossTermForLocation
-from extreme_fit.model.margin_model.polynomial_margin_model.utils import ALTITUDINAL_MODELS, \
-    MODELS_THAT_SHOULD_RAISE_AN_ASSERTION_ERROR, VARIOUS_SPATIO_TEMPORAL_MODELS
+from extreme_fit.model.margin_model.polynomial_margin_model.utils import ALTITUDINAL_GEV_MODELS, \
+    MODELS_THAT_SHOULD_RAISE_AN_ASSERTION_ERROR, VARIOUS_SPATIO_TEMPORAL_MODELS, ALTITUDINAL_GUMBEL_MODELS
 from extreme_fit.model.margin_model.utils import \
     MarginFitMethod
 from projects.altitude_spatial_model.altitudes_fit.altitudes_studies import AltitudesStudies
@@ -52,8 +52,13 @@ class TestGevTemporalQuadraticExtremesMle(unittest.TestCase):
         for model_class in VARIOUS_SPATIO_TEMPORAL_MODELS:
             self.common_test(model_class)
 
-    def test_altitudinal_models(self):
-        for model_class in ALTITUDINAL_MODELS:
+    def test_altitudinal_gev_models(self):
+        for model_class in ALTITUDINAL_GEV_MODELS:
+            self.common_test(model_class)
+
+    def test_altitudinal_gumbel_models(self):
+        for model_class in ALTITUDINAL_GUMBEL_MODELS[:3]:
+            # print(model_class)
             self.common_test(model_class)