From 489a0a94d38f196491c37b41b4f4b30ae537f40f Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Fri, 28 Feb 2020 17:58:35 +0100
Subject: [PATCH] [paper 1] add extension of intensity plot with gev params.
add inverse gumbel transformation.
---
.../abstract_gev_trend_test.py | 92 +++++++++++--------
extreme_fit/distribution/gev/gev_params.py | 11 +++
.../qqplot/main_qqplot_for_big_shapes.py | 27 ++++--
...dy_visualizer_for_non_stationary_trends.py | 2 +-
4 files changed, 87 insertions(+), 45 deletions(-)
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 792dcb8e..6111998b 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
@@ -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()
diff --git a/extreme_fit/distribution/gev/gev_params.py b/extreme_fit/distribution/gev/gev_params.py
index df722015..f27beb70 100644
--- a/extreme_fit/distribution/gev/gev_params.py
+++ b/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)
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 bccf350a..f8efbe7a 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
@@ -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)
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 e2b9eb3d..42b87600 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
@@ -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()
--
GitLab