From 756977ccec241e381570af9ea362c0ea2085d9ca Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Thu, 27 Feb 2020 13:33:26 +0100
Subject: [PATCH] [paper 1] intensity plots. fit ylabel pour selection curves.
---
.../abstract_gev_trend_test.py | 106 +++++++++++++++++-
.../qqplot/plot_qqplot.py | 22 +++-
...study_visualizer_for_fit_witout_maximum.py | 2 +-
.../main_result_trends_and_return_levels.py | 16 +--
.../plot_selection_curves.py | 4 +-
...dy_visualizer_for_non_stationary_trends.py | 4 +
6 files changed, 134 insertions(+), 20 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 8d4fde60..d6a67436 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
@@ -6,6 +6,7 @@ from cached_property import cached_property
from scipy.stats import chi2
from experiment.eurocode_data.utils import EUROCODE_QUANTILE
+from experiment.meteo_france_data.scm_models_data.crocus.crocus_variables import AbstractSnowLoadVariable
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
@@ -209,13 +210,87 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
# Some display function
+ def intensity_plot_wrt_standard_gumbel(self, massif_name, altitude, psnow):
+ 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')
+
+ # 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, 100]
+ 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 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)
+
+ all_quantiles = standard_gumbel_quantiles + unconstrained_empirical_quantiles + quantiles
+ epsilon = 0.5
+ ax_lim = [min(all_quantiles) - epsilon, max(all_quantiles) + epsilon]
+ ax.set_xlim(ax_lim)
+ ax_twiny.set_xlim(ax_lim)
+
+ ax.legend()
+
+ 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})
+ plt.show()
+
def qqplot_wrt_standard_gumbel(self, massif_name, altitude):
ax = plt.gca()
size = 15
- # Standard Gumbel quantiles
- standard_gumbel_distribution = GevParams(loc=0, scale=1, shape=0)
- n = len(self.years)
- standard_gumbel_quantiles = [standard_gumbel_distribution.quantile(i / (n + 1)) for i in range(1, n + 1)]
+ standard_gumbel_quantiles = self.get_standard_gumbel_quantiles()
unconstrained_empirical_quantiles = self.compute_empirical_quantiles(self.unconstrained_estimator)
constrained_empirical_quantiles = self.compute_empirical_quantiles(self.constrained_estimator)
all_quantiles = standard_gumbel_quantiles + unconstrained_empirical_quantiles + constrained_empirical_quantiles
@@ -230,7 +305,14 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
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]
+ 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
@@ -251,6 +333,20 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
plt.show()
+ def get_standard_gumbel_quantiles(self):
+ # Standard Gumbel quantiles
+ standard_gumbel_distribution = GevParams(loc=0, scale=1, shape=0)
+ n = len(self.years)
+ standard_gumbel_quantiles = [standard_gumbel_distribution.quantile(i / (n + 1)) for i in range(1, n + 1)]
+ 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]
+ return corresponding_quantiles
+
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/papers/exceeding_snow_loads/check_mle_convergence_for_trends/qqplot/plot_qqplot.py b/papers/exceeding_snow_loads/check_mle_convergence_for_trends/qqplot/plot_qqplot.py
index 48e1b663..4fc33b67 100644
--- a/papers/exceeding_snow_loads/check_mle_convergence_for_trends/qqplot/plot_qqplot.py
+++ b/papers/exceeding_snow_loads/check_mle_convergence_for_trends/qqplot/plot_qqplot.py
@@ -16,9 +16,7 @@ from papers.exceeding_snow_loads.study_visualizer_for_non_stationary_trends impo
from extreme_fit.distribution.gev.gev_params import GevParams
-def plot_qqplot_for_time_series_with_missing_zeros(
- altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends],
- nb_worst_examples=3):
+def extract_time_serimes_with_worst_number_of_zeros(altitude_to_visualizer, nb_worst_examples):
# Extract all the values
l = []
for a, v in altitude_to_visualizer.items():
@@ -29,9 +27,25 @@ def plot_qqplot_for_time_series_with_missing_zeros(
for a, v, m, p in l:
print(a, m, p)
print('Last standard quantile (depends on the number of data):', last_quantile(p))
+ return l
+
+
+def plot_qqplot_for_time_series_with_missing_zeros(
+ altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends],
+ nb_worst_examples=3):
+ l = extract_time_serimes_with_worst_number_of_zeros(altitude_to_visualizer, nb_worst_examples)
+ for a, v, m, p in l:
v.qqplot(m)
+def plot_intensity_against_gumbel_quantile_for_time_series_with_missing_zeros(
+ altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends],
+ nb_worst_examples=3):
+ l = extract_time_serimes_with_worst_number_of_zeros(altitude_to_visualizer, nb_worst_examples)
+ for a, v, m, p in l:
+ v.intensity_plot(m, p)
+
+
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:
@@ -109,4 +123,4 @@ if __name__ == '__main__':
# plot_qqplot_wrt_standard_gumbel(altitude_to_visualizer)
# plot_hist_psnow(altitude_to_visualizer)
# plot_qqplot_for_time_series_examples(altitude_to_visualizer)
- plot_qqplot_for_time_series_with_missing_zeros(altitude_to_visualizer, nb_worst_examples=3)
+ plot_intensity_against_gumbel_quantile_for_time_series_with_missing_zeros(altitude_to_visualizer, nb_worst_examples=3)
diff --git a/papers/exceeding_snow_loads/check_mle_convergence_for_trends/without_maximum/study_visualizer_for_fit_witout_maximum.py b/papers/exceeding_snow_loads/check_mle_convergence_for_trends/without_maximum/study_visualizer_for_fit_witout_maximum.py
index 6a929d1f..ea904aba 100644
--- a/papers/exceeding_snow_loads/check_mle_convergence_for_trends/without_maximum/study_visualizer_for_fit_witout_maximum.py
+++ b/papers/exceeding_snow_loads/check_mle_convergence_for_trends/without_maximum/study_visualizer_for_fit_witout_maximum.py
@@ -6,7 +6,7 @@ from cached_property import cached_property
from experiment.meteo_france_data.plot.create_shifted_cmap import get_shifted_map
from experiment.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
-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
from experiment.trend_analysis.abstract_score import MeanScore
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 3160bed2..4ee5af9c 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
@@ -68,10 +68,10 @@ def intermediate_result(altitudes, massif_names=None,
# Plots
# plot_trend_map(altitude_to_visualizer)
# plot_diagnosis_risk(altitude_to_visualizer)
- # 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_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)
def major_result():
@@ -91,11 +91,11 @@ def major_result():
if __name__ == '__main__':
- # major_result()
- intermediate_result(altitudes=paper_altitudes, massif_names=['Beaufortain', 'Vercors'],
- uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
- ConfidenceIntervalMethodFromExtremes.ci_mle][1:],
- multiprocessing=True)
+ major_result()
+ # intermediate_result(altitudes=[900], massif_names=['Beaufortain', 'Vercors'],
+ # uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
+ # ConfidenceIntervalMethodFromExtremes.ci_mle][1:],
+ # multiprocessing=True)
# intermediate_result(altitudes=[900, 1200], massif_names=['Maurienne'],
# uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
# ConfidenceIntervalMethodFromExtremes.ci_mle][1:],
diff --git a/papers/exceeding_snow_loads/result_trends_and_return_levels/plot_selection_curves.py b/papers/exceeding_snow_loads/result_trends_and_return_levels/plot_selection_curves.py
index 910cdb99..cb17d737 100644
--- a/papers/exceeding_snow_loads/result_trends_and_return_levels/plot_selection_curves.py
+++ b/papers/exceeding_snow_loads/result_trends_and_return_levels/plot_selection_curves.py
@@ -51,8 +51,8 @@ def plot_selection_curves(altitude_to_visualizer: Dict[int, StudyVisualizerForNo
label='Significative model',
linewidth=linewidth)
ax.legend(loc='upper right', prop={'size': size})
- ax.set_ylabel('Time series where the model is selected\n'
- 'i.e. where it minimizes the AIC (\%)', fontsize=legend_fontsize)
+ ax.set_ylabel(' Frequency of selected model w.r.t all time series \n '
+ 'i.e. for all massifs and altitudes (\%)', fontsize=legend_fontsize)
ax.set_xlabel('Models', fontsize=legend_fontsize)
ax.tick_params(axis='both', which='major', labelsize=labelsize)
ax.set_xticks(x)
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 52af41d3..7885689d 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
@@ -392,6 +392,10 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
# Part 3 - QQPLOT
+ 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)
+
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)
--
GitLab