[exceeding] add full diagnositic plot (qqplot & intensity plots) for gumbel...

[exceeding] add full diagnositic plot (qqplot & intensity plots) for gumbel stationary & add gumbel gof tests from r package

extreme_fit/distribution/gumbel/gumbel_gof.R

+# Title     : TODO
+# Objective : TODO
+# Created by: erwan
+# Created on: 15/05/2020
+# library(rmutil)
+library(evd)
+library(gnFit)
+data = rgumbel(10)
+res = gnfit(data, "gum")
+print(res)
+# r = rnorm(0)
+# disp(r)
+# rGumbel(20)

extreme_fit/distribution/gumbel/gumbel_gof.py

+import numpy as np
+
+from extreme_fit.model.result_from_model_fit.abstract_result_from_model_fit import AbstractResultFromModelFit
+from extreme_fit.model.utils import r
+
+
+def cramer_von_mises_and_anderson_darling_tests_pvalues_for_gumbel_distribution(data):
+    res = r.gnfit(data, "gum")
+    res = AbstractResultFromModelFit.get_python_dictionary(res)
+    res = {k: np.array(v)[0] for k, v in res.items()}
+    return res['Wpval'], res['Apval']
+
+
+if __name__ == '__main__':
+    cramer_von_mises_and_anderson_darling_tests_pvalues_for_gumbel_distribution(np.array([2.0, 3.0]))

extreme_fit/model/utils.py

@@ -25,6 +25,7 @@ pandas2ri.activate()
 r.library('SpatialExtremes')
 r.library('data.table')
 r.library('quantreg')
+r.library('gnFit')
 # Desactivate temporarily warnings
 default_filters = warnings.filters.copy()
 warnings.filterwarnings("ignore")

extreme_trend/abstract_gev_trend_test.py

@@ -63,18 +63,17 @@ class AbstractGevTrendTest(object):
     @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
+        test_res = kstest(rvs=quantiles, cdf="gumbel_r")  # 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
+        index_for_significance_level_5_percent = 2
         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]
+        test_res = anderson(quantiles, dist='gumbel_r')  # type: AndersonResult
+        return test_res.statistic < test_res.critical_values[index_for_significance_level_5_percent]
 
     @property
     def degree_freedom_chi2(self) -> int:
@@ -204,7 +203,6 @@ class AbstractGevTrendTest(object):
                     label += '\n' + year_suffix
             else:
                 label = year_suffix
-            print(label)
             self.plot_model(ax, year, end_proba=end_real_proba, label=label, color=color)
 
         ax_lim = [-1.5, 4]

projects/exceeding_snow_loads/checks/qqplot/plot_qqplot.py

@@ -57,6 +57,13 @@ def plot_intensity_against_gumbel_quantile_for_3_examples(
         v.qqplot(m, color)
 
 
+def plot_full_diagnostic(altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends]):
+    for v in altitude_to_visualizer.values():
+        for m in v.massif_name_to_trend_test_that_minimized_aic.keys():
+            v.intensity_plot(m)
+            v.qqplot(m)
+
+
 def plot_qqplot_for_time_series_examples(altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends]):
     marker_altitude_massif_name_for_paper1 = tuples_for_examples_paper1()
     for color, a, m in marker_altitude_massif_name_for_paper1:

projects/exceeding_snow_loads/section_results/main_result_trends_and_return_levels.py

@@ -3,7 +3,7 @@ from multiprocessing.pool import Pool
 import matplotlib as mpl
 
 from projects.exceeding_snow_loads.checks.qqplot.plot_qqplot import \
-    plot_intensity_against_gumbel_quantile_for_3_examples
+    plot_intensity_against_gumbel_quantile_for_3_examples, plot_full_diagnostic
 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
 
@@ -71,9 +71,10 @@ def intermediate_result(altitudes, massif_names=None,
     # 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)
-    plot_intensity_against_gumbel_quantile_for_3_examples(altitude_to_visualizer)
+    # plot_intensity_against_gumbel_quantile_for_3_examples(altitude_to_visualizer)
+    # plot_full_diagnostic(altitude_to_visualizer)
 
 
 def major_result():
@@ -88,8 +89,8 @@ 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, 1800, 2700]
+    # altitudes = paper_altitudes
+    altitudes = [900, 1800, 2700][:1]
     for study_class in study_classes:
         intermediate_result(altitudes, massif_names, model_subsets_for_uncertainty,
                             uncertainty_methods, study_class,

test/test_extreme_fit/test_distribution/test_gumbel.py

+import unittest
+
+from numpy.random.mtrand import gumbel
+
+from extreme_fit.distribution.gumbel.gumbel_gof import \
+    cramer_von_mises_and_anderson_darling_tests_pvalues_for_gumbel_distribution
+from extreme_fit.model.utils import set_seed_for_test
+
+
+class TestGumbel(unittest.TestCase):
+
+    def test_gof_tests(self):
+        set_seed_for_test(seed=42)
+        data = gumbel(size=60)
+        res = cramer_von_mises_and_anderson_darling_tests_pvalues_for_gumbel_distribution(data)
+        cramer_pvalue, anderson_pvalue = res
+        self.assertGreater(cramer_pvalue, 0.05)
+        self.assertGreater(anderson_pvalue, 0.05)
+
+
+if __name__ == '__main__':
+    unittest.main()