[paper 1] add extension of intensity plot with gev params. add inverse gumbel transformation.

experiment/trend_analysis/univariate_test/extreme_trend_test/abstract_gev_trend_test.py

@@ -228,7 +228,7 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
         print(unconstrained_empirical_quantiles)
 
         ax_twiny = ax.twiny()
-        return_periods = [10, 25, 50, 100]
+        return_periods = [10, 25, 50]
         quantiles = self.get_standard_quantiles_for_return_periods(return_periods, psnow)
         print(quantiles)
         ax_twiny.plot(quantiles, [0 for _ in quantiles], linewidth=0)
@@ -237,42 +237,28 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
         ax_twiny.set_xlabel('Return period w.r.t all annual maxima of GSL (years)', fontsize=size)
 
         # Plot for the discarded model
-        if 'Verdon' in massif_name and altitude == 300:
-            q = [-1.4541688117485054, -1.2811308174310914, -1.216589300814509, -0.7635793791201918, -0.6298883422064275,
-                 -0.5275954855697504, -0.4577268043676126, -0.4497570331795861, -0.1647955002136654,
-                 -0.14492222503785876, -0.139173823298689, -0.11945617994263039, -0.07303100174657867,
-                 -5.497308509286266e-05, 0.13906416388625908, 0.15274793441408543, 0.1717763342727519,
-                 0.17712605315013535, 0.17900143646245203, 0.371986176207554, 0.51640780422156, 0.7380550963951035,
-                 0.7783015252180445, 0.887836077295502, 0.917853338231094, 0.9832396811506262, 1.0359396416309927,
-                 1.1892663813729711, 1.2053261113817888, 1.5695111391491652, 2.3223652143938476, 2.674882764437432,
-                 2.6955728524900406, 2.8155882785356896, 3.282838470153471, 3.2885313947906765]
-            color = 'red'
-            ax.plot(q, sorted_maxima,
-                    label='Discarded model, which is ${}$\n'.format('\mathcal{M}_{\zeta_0, \sigma_1}')
-                          + 'with $\zeta_0=0.84$',
-                    color=color)
-
-            # Extend the curve linear a bit if the return period 50 is not in the quantiles
-
-            def compute_slope_intercept(x, y):
-                x1, x2 = x[-2:]
-                y1, y2 = y[-2:]
-                a = (y2 - y1) / (x2 - x1)
-                b = y1 - a * x1
-                return a, b
-
-            def compute_maxima_corresponding_to_return_period(return_period_quantiles, quantiles, model_maxima):
-                a, b = compute_slope_intercept(quantiles, model_maxima)
-                return a * return_period_quantiles + b
-
-            quantile_return_period_50 = quantiles[-1]
-            if max(q) < quantile_return_period_50:
-                maxima_extended = compute_maxima_corresponding_to_return_period(quantile_return_period_50,
-                                                                                q,
-                                                                                sorted_maxima)
-                ax.plot([q[-1], quantile_return_period_50],
-                        [sorted_maxima[-1], maxima_extended], linestyle='--',
-                        color=color)
+        # if 'Verdon' in massif_name and altitude == 300:
+        #     q = [-1.4541688117485054, -1.2811308174310914, -1.216589300814509, -0.7635793791201918, -0.6298883422064275,
+        #          -0.5275954855697504, -0.4577268043676126, -0.4497570331795861, -0.1647955002136654,
+        #          -0.14492222503785876, -0.139173823298689, -0.11945617994263039, -0.07303100174657867,
+        #          -5.497308509286266e-05, 0.13906416388625908, 0.15274793441408543, 0.1717763342727519,
+        #          0.17712605315013535, 0.17900143646245203, 0.371986176207554, 0.51640780422156, 0.7380550963951035,
+        #          0.7783015252180445, 0.887836077295502, 0.917853338231094, 0.9832396811506262, 1.0359396416309927,
+        #          1.1892663813729711, 1.2053261113817888, 1.5695111391491652, 2.3223652143938476, 2.674882764437432,
+        #          2.6955728524900406, 2.8155882785356896, 3.282838470153471, 3.2885313947906765]
+        #     color = 'red'
+        #     ax.plot(q, sorted_maxima,
+        #             label='Discarded model, which is ${}$\n'.format('\mathcal{M}_{\zeta_0, \sigma_1}')
+        #                   + 'with $\zeta_0=0.84$',
+        #             color=color)
+
+        if self.unconstrained_estimator_gev_params.shape > 0.5:
+            # self.linear_extension(ax, unconstrained_empirical_quantiles, quantiles, sorted_maxima)
+            max_model_quantile = max(unconstrained_empirical_quantiles)
+            # print(sorted(unconstrained_empirical_quantiles))
+            # self.gev_param_extension(ax, 1958, max_model_quantile)
+            self.gev_param_extension(ax, 2017, max_model_quantile)
+
 
         all_quantiles = standard_gumbel_quantiles + unconstrained_empirical_quantiles + quantiles
         epsilon = 0.5
@@ -281,11 +267,41 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
         ax_twiny.set_xlim(ax_lim)
 
         ax.legend()
+        ax.yaxis.grid()
 
         ax.set_xlabel("Standard Gumbel quantile", fontsize=size)
         ax.set_ylabel("Non-zero annual maxima of GSL ({})".format(AbstractSnowLoadVariable.UNIT), fontsize=size)
-        ax.legend(loc='lower right', prop={'size': 10})
+        ax.legend(loc='upper left', prop={'size': 10})
 
+    def gev_param_extension(self, ax, year, max_model_quantile):
+        quantile_standard_gumbel = GevParams(0, 1, 0).quantile(0.98)
+        extended_quantiles = np.linspace(max_model_quantile, quantile_standard_gumbel, 50)
+        gev_params_year = self.unconstrained_estimator.margin_function_from_fit.get_gev_params(
+                coordinate=np.array([year]),
+                is_transformed=False)
+        extended_maxima = [gev_params_year.gumbel_inverse_standardization(q) for q in extended_quantiles]
+        ax.plot(extended_quantiles, extended_maxima, linestyle='--', label='{} extension'.format(year))
+
+    def linear_extension(self, ax, q, quantiles, sorted_maxima):
+        # Extend the curve linear a bit if the return period 50 is not in the quantiles
+        def compute_slope_intercept(x, y):
+            x1, x2 = x[-2:]
+            y1, y2 = y[-2:]
+            a = (y2 - y1) / (x2 - x1)
+            b = y1 - a * x1
+            return a, b
+
+        def compute_maxima_corresponding_to_return_period(return_period_quantiles, quantiles, model_maxima):
+            a, b = compute_slope_intercept(quantiles, model_maxima)
+            return a * return_period_quantiles + b
+
+        quantile_return_period_50 = quantiles[-1]
+        if max(q) < quantile_return_period_50:
+            maxima_extended = compute_maxima_corresponding_to_return_period(quantile_return_period_50,
+                                                                            q,
+                                                                            sorted_maxima)
+            ax.plot([q[-1], quantile_return_period_50],
+                    [sorted_maxima[-1], maxima_extended], linestyle='--', label='linear extension')
 
     def qqplot_wrt_standard_gumbel(self, massif_name, altitude):
         ax = plt.gca()

extreme_fit/distribution/gev/gev_params.py

@@ -133,6 +133,17 @@ class GevParams(AbstractParams):
         else:
             return np.log(1 + self.shape * x) / self.shape
 
+    def gumbel_inverse_standardization(self, x):
+        if self.shape == 0:
+            x = x
+        else:
+            x = (np.exp(self.shape * x) - 1) / self.shape
+        x *= self.scale
+        x += self.location
+        return x
+
+
+
     @property
     def bound(self):
         return self.location - (self.scale / self.shape)

papers/exceeding_snow_loads/check_mle_convergence_for_trends/qqplot/main_qqplot_for_big_shapes.py

@@ -8,7 +8,6 @@ from papers.exceeding_snow_loads.study_visualizer_for_non_stationary_trends impo
     StudyVisualizerForNonStationaryTrends
 
 
-
 def get_tuple_ordered_by_shape(fast=False):
     if fast:
         altitudes = [300]
@@ -16,7 +15,9 @@ def get_tuple_ordered_by_shape(fast=False):
         altitudes = ALL_ALTITUDES_WITHOUT_NAN
     altitude_to_visualizer = {altitude: StudyVisualizerForNonStationaryTrends(CrocusSnowLoadTotal(altitude=altitude),
                                                                               select_only_acceptable_shape_parameter=False,
-                                                                              multiprocessing=True)
+                                                                              multiprocessing=True,
+                                                                              save_to_file=True,
+                                                                              show=False)
                               for altitude in altitudes}
     # Extract all the values
     l = []
@@ -35,9 +36,22 @@ def plot_qqplot_for_time_series_with_worst_shape_parameters(tuple_ordered_by_sha
         print(a, m, shape)
         v.qqplot(m)
     print('Lowest examples:')
-    for a, v, m, shape in l[:1]:
+    for a, v, m, shape in l[:5]:
         print(a, m, shape)
-        v.qqplot(m)
+        # v.qqplot(m)
+
+
+def plot_intensity_for_time_series_with_worst_shape_parameters(tuple_ordered_by_shape, nb_worst_examples=5):
+    l = tuple_ordered_by_shape
+    print('Highest examples:')
+    for a, v, m, shape in l[-nb_worst_examples:][::-1]:
+        print(a, m, shape)
+        v.intensity_plot(m, v.massif_name_to_psnow[m])
+    print('Lowest examples:')
+    for a, v, m, shape in l[:5]:
+        print(a, m, shape)
+        # v.qqplot(m)
+    #     v.intensity_plot(m, v.massif_name_to_psnow[m])
 
 
 def plot_return_level_for_time_series_with_big_shape_parameters(tuple_ordered_by_shape, nb_worst_examples=5):
@@ -85,7 +99,8 @@ for the worst example for -shape
 """
 
 if __name__ == '__main__':
-    fast = False
+    fast = True
     nb = 1 if fast else 5
     tuple_ordered_by_shape = get_tuple_ordered_by_shape(fast=fast)
-    plot_return_level_for_time_series_with_big_shape_parameters(tuple_ordered_by_shape, nb_worst_examples=nb)
+    # plot_return_level_for_time_series_with_big_shape_parameters(tuple_ordered_by_shape, nb_worst_examples=nb)
+    plot_intensity_for_time_series_with_worst_shape_parameters(tuple_ordered_by_shape, nb_worst_examples=nb)

papers/exceeding_snow_loads/study_visualizer_for_non_stationary_trends.py

@@ -399,7 +399,7 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
     def intensity_plot(self, massif_name, psnow, color=None):
         trend_test = self.massif_name_to_trend_test_that_minimized_aic[massif_name]
         trend_test.intensity_plot_wrt_standard_gumbel(massif_name, self.altitude, psnow)
-        self.plot_name = 'intensity_plot_{}_{}'.format(self.altitude, psnow)
+        self.plot_name = 'intensity_plot_{}_{}_{}_{}'.format(self.altitude, massif_name, psnow, trend_test.unconstrained_estimator_gev_params.shape)
         self.show_or_save_to_file(add_classic_title=False, no_title=True)
         plt.close()