[paper 1] add qqplots and standard gumbel transformation

experiment/paper_past_snow_loads/check_mle_convergence_for_trends/qqplot/plot_qqplot.py

+from typing import Dict
+
+from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal
+from experiment.paper_past_snow_loads.data.main_example_swe_total_plot import marker_altitude_massif_name_for_paper1
+from experiment.paper_past_snow_loads.paper_main_utils import load_altitude_to_visualizer
+from experiment.paper_past_snow_loads.study_visualizer_for_non_stationary_trends import \
+    StudyVisualizerForNonStationaryTrends
+from experiment.trend_analysis.univariate_test.extreme_trend_test.abstract_gev_trend_test import AbstractGevTrendTest
+
+
+def plot_qqplot_wrt_standard_gumbel(altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends],
+                                    plot_all=False):
+    if plot_all:
+        pass
+    else:
+        # Plot only some examples
+        plot_qqplot_for_time_series_examples(altitude_to_visualizer)
+        plot_qqplot_for_time_series_with_missing_zeros(altitude_to_visualizer)
+
+
+def plot_qqplot_for_time_series_with_missing_zeros(altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends],
+                                                   nb_worst_examples=3):
+    # Extract all the values
+    l = []
+    for a, v in altitude_to_visualizer.items():
+        l.extend([(a, v, m, p) for m, p in v.massif_name_to_psnow.items()])
+    # Sort them and keep the worst examples
+    l = sorted(l, key=lambda t: t[-1])[:nb_worst_examples]
+    print('Worst examples:')
+    for a, v, m, p in l:
+        print(a, m, p)
+        v.qqplot(m)
+
+
+def plot_qqplot_for_time_series_examples(altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends]):
+    for color, a, m in marker_altitude_massif_name_for_paper1:
+        v = altitude_to_visualizer[a]
+        v.qqplot(m, color)
+
+
+if __name__ == '__main__':
+    # for the five worst, 300 is interesti
+    altitudes = [300, 900, 1800, 2700]
+    altitude_to_visualizer = {altitude:  StudyVisualizerForNonStationaryTrends(CrocusSnowLoadTotal(altitude=altitude),
+                                                                               multiprocessing=True)
+                              for altitude in altitudes}
+    plot_qqplot_wrt_standard_gumbel(altitude_to_visualizer)
+

experiment/paper_past_snow_loads/check_mle_convergence_for_trends/main_shape_repartition.py → experiment/paper_past_snow_loads/check_mle_convergence_for_trends/shape/main_shape_repartition.py

 from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal
-from experiment.paper_past_snow_loads.check_mle_convergence_for_trends.study_visualizer_for_shape_repartition import \
+from experiment.paper_past_snow_loads.check_mle_convergence_for_trends.shape.study_visualizer_for_shape_repartition import \
     StudyVisualizerForShape
 from experiment.paper_past_snow_loads.paper_main_utils import load_altitude_to_visualizer
 

experiment/paper_past_snow_loads/data/main_example_swe_total_plot.py

@@ -7,6 +7,12 @@ from experiment.meteo_france_data.scm_models_data.visualization.study_visualizat
     StudyVisualizer
 from experiment.paper_past_snow_loads.paper_utils import dpi_paper1_figure
 
+marker_altitude_massif_name_for_paper1 = [
+    ('magenta', 900, 'Ubaye'),
+    ('darkmagenta', 1800, 'Vercors'),
+    ('mediumpurple', 2700, 'Beaufortain'),
+]
+
 
 def max_graph_annual_maxima_poster():
     """
@@ -16,15 +22,11 @@ def max_graph_annual_maxima_poster():
     """
     save_to_file = True
     study_class = CrocusSnowLoadTotal
-    marker_altitude_massif_name = [
-        ('magenta', 900, 'Ubaye'),
-        ('darkmagenta', 1800, 'Vercors'),
-        ('mediumpurple', 2700, 'Beaufortain'),
-    ]
+
     ax = plt.gca()
     ax.set_ylim([0, 20])
     ax.set_yticks(list(range(0, 21, 2)))
-    for color, altitude, massif_name in marker_altitude_massif_name:
+    for color, altitude, massif_name in marker_altitude_massif_name_for_paper1:
         for study in study_iterator_global([study_class], altitudes=[altitude]):
             study_visualizer = StudyVisualizer(study, save_to_file=save_to_file,
                                                verbose=True,

experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py

@@ -25,6 +25,8 @@ from experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_two
     GevLocationAgainstGumbel, GevScaleAgainstGumbel
 from experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_two_parameters.gumbel_test_two_parameters import \
     GumbelLocationAndScaleTrendTest
+from extreme_fit.model.margin_model.linear_margin_model.abstract_temporal_linear_margin_model import \
+    TemporalMarginFitMethod
 from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import GumbelTemporalModel
 from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import \
     ConfidenceIntervalMethodFromExtremes
@@ -46,7 +48,7 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
                  effective_temporal_covariate=2017,
                  relative_change_trend_plot=True,
                  non_stationary_trend_test_to_marker=None,
-                 fit_method=None,
+                 fit_method=TemporalMarginFitMethod.extremes_fevd_mle,
                  select_only_acceptable_shape_parameter=False):
         super().__init__(study, show, save_to_file, only_one_graph, only_first_row, vertical_kde_plot,
                          year_for_kde_plot, plot_block_maxima_quantiles, temporal_non_stationarity,
@@ -218,7 +220,9 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
 
     @property
     def label_tdrl_bar(self):
-        return 'Change in {} years'.format(AbstractGevTrendTest.nb_years_for_quantile_evolution)
+        nb_years = AbstractGevTrendTest.nb_years_for_quantile_evolution
+        suffix = 'per decade' if nb_years == 10 else 'in {} years'.format(nb_years)
+        return 'Change {}'.format(suffix)
 
     @property
     def ticks_values_and_labels(self):
@@ -361,3 +365,12 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
                        uncertainty.confidence_interval[1] > eurocode)
                       for eurocode, uncertainty in eurocode_and_uncertainties])
         return 100 * np.mean(a, axis=0)
+
+    # Part 3 - Zeros
+
+    # Part 4 - QQPLOT
+
+    def qqplot(self, massif_name, color=None):
+        trend_test = self.massif_name_to_trend_test_that_minimized_aic[massif_name]
+        marker = self.massif_name_to_marker_style[massif_name]
+        trend_test.qqplot_wrt_standard_gumbel(marker, color)
\ No newline at end of file

experiment/trend_analysis/univariate_test/extreme_trend_test/abstract_gev_trend_test.py

 import numpy as np
+import matplotlib.pyplot as plt
 import pandas as pd
 from cached_property import cached_property
 from scipy.stats import chi2
@@ -161,8 +162,9 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
     def relative_change_in_return_level(self, initial_year, final_year):
         return_level_values = []
         for year in [initial_year, final_year]:
-            gev_params = self.unconstrained_estimator.margin_function_from_fit.get_gev_params(coordinate=np.array([year]),
-                                                                                            is_transformed=False)
+            gev_params = self.unconstrained_estimator.margin_function_from_fit.get_gev_params(
+                coordinate=np.array([year]),
+                is_transformed=False)
             return_level_values.append(gev_params.quantile(self.quantile_level))
         change_until_final_year = self.time_derivative_times_years(nb_years=final_year - initial_year)
         change_in_between = (return_level_values[1] - return_level_values[0])
@@ -204,3 +206,27 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
     @classproperty
     def label(self):
         return '\\mathcal{M}_{%s}'
+
+    # Some display function
+
+    def qqplot_wrt_standard_gumbel(self, marker, color=None):
+        # Standard Gumbel quantiles
+        standard_gumbel_distribution = GevParams(loc=0, scale=1, shape=0)
+        n = len(self.years)
+        standard_gumbel_quantiles = [standard_gumbel_distribution.quantile(i / (n + 1)) for i in range(1, n + 1)]
+        unconstrained_empirical_quantiles = self.compute_empirical_quantiles(self.unconstrained_estimator)
+        constrained_empirical_quantiles = self.compute_empirical_quantiles(self.constrained_estimator)
+        plt.plot(standard_gumbel_quantiles, standard_gumbel_quantiles, color=color)
+        plt.plot(standard_gumbel_quantiles, constrained_empirical_quantiles, 'x')
+        plt.plot(standard_gumbel_quantiles, unconstrained_empirical_quantiles, linestyle='None', **marker)
+        plt.show()
+
+    def compute_empirical_quantiles(self, estimator):
+        empirical_quantiles = []
+        for year, maximum in sorted(zip(self.years, self.maxima), key=lambda t: t[1]):
+            gev_param = estimator.margin_function_from_fit.get_gev_params(
+                coordinate=np.array([year]),
+                is_transformed=False)
+            maximum_standardized = gev_param.gumbel_standardization(maximum)
+            empirical_quantiles.append(maximum_standardized)
+        return empirical_quantiles

extreme_fit/distribution/gev/gev_params.py

@@ -122,4 +122,12 @@ class GevParams(AbstractParams):
             cls.LOC: 'mu',
             cls.SCALE: 'sigma',
             cls.SHAPE: 'zeta',
-        }[gev_param_name]
\ No newline at end of file
+        }[gev_param_name]
+
+    def gumbel_standardization(self, x):
+        x -= self.location
+        x /= self.scale
+        if self.shape == 0:
+            return x
+        else:
+            return np.log(1 + self.shape * x) / self.shape
\ No newline at end of file

extreme_fit/model/margin_model/linear_margin_model/abstract_temporal_linear_margin_model.py

@@ -25,14 +25,16 @@ class AbstractTemporalLinearMarginModel(LinearMarginModel):
     """Linearity only with respect to the temporal coordinates"""
 
     def __init__(self, coordinates: AbstractCoordinates, use_start_value=False, params_start_fit=None,
-                 params_sample=None, starting_point=None, fit_method=TemporalMarginFitMethod.is_mev_gev_fit,
-                 nb_iterations_for_bayesian_fit=5000, params_start_fit_bayesian=None,
+                 params_sample=None, starting_point=None,
+                 fit_method=TemporalMarginFitMethod.is_mev_gev_fit,
+                 nb_iterations_for_bayesian_fit=5000,
+                 params_start_fit_bayesian=None,
                  type_for_MLE="GEV"):
         super().__init__(coordinates, use_start_value, params_start_fit, params_sample, starting_point)
         self.type_for_mle = type_for_MLE
         self.params_start_fit_bayesian = params_start_fit_bayesian
         self.nb_iterations_for_bayesian_fit = nb_iterations_for_bayesian_fit
-        assert isinstance(fit_method, TemporalMarginFitMethod)
+        assert isinstance(fit_method, TemporalMarginFitMethod), fit_method
         self.fit_method = fit_method
 
     def fitmargin_from_maxima_gev(self, data: np.ndarray, df_coordinates_spat: pd.DataFrame,