[exceeding] add anderson and KS test for the goodness of fit

extreme_trend/abstract_gev_trend_test.py

@@ -3,7 +3,9 @@ from math import ceil, floor
 import matplotlib.pyplot as plt
 import numpy as np
 from cached_property import cached_property
-from scipy.stats import chi2
+from scipy.stats import chi2, kstest, anderson
+from scipy.stats.morestats import AndersonResult
+from scipy.stats.stats import KstestResult
 
 from extreme_data.eurocode_data.utils import EUROCODE_QUANTILE, YEAR_OF_INTEREST_FOR_RETURN_LEVEL
 from extreme_data.meteo_france_data.scm_models_data.crocus.crocus_variables import AbstractSnowLoadVariable
@@ -54,6 +56,22 @@ class AbstractGevTrendTest(object):
     def is_significant(self) -> bool:
         return self.likelihood_ratio > chi2.ppf(q=1 - self.SIGNIFICANCE_LEVEL, df=self.degree_freedom_chi2)
 
+    @property
+    def goodness_of_fit_ks_test(self):
+        quantiles = self.compute_empirical_quantiles(estimator=self.unconstrained_estimator)
+        test_res = kstest(rvs=quantiles, cdf="gumbel_l")  # type: KstestResult
+        return test_res.pvalue < self.SIGNIFICANCE_LEVEL
+
+    @property
+    def goodness_of_fit_anderson_test(self):
+        assert self.SIGNIFICANCE_LEVEL == 0.05
+        # significance_level=array([25. , 10. ,  5. ,  2.5,  1. ]))
+        index_for_significance_level_5_percent = 4
+        quantiles = self.compute_empirical_quantiles(estimator=self.unconstrained_estimator)
+        test_res = anderson(quantiles, dist='gumbel')  # type: AndersonResult
+        print(test_res)
+        return test_res.statistic > test_res.critical_values[index_for_significance_level_5_percent]
+
     @property
     def degree_freedom_chi2(self) -> int:
         raise NotImplementedError

extreme_trend/visualizers/study_visualizer_for_non_stationary_trends.py

@@ -305,6 +305,18 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
         return Counter(
             [type(t) for t in self.massif_name_to_trend_test_that_minimized_aic.values() if t.is_significant])
 
+    @cached_property
+    def selected_and_anderson_goodness_of_fit_trend_test_class_counter(self):
+        return Counter(
+            [type(t) for t in self.massif_name_to_trend_test_that_minimized_aic.values()
+             if t.goodness_of_fit_anderson_test])
+
+    @cached_property
+    def selected_and_kstest_goodness_of_fit_trend_test_class_counter(self):
+        return Counter(
+            [type(t) for t in self.massif_name_to_trend_test_that_minimized_aic.values()
+             if t.goodness_of_fit_ks_test])
+
     @cached_property
     def massif_name_to_marker_style(self):
         d = {}

projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/main_snowfall_article.py

@@ -74,8 +74,8 @@ def intermediate_result(altitudes, massif_names=None,
     # Plots
     validation_plot(altitude_to_visualizer, order_derivative=0)
     validation_plot(altitude_to_visualizer, order_derivative=1)
-    # plot_snowfall_mean(altitude_to_visualizer)
-    # plot_snowfall_time_derivative_mean(altitude_to_visualizer)
+    plot_snowfall_mean(altitude_to_visualizer)
+    plot_snowfall_time_derivative_mean(altitude_to_visualizer)
 
 
 def major_result():

projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/validation_plot.py

@@ -24,20 +24,23 @@ def validation_plot(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanValu
         altitude_to_relative_differences[altitude] = plot_function(visualizer)
         study_visualizer.show_or_save_to_file(add_classic_title=False, dpi=500)
     # Shoe plot with respect to the altitude.
-    plot_shoe_relative_differences_distribution(altitude_to_relative_differences, altitudes, study_visualizer, order_derivative)
+    plot_shoe_relative_differences_distribution(altitude_to_relative_differences, altitudes, study_visualizer,
+                                                order_derivative)
     study_visualizer.show_or_save_to_file(add_classic_title=False, dpi=500)
     plt.close()
 
 
-def plot_shoe_relative_differences_distribution(altitude_to_relative_differences, altitudes, visualizer, order_derivative):
+def plot_shoe_relative_differences_distribution(altitude_to_relative_differences, altitudes, visualizer,
+                                                order_derivative):
     study = visualizer.study
     ax = plt.gca()
     width = 150
     ax.boxplot([altitude_to_relative_differences[a] for a in altitudes], positions=altitudes, widths=width)
     ax.set_xlim([min(altitudes) - width, max(altitudes) + width])
     moment = '' if order_derivative == 0 else 'time derivative of '
-    ylabel = 'Relative difference of the {} model mean \n' \
-             'w.r.t. the {}empirical mean of {} (\%)'.format(moment, moment, SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)])
+    ylabel = 'Global relative difference of the {} model mean \n' \
+             'w.r.t. the {}empirical mean of {} (\%)'.format(moment, moment,
+                                                             SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)])
     visualizer.plot_trends = ylabel
     ax.set_ylabel(ylabel)
     ax.set_xlabel('Altitude (m)')
@@ -57,14 +60,17 @@ def plot_relative_difference_map_order_zero(visualizer: StudyVisualizerForMeanVa
                                         'for the ' + label, graduation=1)
     return list(visualizer.massif_name_to_relative_difference_for_mean.values())
 
+
 def plot_relative_difference_map_order_one(visualizer: StudyVisualizerForMeanValues):
     study = visualizer.study
-    label = ' time derivative of mean annual maxima of {} ({})'.format(SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)], study.variable_unit)
+    label = ' time derivative of mean annual maxima of {} ({})'.format(SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)],
+                                                                       study.variable_unit)
     visualizer.plot_abstract_fast(massif_name_to_value=visualizer.massif_name_to_change_ratio_in_empirical_mean,
                                   label='Empirical' + label, negative_and_positive_values=False, graduation=0.5)
     visualizer.plot_abstract_fast(massif_name_to_value=visualizer.massif_name_to_change_ratio_in_model_mean,
                                   label='Model' + label, negative_and_positive_values=False, graduation=0.5)
-    visualizer.plot_abstract_fast(massif_name_to_value=visualizer.massif_name_to_relative_difference_for_change_ratio_in_mean,
-                                  label='Relative difference of the model mean w.r.t. the empirical mean \n'
-                                        'for the ' + label, graduation=1)
-    return list(visualizer.massif_name_to_relative_difference_for_mean.values())
+    visualizer.plot_abstract_fast(
+        massif_name_to_value=visualizer.massif_name_to_relative_difference_for_change_ratio_in_mean,
+        label='Relative difference of the model mean w.r.t. the empirical mean \n'
+              'for the ' + label, graduation=5)
\ No newline at end of file
+    return list(visualizer.massif_name_to_relative_difference_for_change_ratio_in_mean.values())

projects/exceeding_snow_loads/section_results/main_result_trends_and_return_levels.py

@@ -2,6 +2,7 @@ from multiprocessing.pool import Pool
 
 import matplotlib as mpl
 
+from projects.exceeding_snow_loads.section_results.plot_selection_curves import plot_selection_curves
 from projects.exceeding_snow_loads.section_results.plot_trend_curves import plot_trend_map
 
 mpl.rcParams['text.usetex'] = True
@@ -64,11 +65,11 @@ def intermediate_result(altitudes, massif_names=None,
             _ = compute_minimized_aic(visualizer)
 
     # Plots
-    plot_trend_map(altitude_to_visualizer)
+    # plot_trend_map(altitude_to_visualizer)
     # plot_trend_curves(altitude_to_visualizer={a: v for a, v in altitude_to_visualizer.items() if a >= 900})
     # plot_uncertainty_massifs(altitude_to_visualizer)
     # plot_uncertainty_histogram(altitude_to_visualizer)
-    # plot_selection_curves(altitude_to_visualizer)
+    plot_selection_curves(altitude_to_visualizer)
     # uncertainty_interval_size(altitude_to_visualizer)
 
 
@@ -84,9 +85,12 @@ def major_result():
     #                                  ModelSubsetForUncertainty.non_stationary_gumbel_and_gev]
     model_subsets_for_uncertainty = None
     # study_classes = [CrocusSnowLoad3Days, CrocusSnowLoad5Days, CrocusSnowLoad7Days][::-1]
+    altitudes = paper_altitudes
+    # altitudes = [900]
     for study_class in study_classes:
-        intermediate_result(paper_altitudes, massif_names, model_subsets_for_uncertainty,
-                            uncertainty_methods, study_class)
+        intermediate_result(altitudes, massif_names, model_subsets_for_uncertainty,
+                            uncertainty_methods, study_class,
+                            multiprocessing=True)
 
 
 if __name__ == '__main__':

projects/exceeding_snow_loads/section_results/plot_selection_curves.py

@@ -21,6 +21,8 @@ def plot_selection_curves(altitude_to_visualizer: Dict[int, StudyVisualizerForNo
 
     selected_counter = merge_counter([v.selected_trend_test_class_counter for v in altitude_to_visualizer.values()])
     selected_and_significative_counter = merge_counter([v.selected_and_significative_trend_test_class_counter for v in altitude_to_visualizer.values()])
+    # selected_and_significative_counter = merge_counter([v.selected_and_anderson_goodness_of_fit_trend_test_class_counter for v in altitude_to_visualizer.values()])
+    # selected_and_significative_counter = merge_counter([v.selected_and_kstest_goodness_of_fit_trend_test_class_counter for v in altitude_to_visualizer.values()])
     total_of_selected_models = sum(selected_counter.values())
     l = sorted(enumerate(visualizer.non_stationary_trend_test), key=lambda e: selected_counter[e[1]])
     permutation = [i for i, v in l][::-1]