From c299e2b0b30becd6ae0a54bc9b41623a33a7e937 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Fri, 15 May 2020 11:34:58 +0200
Subject: [PATCH] [exceeding] add full diagnositic plot (qqplot & intensity
plots) for gumbel stationary & add gumbel gof tests from r package
---
extreme_fit/distribution/gumbel/__init__.py | 0
extreme_fit/distribution/gumbel/gumbel_gof.R | 13 +++++++++++
extreme_fit/distribution/gumbel/gumbel_gof.py | 15 +++++++++++++
extreme_fit/model/utils.py | 1 +
extreme_trend/abstract_gev_trend_test.py | 12 +++++-----
.../checks/qqplot/plot_qqplot.py | 7 ++++++
.../main_result_trends_and_return_levels.py | 11 +++++-----
.../test_distribution/test_gumbel.py | 22 +++++++++++++++++++
8 files changed, 69 insertions(+), 12 deletions(-)
create mode 100644 extreme_fit/distribution/gumbel/__init__.py
create mode 100644 extreme_fit/distribution/gumbel/gumbel_gof.R
create mode 100644 extreme_fit/distribution/gumbel/gumbel_gof.py
create mode 100644 test/test_extreme_fit/test_distribution/test_gumbel.py
diff --git a/extreme_fit/distribution/gumbel/__init__.py b/extreme_fit/distribution/gumbel/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/extreme_fit/distribution/gumbel/gumbel_gof.R b/extreme_fit/distribution/gumbel/gumbel_gof.R
new file mode 100644
index 00000000..847a8246
--- /dev/null
+++ b/extreme_fit/distribution/gumbel/gumbel_gof.R
@@ -0,0 +1,13 @@
+# 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)
diff --git a/extreme_fit/distribution/gumbel/gumbel_gof.py b/extreme_fit/distribution/gumbel/gumbel_gof.py
new file mode 100644
index 00000000..39aded04
--- /dev/null
+++ b/extreme_fit/distribution/gumbel/gumbel_gof.py
@@ -0,0 +1,15 @@
+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]))
diff --git a/extreme_fit/model/utils.py b/extreme_fit/model/utils.py
index c04eb802..968b877d 100644
--- a/extreme_fit/model/utils.py
+++ b/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")
diff --git a/extreme_trend/abstract_gev_trend_test.py b/extreme_trend/abstract_gev_trend_test.py
index 66b7b889..2f79e7b3 100644
--- a/extreme_trend/abstract_gev_trend_test.py
+++ b/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]
diff --git a/projects/exceeding_snow_loads/checks/qqplot/plot_qqplot.py b/projects/exceeding_snow_loads/checks/qqplot/plot_qqplot.py
index 52951672..762f45f5 100644
--- a/projects/exceeding_snow_loads/checks/qqplot/plot_qqplot.py
+++ b/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:
diff --git a/projects/exceeding_snow_loads/section_results/main_result_trends_and_return_levels.py b/projects/exceeding_snow_loads/section_results/main_result_trends_and_return_levels.py
index 6e06ab37..c22c2832 100644
--- a/projects/exceeding_snow_loads/section_results/main_result_trends_and_return_levels.py
+++ b/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,
diff --git a/test/test_extreme_fit/test_distribution/test_gumbel.py b/test/test_extreme_fit/test_distribution/test_gumbel.py
new file mode 100644
index 00000000..a1a579ed
--- /dev/null
+++ b/test/test_extreme_fit/test_distribution/test_gumbel.py
@@ -0,0 +1,22 @@
+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()
--
GitLab