diff --git a/experiment/trend_analysis/univariate_test/extreme_trend_test/abstract_gev_trend_test.py b/experiment/trend_analysis/univariate_test/extreme_trend_test/abstract_gev_trend_test.py
index e3084effa6c0ab5720d3c0a4f4b86fde83729e4f..ff001c28551994496658a3315312ea2ee81dc778 100644
--- a/experiment/trend_analysis/univariate_test/extreme_trend_test/abstract_gev_trend_test.py
+++ b/experiment/trend_analysis/univariate_test/extreme_trend_test/abstract_gev_trend_test.py
@@ -223,9 +223,10 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
epsilon = 0.5
ax_lim = [min(all_quantiles) - epsilon, max(all_quantiles) + epsilon]
ax.plot(standard_gumbel_quantiles, standard_gumbel_quantiles, color='k')
- ax.plot(standard_gumbel_quantiles, constrained_empirical_quantiles, 'x', label='Stationary Gumbel model $\mathcal{M}_0$')
+ ax.plot(standard_gumbel_quantiles, constrained_empirical_quantiles, 'x',
+ label='Stationary Gumbel model $\mathcal{M}_0$')
ax.plot(standard_gumbel_quantiles, unconstrained_empirical_quantiles, linestyle='None',
- label='Selected model $\mathcal{M}_N$', **marker)
+ label='Selected model $\mathcal{M}_N$', **marker)
ax.set_xlabel("Standard Gumbel quantile", fontsize=size)
ax.set_ylabel("Standard Empirical quantile", fontsize=size)
ax.legend(loc='upper left', prop={'size': 10})
@@ -239,6 +240,19 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
plt.show()
+ def return_level_plot_comparison(self, ax, label, color=None):
+ # ax = plt.gca()
+ size = 15
+ # Load Gev parameter in 2017 for the unconstrained estimator
+ gev_params = self.unconstrained_estimator.margin_function_from_fit.get_gev_params(coordinate=np.array([2017]),
+ is_transformed=False) # type: GevParams
+ gev_params_with_corrected_shape = GevParams(loc=gev_params.location,
+ scale=gev_params.scale,
+ shape=0.5)
+ suffix = 'in 2017'
+ gev_params.return_level_plot_against_return_period(ax, color, linestyle='-', label=label, suffix_return_level_label=suffix)
+ gev_params_with_corrected_shape.return_level_plot_against_return_period(ax, color=color, linestyle='--', suffix_return_level_label=suffix)
+
def compute_empirical_quantiles(self, estimator):
empirical_quantiles = []
for year, maximum in sorted(zip(self.years, self.maxima), key=lambda t: t[1]):
diff --git a/extreme_fit/distribution/gev/gev_params.py b/extreme_fit/distribution/gev/gev_params.py
index eb957c3564ba25d5c70edbe2261910caa5a62cf0..f5bab17e1a9d91e8d47e38a3b4b223e8d8bf71eb 100644
--- a/extreme_fit/distribution/gev/gev_params.py
+++ b/extreme_fit/distribution/gev/gev_params.py
@@ -1,4 +1,5 @@
from collections import OrderedDict
+import matplotlib.pyplot as plt
from typing import List
from cached_property import cached_property
@@ -148,4 +149,23 @@ class GevParams(AbstractParams):
if self.shape <= 0:
return np.inf
else:
- return self.bound
\ No newline at end of file
+ return self.bound
+
+ def return_level_plot_against_return_period(self, ax=None, color=None, linestyle=None, label=None, show=False,
+ suffix_return_level_label=''):
+ if ax is None:
+ ax = plt.gca()
+ # Plot return level against return period
+ return_periods = list(range(2, 61))
+ quantiles = [self.quantile(1 - 1 / return_period) for return_period in return_periods]
+ return_period_to_quantile = dict(zip(return_periods, quantiles))
+ ax.vlines(50, 0, return_period_to_quantile[50])
+ ax.plot(return_periods, quantiles, color=color, linestyle=linestyle, label=label)
+ ax.set_xlabel('Return period')
+ ax.legend()
+
+ ax.set_xticks([10 * i for i in range(1, 7)])
+ ax.set_ylabel('Return level {}'.format(suffix_return_level_label))
+ plt.gca().set_ylim(bottom=0)
+ if show:
+ plt.show()
diff --git a/papers/exceeding_snow_loads/check_mle_convergence_for_trends/qqplot/main_qqplot_for_big_shapes.py b/papers/exceeding_snow_loads/check_mle_convergence_for_trends/qqplot/main_qqplot_for_big_shapes.py
index 330cef7aaef0096623a20af1ea5384086f9b8951..291851f2d56f1a29a7ddcb8bf8aaf20c9c698796 100644
--- a/papers/exceeding_snow_loads/check_mle_convergence_for_trends/qqplot/main_qqplot_for_big_shapes.py
+++ b/papers/exceeding_snow_loads/check_mle_convergence_for_trends/qqplot/main_qqplot_for_big_shapes.py
@@ -1,30 +1,35 @@
from typing import Dict
+import matplotlib.pyplot as plt
from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal
from experiment.meteo_france_data.scm_models_data.visualization.study_visualization.main_study_visualizer import \
ALL_ALTITUDES_WITHOUT_NAN
-from experiment.exceeding_snow_loads.study_visualizer_for_non_stationary_trends import \
+from papers.exceeding_snow_loads.study_visualizer_for_non_stationary_trends import \
StudyVisualizerForNonStationaryTrends
-def qqplots_for_biggest_shape_parameters_before_selection():
- altitudes = ALL_ALTITUDES_WITHOUT_NAN
+
+def get_tuple_ordered_by_shape(fast=False):
+ if fast:
+ altitudes = [300]
+ else:
+ altitudes = ALL_ALTITUDES_WITHOUT_NAN
altitude_to_visualizer = {altitude: StudyVisualizerForNonStationaryTrends(CrocusSnowLoadTotal(altitude=altitude),
select_only_acceptable_shape_parameter=False,
multiprocessing=True)
for altitude in altitudes}
- plot_qqplot_for_time_series_with_worst_shape_parameters(altitude_to_visualizer, nb_worst_examples=5)
-
-
-def plot_qqplot_for_time_series_with_worst_shape_parameters(
- 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, t.unconstrained_estimator_gev_params.shape) for m, t in v.massif_name_to_trend_test_that_minimized_aic.items()])
+ l.extend([(a, v, m, t.unconstrained_estimator_gev_params.shape) for m, t in
+ v.massif_name_to_trend_test_that_minimized_aic.items()])
# Sort them and keep the highest examples
l = sorted(l, key=lambda t: t[-1])
+ return l
+
+
+def plot_qqplot_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)
@@ -34,6 +39,19 @@ def plot_qqplot_for_time_series_with_worst_shape_parameters(
print(a, m, shape)
v.qqplot(m)
+
+def plot_return_level_for_time_series_with_big_shape_parameters(tuple_ordered_by_shape, nb_worst_examples=5):
+ # Extract all the values
+ l = tuple_ordered_by_shape
+ print('Highest examples:')
+ ax = plt.gca()
+ colors = ['orange', 'red', 'blue', 'green', 'yellow']
+ for (a, v, m, shape), color in zip(l[-nb_worst_examples:][::-1], colors):
+ print(a, m, shape, color)
+ v.return_level_plot(ax, m, color)
+ plt.show()
+
+
"""
10 Worst examples:
300 Oisans 1.1070477650581747
@@ -50,7 +68,6 @@ def plot_qqplot_for_time_series_with_worst_shape_parameters(
300 Bauges 0.39291367817905504
"""
-
"""
for the worst example for -shape
@@ -66,6 +83,8 @@ for the worst example for -shape
2700 Mont-Blanc 0.2712596135797868
"""
-
if __name__ == '__main__':
- qqplots_for_biggest_shape_parameters_before_selection()
+ fast = False
+ 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)
diff --git a/papers/exceeding_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py b/papers/exceeding_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py
index 043da14ed538264a5eef132035deee569cd24551..d0a7ab8d9070229f96d32b8dd3d582100ad1684b 100644
--- a/papers/exceeding_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py
+++ b/papers/exceeding_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py
@@ -2,7 +2,7 @@ from multiprocessing.pool import Pool
import matplotlib as mpl
-from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal
+from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal, CrocusSnowLoad3Days
from papers.exceeding_snow_loads.paper_main_utils import load_altitude_to_visualizer
from papers.exceeding_snow_loads.paper_utils import paper_study_classes, paper_altitudes
from papers.exceeding_snow_loads.result_trends_and_return_levels.plot_diagnosis_risk import plot_diagnosis_risk
diff --git a/papers/exceeding_snow_loads/study_visualizer_for_non_stationary_trends.py b/papers/exceeding_snow_loads/study_visualizer_for_non_stationary_trends.py
index cd7b174888a0892da556ca242bbddd6bf6bcc05e..df7703e97fc90d09d9b716096fb306acbc7ea523 100644
--- a/papers/exceeding_snow_loads/study_visualizer_for_non_stationary_trends.py
+++ b/papers/exceeding_snow_loads/study_visualizer_for_non_stationary_trends.py
@@ -401,6 +401,12 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
marker = self.massif_name_to_marker_style[massif_name]
trend_test.qqplot_wrt_standard_gumbel(marker, color)
+ def return_level_plot(self, ax, massif_name, color=None):
+ trend_test = self.massif_name_to_trend_test_that_minimized_aic[massif_name]
+ model_name = '${}$'.format(trend_test.label)
+ label = '{} at {}m with {}'.format(massif_name, self.altitude, model_name)
+ trend_test.return_level_plot_comparison(ax, label, color)
+
# Part 4 - Trend plot
@property
diff --git a/test/test_extreme_fit/test_distribution/test_gev/test_gev_params.py b/test/test_extreme_fit/test_distribution/test_gev/test_gev_params.py
index 41e9d2f24b10eb177fd846a9f8ce94e30685fa61..c20b5782359dd03f49baac6903e3831629806c7b 100644
--- a/test/test_extreme_fit/test_distribution/test_gev/test_gev_params.py
+++ b/test/test_extreme_fit/test_distribution/test_gev/test_gev_params.py
@@ -54,6 +54,10 @@ class TestGevParams(unittest.TestCase):
gev_params = GevParams(loc=mu, shape=chi, scale=sigma)
self.assertEqual(gev_params.variance, ((sigma / chi) ** 2) * (gamma(1 - 2 * chi) - (gamma(1 - chi) ** 2)))
+ def test_return_level_plot(self):
+ gev_params = GevParams(loc=0.0, shape=0.0, scale=1.0)
+ gev_params.return_level_plot_against_return_period(show=False)
+
if __name__ == '__main__':
unittest.main()