[contrasting] add likelihood & improve plot of the model moments

extreme_fit/estimator/margin_estimator/abstract_margin_estimator.py

@@ -38,7 +38,6 @@ class LinearMarginEstimator(AbstractMarginEstimator):
         return self.dataset.coordinates.df_spatial_coordinates(split=self.train_split,
                                                                drop_duplicates=self.margin_model.drop_duplicates)
 
-
     @property
     def df_coordinates_temp(self):
         return self.dataset.coordinates.df_temporal_coordinates_for_fit(split=self.train_split,
@@ -67,6 +66,9 @@ class LinearMarginEstimator(AbstractMarginEstimator):
             assert not np.isinf(nllh)
         return nllh
 
+    def deviance(self, split=Split.all):
+        return 2 * self.nllh(split=split)
+
     def aic(self, split=Split.all):
         aic = 2 * self.margin_model.nb_params + 2 * self.nllh(split=split)
         npt.assert_almost_equal(self.result_from_model_fit.aic, aic, decimal=5)
@@ -75,4 +77,3 @@ class LinearMarginEstimator(AbstractMarginEstimator):
     def bic(self, split=Split.all):
         n = len(self.dataset.maxima_gev(split=split))
         return np.log(n) * self.margin_model.nb_params + 2 * self.nllh(split=split)
-

projects/altitude_spatial_model/altitudes_fit/altitudes_studies.py

@@ -45,17 +45,25 @@ class AltitudesStudies(object):
     def spatio_temporal_dataset(self, massif_name, s_split_spatial: pd.Series = None,
                                 s_split_temporal: pd.Series = None):
         coordinate_values_to_maxima = {}
-        massif_altitudes = []
-        for altitude in self.altitudes:
+        massif_altitudes = self.massif_name_to_altitudes[massif_name]
+        for altitude in massif_altitudes:
             study = self.altitude_to_study[altitude]
-            if massif_name in study.study_massif_names:
-                massif_altitudes.append(altitude)
-                for year, maxima in zip(study.ordered_years, study.massif_name_to_annual_maxima[massif_name]):
-                    coordinate_values_to_maxima[(altitude, year)] = [maxima]
+            for year, maxima in zip(study.ordered_years, study.massif_name_to_annual_maxima[massif_name]):
+                coordinate_values_to_maxima[(altitude, year)] = [maxima]
         coordinates = self.spatio_temporal_coordinates(s_split_spatial, s_split_temporal, massif_altitudes)
         observations = AnnualMaxima.from_coordinates(coordinates, coordinate_values_to_maxima)
         return AbstractDataset(observations=observations, coordinates=coordinates)
 
+    @cached_property
+    def massif_name_to_altitudes(self):
+        massif_names = self.study.all_massif_names()
+        massif_name_to_altitudes = {massif_name: [] for massif_name in massif_names}
+        for altitude in self.altitudes:
+            study = self.altitude_to_study[altitude]
+            for massif_name in study.study_massif_names:
+                massif_name_to_altitudes[massif_name].append(altitude)
+        return massif_name_to_altitudes
+
     # Coordinates Loader
 
     def spatio_temporal_coordinates(self, s_split_spatial: pd.Series = None, s_split_temporal: pd.Series = None,

projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py

+import matplotlib as mpl
+mpl.rcParams['text.usetex'] = True
+mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
+
 from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day, SafranSnowfall3Days, \
     SafranSnowfall5Days, SafranSnowfall7Days, SafranPrecipitation1Day, SafranPrecipitation3Days, \
     SafranPrecipitation5Days, SafranPrecipitation7Days
@@ -13,8 +17,7 @@ def plot_altitudinal_fit(studies, massif_names=None):
                                                                   model_classes=ALTITUDINAL_MODELS,
                                                                   massif_names=massif_names,
                                                                   show=False)
-    visualizer.plot_mean()
-    visualizer.plot_relative_change()
+    visualizer.plot_moments()
     visualizer.plot_shape_map()
 
 
@@ -29,8 +32,9 @@ def plot_moments(studies, massif_names=None):
 
 
 def main():
-    # 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][4:7]
+    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][:]
@@ -42,7 +46,7 @@ def main():
     for study_class in study_classes:
         studies = AltitudesStudies(study_class, altitudes, season=Season.winter_extended)
         # plot_time_series(studies, massif_names)
-        plot_moments(studies, massif_names)
+        # plot_moments(studies, massif_names)
         plot_altitudinal_fit(studies, massif_names)
 
 

projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py

@@ -34,23 +34,26 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
             old_fold_fit = OneFoldFit(massif_name, dataset, model_classes, self.fit_method)
             self.massif_name_to_one_fold_fit[massif_name] = old_fold_fit
 
-    def plot_mean(self):
-        self.plot_general('mean')
+    def plot_moments(self):
+        for method_name in ['moment', 'changes_in_the_moment', 'relative_changes_in_the_moment']:
+            for order in [1, 2]:
+                self.plot_general(method_name, order)
 
-    def plot_relative_change(self):
-        self.plot_general('relative_changes_in_the_mean')
-
-    def plot_general(self, method_name):
+    def plot_general(self, method_name, order):
         ax = plt.gca()
         min_altitude, *_, max_altitude = self.studies.altitudes
-        altitudes = np.linspace(min_altitude, max_altitude, num=50)
+        altitudes_plot = np.linspace(min_altitude, max_altitude, num=50)
         for massif_id, massif_name in enumerate(self.massif_names):
+            massif_altitudes = self.studies.massif_name_to_altitudes[massif_name]
+            ind = (min(massif_altitudes) <= altitudes_plot) & (altitudes_plot <= max(massif_altitudes))
+            massif_altitudes_plot = altitudes_plot[ind]
             old_fold_fit = self.massif_name_to_one_fold_fit[massif_name]
-            values = old_fold_fit.__getattribute__(method_name)(altitudes)
-            plot_against_altitude(altitudes, ax, massif_id, massif_name, values)
+            values = old_fold_fit.__getattribute__(method_name)(massif_altitudes_plot, order=order)
+            plot_against_altitude(massif_altitudes_plot, ax, massif_id, massif_name, values)
         # Plot settings
         ax.legend(prop={'size': 7}, ncol=3)
         moment = ' '.join(method_name.split('_'))
+        moment = moment.replace('moment', 'mean' if order == 1 else 'std')
         plot_name = '{} annual maxima of {}'.format(moment, SCM_STUDY_CLASS_TO_ABBREVIATION[self.studies.study_class])
         ax.set_ylabel('{} ({})'.format(plot_name, self.study.variable_unit), fontsize=15)
         ax.set_xlabel('altitudes', fontsize=15)

projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py

@@ -2,12 +2,15 @@ import numpy.testing as npt
 import numpy as np
 import rpy2
 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.utils import MarginFitMethod
 
 
 class OneFoldFit(object):
+    SIGNIFICANCE_LEVEL = 0.05
 
     def __init__(self, massif_name, dataset, models_classes, fit_method=MarginFitMethod.extremes_fevd_mle):
         self.massif_name = massif_name
@@ -22,22 +25,37 @@ class OneFoldFit(object):
                                                                                               dataset=self.dataset,
                                                                                               fit_method=self.fit_method)
 
-    def mean(self, altitudes, year=2019):
-        return [self.get_mean(altitude, year) for altitude in altitudes]
-
-    def get_mean(self, altitude, year):
-        return self.get_gev_params(altitude, year).mean
+    def get_moment(self, altitude, year, order=1):
+        gev_params = self.get_gev_params(altitude, year)
+        if order == 1:
+            return gev_params.mean
+        elif order == 2:
+            return gev_params.std
+        else:
+            raise NotImplementedError
 
     def get_gev_params(self, altitude, year):
         coordinate = np.array([altitude, year])
         gev_params = self.best_function_from_fit.get_params(coordinate, is_transformed=False)
         return gev_params
 
-    def relative_changes_in_the_mean(self, altitudes, year=2019, nb_years=50):
+    def moment(self, altitudes, year=2019, order=1):
+        return [self.get_moment(altitude, year, order) for altitude in altitudes]
+
+    def changes_in_the_moment(self, altitudes, year=2019, nb_years=50, order=1):
+        changes = []
+        for altitude in altitudes:
+            mean_after = self.get_moment(altitude, year, order)
+            mean_before = self.get_moment(altitude, year - nb_years, order)
+            change = mean_after - mean_before
+            changes.append(change)
+        return changes
+
+    def relative_changes_in_the_moment(self, altitudes, year=2019, nb_years=50, order=1):
         relative_changes = []
         for altitude in altitudes:
-            mean_after = self.get_mean(altitude, year)
-            mean_before = self.get_mean(altitude, year - nb_years)
+            mean_after = self.get_moment(altitude, year, order)
+            mean_before = self.get_moment(altitude, year - nb_years, order)
             relative_change = 100 * (mean_after - mean_before) / mean_before
             relative_changes.append(relative_change)
         return relative_changes
@@ -75,4 +93,29 @@ class OneFoldFit(object):
 
     @property
     def best_name(self):
-        return self.best_estimator.margin_model.name_str
+        name = self.best_estimator.margin_model.name_str
+        latex_command = 'textbf' if self.is_significant else 'textrm'
+        return '$\\' + latex_command + '{' + name + '}$'
+
+    # Significant
+
+    @property
+    def stationary_estimator(self):
+        for estimator in self.sorted_estimators_with_finite_aic:
+            if isinstance(estimator.margin_model, StationaryAltitudinal):
+                return estimator
+
+    @property
+    def likelihood_ratio(self):
+        return self.stationary_estimator.deviance() - self.best_estimator.deviance()
+
+    @property
+    def degree_freedom_chi2(self):
+        return self.best_estimator.margin_model.nb_params - self.stationary_estimator.margin_model.nb_params
+
+    @property
+    def is_significant(self) -> bool:
+        if isinstance(self.best_estimator.margin_model, StationaryAltitudinal):
+            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

@@ -84,10 +84,11 @@ class TestGevTemporalQuadraticExtremesMle(unittest.TestCase):
 #         self.assertAlmostEqual(estimator.result_from_model_fit.aic, estimator.aic(split=estimator.train_split))
 #         self.assertAlmostEqual(estimator.result_from_model_fit.bic, estimator.bic(split=estimator.train_split))
 #
-#     # def test_altitudinal_models(self):
-#     #     for model_class in ALTITUDINAL_MODELS:
-#     #         self.common_test(model_class)
-#
+#     #
+#     #     # def test_altitudinal_models(self):
+#     #     #     for model_class in ALTITUDINAL_MODELS:
+#     #     #         self.common_test(model_class)
+#     #
 #     def test_wrong(self):
 #         self.common_test(NonStationaryAltitudinalLocationQuadraticCrossTermForLocation)
 #         # self.common_test(NonStationaryAltitudinalLocationQuadraticScaleLinearCrossTermForLocation)