diff --git a/extreme_trend/abstract_gev_trend_test.py b/extreme_trend/abstract_gev_trend_test.py
index 57f8921bc81ffd34f4c11a7696f83b36faa9e927..66b7b889626c49bcf21c2ecde4a3808898593c31 100644
--- a/extreme_trend/abstract_gev_trend_test.py
+++ b/extreme_trend/abstract_gev_trend_test.py
@@ -14,7 +14,7 @@ from extreme_fit.distribution.gev.gev_params import GevParams
from extreme_fit.model.margin_model.utils import \
MarginFitMethod
from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import \
- StationaryTemporalModel
+ StationaryTemporalModel, GumbelTemporalModel
from extreme_fit.model.utils import SafeRunException
from root_utils import classproperty
from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
@@ -50,6 +50,10 @@ class AbstractGevTrendTest(object):
return fitted_linear_margin_estimator(self.unconstrained_model_class, self.coordinates, self.dataset,
self.starting_year, self.fit_method)
+ @property
+ def unconstrained_model_is_stationary(self):
+ return self.unconstrained_model_class in [StationaryTemporalModel, GumbelTemporalModel]
+
# Likelihood ratio test
@property
@@ -168,6 +172,58 @@ class AbstractGevTrendTest(object):
# Some display function
+ def intensity_plot_wrt_standard_gumbel_simple(self, massif_name, altitude, color='grey', psnow=1):
+ ax = plt.gca()
+ sorted_maxima = sorted(self.maxima)
+ label_generic = '{} massif \nat {} m '.format(massif_name, altitude)
+ size = 15
+
+ # Plot for the empirical
+ standard_gumbel_quantiles = self.get_standard_gumbel_quantiles()
+ ax.plot(standard_gumbel_quantiles, sorted_maxima, linestyle='None',
+ label='Empirical model', marker='o', color='black')
+
+ ax_twiny = ax.twiny()
+
+ # 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))
+
+ end_real_proba = 1 - (0.02 / psnow)
+ if self.unconstrained_model_is_stationary:
+ years = [None]
+ else:
+ color = None
+ years = [1959, 1989, 2019]
+ for j, year in enumerate(years):
+ year_suffix = '' if year is None else 'for the year t={}'.format(year)
+ if j == 0:
+ label = 'Selected model i.e. ${}$'.format(self.label)
+ if year is not None:
+ 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]
+ ax.set_xlim(ax_lim)
+ ax_twiny.set_xlim(ax_lim)
+ ax.set_xticks([-1 + i for i in range(6)])
+ epsilon = 0.005
+ ax.set_ylim(bottom=-epsilon)
+ lalsize = 13
+ ax.tick_params(axis='both', which='major', labelsize=lalsize)
+ ax_twiny.tick_params(axis='both', which='major', labelsize=lalsize)
+
+ 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(prop={'size': 17})
+
+
def intensity_plot_wrt_standard_gumbel(self, massif_name, altitude, psnow):
ax = plt.gca()
sorted_maxima = sorted(self.maxima)
@@ -238,13 +294,14 @@ class AbstractGevTrendTest(object):
ax.legend(prop={'size': 17})
def plot_model(self, ax, year, start_proba=0.02, end_proba=0.98, label='', color=None):
+ if year is None:
+ # for the stationary model
+ year = 2019
+
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)
- label = 'Y({})'.format(year) if year is not None else label
- if year is None:
- year = 2019
gev_params_year = self.get_unconstrained_gev_params(year)
extended_maxima = [gev_params_year.gumbel_inverse_standardization(q) for q in extended_quantiles]
@@ -301,21 +358,6 @@ class AbstractGevTrendTest(object):
label_generic = '{} massif \nat {} m '.format(massif_name, altitude)
ax.plot(standard_gumbel_quantiles, unconstrained_empirical_quantiles, linestyle='None',
label=label_generic + '(selected model is ${}$)'.format(self.label), marker='o')
- 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]
- print(len(q), len(standard_gumbel_quantiles))
- ax.plot(standard_gumbel_quantiles, q, linestyle='None',
- label=label_generic
- + '(discarded model is ${}$\n'.format('\mathcal{M}_{\zeta_0, \sigma_1}')
- + 'with $\zeta_0=0.84$)',
- marker='o')
ax.set_xlabel("Standard Gumbel quantile", fontsize=size)
ax.set_ylabel("Standard Empirical quantile", fontsize=size)
@@ -328,8 +370,6 @@ class AbstractGevTrendTest(object):
ax.grid()
ax.tick_params(labelsize=size)
- plt.show()
-
def get_standard_gumbel_quantiles(self):
# Standard Gumbel quantiles
standard_gumbel_distribution = GevParams(loc=0, scale=1, shape=0)
@@ -338,7 +378,6 @@ class AbstractGevTrendTest(object):
return standard_gumbel_quantiles
def get_standard_quantiles_for_return_periods(self, return_periods, psnow):
- n = len(self.years)
p_list = [1 - ((1 / return_period) / psnow) for return_period in return_periods]
standard_gumbel_distribution = GevParams(loc=0, scale=1, shape=0)
corresponding_quantiles = [standard_gumbel_distribution.quantile(p) for p in p_list]
diff --git a/extreme_trend/visualizers/study_visualizer_for_non_stationary_trends.py b/extreme_trend/visualizers/study_visualizer_for_non_stationary_trends.py
index 419cca60583e5954ef8fc76357f8127dc736c4a5..b91631b3304b77e614d4733f6c9a84b7e553fb84 100644
--- a/extreme_trend/visualizers/study_visualizer_for_non_stationary_trends.py
+++ b/extreme_trend/visualizers/study_visualizer_for_non_stationary_trends.py
@@ -420,9 +420,10 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
# Part 3 - QQPLOT
- def intensity_plot(self, massif_name, psnow, color=None):
+ def intensity_plot(self, massif_name, color=None):
+ psnow = self.massif_name_to_psnow[massif_name]
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)
+ trend_test.intensity_plot_wrt_standard_gumbel_simple(massif_name, self.altitude, color, 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()
@@ -430,6 +431,9 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
def qqplot(self, massif_name, color=None):
trend_test = self.massif_name_to_trend_test_that_minimized_aic[massif_name]
trend_test.qqplot_wrt_standard_gumbel(massif_name, self.altitude)
+ self.plot_name = 'qpplot_plot_{}_{}_{}'.format(self.altitude, massif_name, trend_test.unconstrained_estimator_gev_params.shape)
+ self.show_or_save_to_file(add_classic_title=False, no_title=True)
+ plt.close()
def return_level_plot(self, ax, ax2, massif_name, color=None):
trend_test = self.massif_name_to_trend_test_that_minimized_aic[massif_name]
diff --git a/projects/exceeding_snow_loads/checks/qqplot/plot_qqplot.py b/projects/exceeding_snow_loads/checks/qqplot/plot_qqplot.py
index b0e9901fa30c9f6cac5b3d3ec17a5e3b94cb8fb8..529516725d70a0ae895202fcdc7e506ca13bc04c 100644
--- a/projects/exceeding_snow_loads/checks/qqplot/plot_qqplot.py
+++ b/projects/exceeding_snow_loads/checks/qqplot/plot_qqplot.py
@@ -11,10 +11,10 @@ from extreme_fit.model.margin_model.utils import \
MarginFitMethod
from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_extract_eurocode_return_level import \
AbstractExtractEurocodeReturnLevel
-from projects.exceeding_snow_loads.data.main_example_swe_total_plot import tuples_for_examples_paper1
from extreme_trend.visualizers.study_visualizer_for_non_stationary_trends import \
StudyVisualizerForNonStationaryTrends
from extreme_fit.distribution.gev.gev_params import GevParams
+from projects.exceeding_snow_loads.section_data.main_example_swe_total_plot import tuples_for_examples_paper1
def extract_time_serimes_with_worst_number_of_zeros(altitude_to_visualizer, nb_worst_examples):
@@ -49,6 +49,14 @@ def plot_intensity_against_gumbel_quantile_for_time_series_with_missing_zeros(
v.intensity_plot(m, p)
+def plot_intensity_against_gumbel_quantile_for_3_examples(
+ altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends]):
+ for color, altitude, m in tuples_for_examples_paper1():
+ v = altitude_to_visualizer[altitude]
+ v.intensity_plot(m, color)
+ v.qqplot(m, color)
+
+
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 d777a99b6eb2ad17d9f912f96a17b312123bce81..6e06ab37e33803a036c11e36d1fb9372275466ed 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
@@ -2,6 +2,8 @@ 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
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
@@ -69,8 +71,9 @@ 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)
def major_result():
@@ -86,7 +89,7 @@ def major_result():
model_subsets_for_uncertainty = None
# study_classes = [CrocusSnowLoad3Days, CrocusSnowLoad5Days, CrocusSnowLoad7Days][::-1]
altitudes = paper_altitudes
- # altitudes = [900]
+ altitudes = [900, 1800, 2700]
for study_class in study_classes:
intermediate_result(altitudes, massif_names, model_subsets_for_uncertainty,
uncertainty_methods, study_class,