From 535d490c2b810366427bfd162ffcc81c3ef546e9 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Fri, 28 Feb 2020 18:39:07 +0100
Subject: [PATCH] [paper 1] add extension of intensity plot with gev params.
add inverse gumbel transformation.
---
.../abstract_gev_trend_test.py | 37 +++++++++++--------
.../qqplot/main_qqplot_for_big_shapes.py | 2 +-
2 files changed, 22 insertions(+), 17 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 6111998b..295336f5 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
@@ -221,21 +221,26 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
ax.plot(standard_gumbel_quantiles, sorted_maxima, linestyle='None',
label='Empirical model', marker='o')
- # Plot for the selected model
- unconstrained_empirical_quantiles = self.compute_empirical_quantiles(self.unconstrained_estimator)
- ax.plot(unconstrained_empirical_quantiles, sorted_maxima,
- label='Selected model, which is ${}$'.format(self.label))
- print(unconstrained_empirical_quantiles)
+
ax_twiny = ax.twiny()
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)
ax_twiny.set_xticks(quantiles)
ax_twiny.set_xticklabels([return_period for return_period in return_periods])
ax_twiny.set_xlabel('Return period w.r.t all annual maxima of GSL (years)', fontsize=size)
+ # Plot for the selected model with line break
+ unconstrained_empirical_quantiles = self.compute_empirical_quantiles(self.unconstrained_estimator)
+ # ax.plot(unconstrained_empirical_quantiles, sorted_maxima,
+ # label='Selected model, which is ${}$'.format(self.label))
+ # Plot tor the selected model for different year
+
+ end_real_proba = 1 - (0.02 / psnow)
+ self.plot_model(ax, 1958, end_proba=end_real_proba, label='Selected model, which is ${}$'.format(self.label))
+ self.plot_model(ax, 2017, end_proba=end_real_proba, label='Selected model, which is ${}$'.format(self.label))
+
# Plot for the discarded model
# if 'Verdon' in massif_name and altitude == 300:
# q = [-1.4541688117485054, -1.2811308174310914, -1.216589300814509, -0.7635793791201918, -0.6298883422064275,
@@ -252,12 +257,7 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
# + '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
@@ -273,14 +273,19 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
ax.set_ylabel("Non-zero annual maxima of GSL ({})".format(AbstractSnowLoadVariable.UNIT), fontsize=size)
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)
+ def plot_model(self, ax, year, start_proba=0.02, end_proba=0.98, label=''):
+ standard_gumbel = GevParams(0, 1, 0)
+ start_quantile = standard_gumbel.quantile(start_proba)
+ end_quantile = standard_gumbel.quantile(end_proba)
+ extended_quantiles = np.linspace(start_quantile, end_quantile, 500)
+ if year is None:
+ year = 2017
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))
+ label = 'Y({})'.format(year) if year is not None else label
+ ax.plot(extended_quantiles, extended_maxima, linestyle='-', label=label)
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
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 f8efbe7a..4f89912a 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
@@ -99,7 +99,7 @@ for the worst example for -shape
"""
if __name__ == '__main__':
- fast = True
+ 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)
--
GitLab