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 c6040667a23f529945ab815d145485a0b950f4df..2aed224f01156565b6cf1ebdda7e2415bb12b302 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
@@ -1,7 +1,7 @@
-import numpy as np
from math import ceil, floor
+
import matplotlib.pyplot as plt
-import pandas as pd
+import numpy as np
from cached_property import cached_property
from scipy.stats import chi2
@@ -9,13 +9,12 @@ from experiment.eurocode_data.utils import EUROCODE_QUANTILE
from experiment.trend_analysis.univariate_test.abstract_univariate_test import AbstractUnivariateTest
from experiment.trend_analysis.univariate_test.utils import load_temporal_coordinates_and_dataset, \
fitted_linear_margin_estimator
-from extreme_fit.estimator.margin_estimator.abstract_margin_estimator import LinearMarginEstimator
+from extreme_fit.distribution.gev.gev_params import GevParams
from extreme_fit.model.margin_model.linear_margin_model.abstract_temporal_linear_margin_model import \
- AbstractTemporalLinearMarginModel, TemporalMarginFitMethod
+ TemporalMarginFitMethod
from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import \
StationaryTemporalModel
from extreme_fit.model.utils import SafeRunException
-from extreme_fit.distribution.gev.gev_params import GevParams
from root_utils import classproperty
from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
@@ -240,25 +239,58 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
plt.show()
+ def compute_empirical_quantiles(self, estimator):
+ empirical_quantiles = []
+ for year, maximum in sorted(zip(self.years, self.maxima), key=lambda t: t[1]):
+ gev_param = estimator.margin_function_from_fit.get_gev_params_with_big_shape_and_correct_shape(
+ coordinate=np.array([year]),
+ is_transformed=False)
+ maximum_standardized = gev_param.gumbel_standardization(maximum)
+ empirical_quantiles.append(maximum_standardized)
+ return empirical_quantiles
+
+ # For some visualizations
+
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, gev_params_with_corrected_shape = self.get_gev_params_with_big_shape_and_correct_shape()
+ 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 return_level_plot_difference(self, ax, ax2, label, color=None):
+ gev_params, gev_params_with_corrected_shape = self.get_gev_params_with_big_shape_and_correct_shape()
+ return_periods = list(range(2, 61))
+ quantile_with_big_shape = gev_params.get_return_level(return_periods)
+ quantile_with_small_shape = gev_params_with_corrected_shape.get_return_level(return_periods)
+ difference_quantile = quantile_with_big_shape - quantile_with_small_shape
+ relative_difference = 100 * difference_quantile / quantile_with_small_shape
+ # Plot the difference on the two axis
+
+ ax.vlines(50, 0, np.max(difference_quantile))
+ ax.plot(return_periods, difference_quantile, color=color, linestyle='-', label=label)
+ ax.legend(loc='upper left')
+ difference_ylabel = 'difference return level in 2017'
+ ax.set_ylabel(difference_ylabel + ' (kPa)')
+
+ ax2.vlines(50, 0, np.max(relative_difference))
+ ax2.plot(return_periods, relative_difference, color=color, linestyle='--', label=label)
+ ax2.legend(loc='upper right')
+ ax2.yaxis.grid()
+ ax2.set_ylabel('relative ' + difference_ylabel + ' (%)')
+
+ ax.set_xlabel('Return period')
+ ax.set_xticks([10 * i for i in range(1, 7)])
+ plt.gca().set_ylim(bottom=0)
+
+ def get_gev_params_with_big_shape_and_correct_shape(self):
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]):
- gev_param = estimator.margin_function_from_fit.get_gev_params(
- coordinate=np.array([year]),
- is_transformed=False)
- maximum_standardized = gev_param.gumbel_standardization(maximum)
- empirical_quantiles.append(maximum_standardized)
- return empirical_quantiles
+ return gev_params, gev_params_with_corrected_shape
diff --git a/extreme_fit/distribution/gev/gev_params.py b/extreme_fit/distribution/gev/gev_params.py
index f5bab17e1a9d91e8d47e38a3b4b223e8d8bf71eb..df722015df3ee9c0f8206ad35e795b833b1a45a8 100644
--- a/extreme_fit/distribution/gev/gev_params.py
+++ b/extreme_fit/distribution/gev/gev_params.py
@@ -157,7 +157,7 @@ class GevParams(AbstractParams):
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]
+ quantiles = self.get_return_level(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)
@@ -169,3 +169,6 @@ class GevParams(AbstractParams):
plt.gca().set_ylim(bottom=0)
if show:
plt.show()
+
+ def get_return_level(self, return_periods):
+ return np.array([self.quantile(1 - 1 / return_period) for return_period in return_periods])
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 291851f2d56f1a29a7ddcb8bf8aaf20c9c698796..bccf350acf08a53ebf1f346b71aea6d1e56c13f7 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
@@ -45,10 +45,11 @@ def plot_return_level_for_time_series_with_big_shape_parameters(tuple_ordered_by
l = tuple_ordered_by_shape
print('Highest examples:')
ax = plt.gca()
+ ax2 = ax.twinx()
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)
+ v.return_level_plot(ax, ax2, m, color)
plt.show()
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 34be06091da90a40bd9ae8408ea2e103da08cf6e..9410781997878dc8d5da72f07825eece5bd68102 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
@@ -397,11 +397,11 @@ 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):
+ def return_level_plot(self, ax, ax2, 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)
+ trend_test.return_level_plot_difference(ax, ax2, label, color)
# Part 4 - Trend plot