diff --git a/experiment/paper_past_snow_loads/check_mle_convergence_for_trends/qqplot/plot_qqplot.py b/experiment/paper_past_snow_loads/check_mle_convergence_for_trends/qqplot/plot_qqplot.py
index 1a743dcb6bab143ddd1f800c45d62206fac7ca05..1d2d859fc9953fd85cd11f8f3883985046031edd 100644
--- a/experiment/paper_past_snow_loads/check_mle_convergence_for_trends/qqplot/plot_qqplot.py
+++ b/experiment/paper_past_snow_loads/check_mle_convergence_for_trends/qqplot/plot_qqplot.py
@@ -11,6 +11,7 @@ from experiment.meteo_france_data.scm_models_data.visualization.study_visualizat
from experiment.paper_past_snow_loads.data.main_example_swe_total_plot import tuples_for_examples_paper1
from experiment.paper_past_snow_loads.study_visualizer_for_non_stationary_trends import \
StudyVisualizerForNonStationaryTrends
+from extreme_fit.distribution.gev.gev_params import GevParams
def plot_qqplot_for_time_series_with_missing_zeros(
@@ -25,6 +26,7 @@ def plot_qqplot_for_time_series_with_missing_zeros(
print('Worst examples:')
for a, v, m, p in l:
print(a, m, p)
+ print(last_quantile(p))
v.qqplot(m)
@@ -37,6 +39,7 @@ def plot_qqplot_for_time_series_examples(altitude_to_visualizer: Dict[int, Study
def plot_hist_psnow(altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends]):
"""Plot an histogram of psnow containing data from all the visualizers given as argument"""
+ ax = plt.gca()
# Gather the data
data = [list(v.massif_name_to_psnow.values()) for v in altitude_to_visualizer.values()]
data = list(chain.from_iterable(data))
@@ -46,27 +49,45 @@ def plot_hist_psnow(altitude_to_visualizer: Dict[int, StudyVisualizerForNonStati
print(percentage_of_one)
data = [d for d in data if d < 1]
# Plot histogram
- nb_bins = 13
+ nb_bins = 7
percentage = False
weights = [1 / len(data) for _ in data] if percentage else None
- plt.hist(data, bins=nb_bins, range=(0.35, 1), weights=weights)
- plt.xticks([0.05 * i + 0.35 for i in range(nb_bins + 1)])
+ count, *_ = ax.hist(data, bins=nb_bins, range=(0.3, 1), weights=weights, rwidth=0.8)
+ # Set x ticks for the histogram
+ ax.set_xticks([0.1 * i + 0.3 for i in range(nb_bins + 1)])
+ ax.set_xlim([0.3, 1])
+ size = 10
if weights:
plt.gca().yaxis.set_major_formatter(PercentFormatter(1))
- plt.xlabel('Distribution of P(Y > 0) when $\\neq 1$')
+ ax.set_xlabel('Probability to have an annual maxima equal to zero, i.e. P(Y > 0)', fontsize=size)
s = '%' if percentage else 'Number'
- plt.ylabel('{} of time series'.format(s))
+ ylabel = '{} of time series with some\nannual maxima equal to zero, i.e. with $P(Y > 0) < 1$'
+ ax.set_ylabel(ylabel.format(s), fontsize=size)
+ if not percentage:
+ # Set y ticks for the histogram
+ max_data = int(max(count))
+ ticks = [i for i in range(max_data + 2)]
+ ax.set_yticks(ticks)
+ ax.set_ylim([min(ticks), max(ticks)])
+ ax.tick_params(labelsize=size)
plt.show()
+def last_quantile(psnow):
+ n = int(60 * psnow)
+ last_proba = n / (n + 1)
+ standard_gumbel = GevParams(0.0, 1.0, 0.0)
+ return standard_gumbel.quantile(last_proba)
+
+
if __name__ == '__main__':
# altitudes = [300, 600, 900, 1200, 1500, 1800][:2]
- # altitudes = ALL_ALTITUDES_WITHOUT_NAN
- altitudes = [900, 1800, 2700]
+ altitudes = ALL_ALTITUDES_WITHOUT_NAN
+ # altitudes = [900, 1800, 2700]
altitude_to_visualizer = {altitude: StudyVisualizerForNonStationaryTrends(CrocusSnowLoadTotal(altitude=altitude),
multiprocessing=True)
for altitude in altitudes}
# plot_qqplot_wrt_standard_gumbel(altitude_to_visualizer)
- # plot_hist_psnow(altitude_to_visualizer)
- plot_qqplot_for_time_series_examples(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)
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 9a07a18e1c5fe518fe1cf04b80ac03ae18b447d3..ba647019d76d6e736f9980c6464e15c16141d4f5 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,4 +1,5 @@
import numpy as np
+from math import ceil, floor
import matplotlib.pyplot as plt
import pandas as pd
from cached_property import cached_property
@@ -210,20 +211,32 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
# Some display function
def qqplot_wrt_standard_gumbel(self, marker, color=None):
- size = 20
+ 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)]
unconstrained_empirical_quantiles = self.compute_empirical_quantiles(self.unconstrained_estimator)
constrained_empirical_quantiles = self.compute_empirical_quantiles(self.constrained_estimator)
- plt.plot(standard_gumbel_quantiles, standard_gumbel_quantiles, color=color)
- plt.plot(standard_gumbel_quantiles, constrained_empirical_quantiles, 'x', label='Gumbel model')
- plt.plot(standard_gumbel_quantiles, unconstrained_empirical_quantiles, linestyle='None',
- label='Selected model', **marker)
- plt.xlabel("Standard Gumbel quantiles", fontsize=15)
- plt.ylabel("Empirical quantiles", fontsize=15)
- plt.legend(loc='upper left', prop={'size': size})
+ all_quantiles = standard_gumbel_quantiles + unconstrained_empirical_quantiles + constrained_empirical_quantiles
+ epsilon = 0.5
+ ax_lim = [min(all_quantiles) - epsilon, max(all_quantiles) + epsilon]
+ ax.plot(standard_gumbel_quantiles, standard_gumbel_quantiles, color='k')
+ ax.plot(standard_gumbel_quantiles, constrained_empirical_quantiles, 'x', label='Stationary Gumbel model $\mathcal{M}_0$')
+ ax.plot(standard_gumbel_quantiles, unconstrained_empirical_quantiles, linestyle='None',
+ label='Selected model $\mathcal{M}_N$', **marker)
+ ax.set_xlabel("Standard Gumbel quantiles", fontsize=size)
+ ax.set_ylabel("Empirical quantiles", fontsize=size)
+ ax.legend(loc='upper left', prop={'size': 10})
+ ax.set_xlim(ax_lim)
+ ax.set_ylim(ax_lim)
+ ticks = [i for i in range(ceil(ax_lim[0]), floor(ax_lim[1]) + 1)]
+ ax.set_xticks(ticks)
+ ax.set_yticks(ticks)
+ ax.grid()
+ ax.tick_params(labelsize=size)
+
plt.show()
def compute_empirical_quantiles(self, estimator):