diff --git a/extreme_data/meteo_france_data/scm_models_data/utils_function.py b/extreme_data/meteo_france_data/scm_models_data/utils_function.py
index 5ee9281a97f9720a13b3e5728773e2c58cf754fc..55440cde70d77ff77ca30c79005a939a1e39b8c2 100644
--- a/extreme_data/meteo_france_data/scm_models_data/utils_function.py
+++ b/extreme_data/meteo_france_data/scm_models_data/utils_function.py
@@ -42,12 +42,8 @@ def fitted_stationary_gev_with_uncertainty_interval(x_gev, fit_method=MarginFitM
return_level_list = ReturnLevelBootstrap(fit_method, model_class, starting_year, x_gev,
quantile_level).compute_all_return_level()
# Remove infinite return levels and return level
- len_before_remove = len(return_level_list)
- threshold = 2000
- return_level_list = [r for r in return_level_list if (not np.isinf(r)) and (r < threshold)]
- len_after_remove = len(return_level_list)
- if len_after_remove < len_before_remove:
- print('Nb of fit removed (inf or > {}:'.format(threshold), len_before_remove - len_after_remove)
+ # percentage_of_inf = 100 * sum([np.isinf(r) for r in return_level_list]) / len(return_level_list)
+ # print('Percentage of fit with inf = {} \%'.format(percentage_of_inf))
confidence_interval = tuple([np.quantile(return_level_list, q)
for q in AbstractExtractEurocodeReturnLevel.bottom_and_upper_quantile])
mean_estimate = gev_param.quantile(quantile_level)
diff --git a/extreme_fit/estimator/margin_estimator/abstract_margin_estimator.py b/extreme_fit/estimator/margin_estimator/abstract_margin_estimator.py
index f233fa0ae67efe0cd7d16f5984298ccdc6c8c93c..168d7a10e4fc0f2ca728cd7618e7ee4a5d626044 100644
--- a/extreme_fit/estimator/margin_estimator/abstract_margin_estimator.py
+++ b/extreme_fit/estimator/margin_estimator/abstract_margin_estimator.py
@@ -75,11 +75,16 @@ class LinearMarginEstimator(AbstractMarginEstimator):
assert not np.isinf(nllh)
return nllh
- def sorted_empirical_standard_gumbel_quantiles(self, split=Split.all):
+ def sorted_empirical_standard_gumbel_quantiles(self, split=Split.all, coordinate_for_filter=None):
sorted_empirical_quantiles = []
maxima_values = self.dataset.maxima_gev(split=split)
coordinate_values = self.dataset.df_coordinates(split=split).values
for maximum, coordinate in zip(maxima_values, coordinate_values):
+ if coordinate_for_filter is not None:
+ assert len(coordinate) == len(coordinate_for_filter)
+ keep = any([(f is not None) and (c == f) for c, f in zip(coordinate, coordinate_for_filter)])
+ if not keep:
+ continue
gev_param = self.function_from_fit.get_params(
coordinate=coordinate,
is_transformed=False)
diff --git a/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py b/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
index 0f660a59e366c246cc22693bceb44ea41df96864..c0199b7ce32d7afeba832305004eda96f096439f 100644
--- a/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
+++ b/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
@@ -29,27 +29,29 @@ from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitude_gr
from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.plot_total_aic import plot_individual_aic
from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day, SafranSnowfall3Days, \
- SafranSnowfall5Days, SafranSnowfall7Days, SafranDateFirstSnowfall, SafranPrecipitation1Day, SafranPrecipitation3Days
+ SafranSnowfall5Days, SafranSnowfall7Days, SafranDateFirstSnowfall, SafranPrecipitation1Day, \
+ SafranPrecipitation3Days, SafranSnowfallCenterOnDay1dayMeanRate, SafranSnowfallNotCenterOnDay1day
from extreme_data.meteo_france_data.scm_models_data.utils import Season
def main():
- study_classes = [SafranSnowfall2020, SafranSnowfall2019,
- SafranSnowfall1Day, SafranSnowfall3Days,
- SafranSnowfall5Days, SafranSnowfall7Days][:3]
+ study_classes = [SafranSnowfall1Day
+ , SafranSnowfall3Days,
+ SafranSnowfall5Days, SafranSnowfall7Days][:1]
seasons = [Season.annual, Season.winter, Season.spring, Season.automn][:1]
set_seed_for_test()
- model_must_pass_the_test = True
+ model_must_pass_the_test = False
+ AbstractExtractEurocodeReturnLevel.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY = 0.2
- fast = None
+ fast = True
if fast is None:
- massif_names = None
- altitudes_list = altitudes_for_groups[2:3]
+ massif_names = ['Vanoise']
+ altitudes_list = altitudes_for_groups[:]
elif fast:
AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP = 10
- massif_names = ['Vanoise', 'Haute-Maurienne', 'Vercors'][:2]
- altitudes_list = altitudes_for_groups[3:]
+ massif_names = ['Vanoise', 'Haute-Maurienne', 'Vercors'][:1]
+ altitudes_list = altitudes_for_groups[1:2]
else:
massif_names = None
altitudes_list = altitudes_for_groups[:]
diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
index 1e4ec4f0022f8f99f66da7405b3fcee796dfb2fa..c5cb88181887cb78171f752c78acbbd732913dda 100644
--- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
+++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
@@ -508,22 +508,22 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
def plot_qqplots(self):
for massif_name, one_fold_fit in self.massif_name_to_one_fold_fit.items():
ax = plt.gca()
- model_name = self.massif_name_to_best_name[massif_name]
altitudes = self.massif_name_to_massif_altitudes[massif_name]
massif_name_corrected = massif_name.replace('_', ' ')
- label = '{} for altitudes {}'.format(massif_name_corrected, ' & '.join([str(a) + 'm' for a in altitudes]))
-
-
all_quantiles = []
- standard_gumbel_quantiles = one_fold_fit.standard_gumbel_quantiles()
-
+ for altitude in self.studies.altitudes:
+ coordinate_for_filter = (altitude, None)
+ unconstrained_empirical_quantiles = one_fold_fit.best_estimator.sorted_empirical_standard_gumbel_quantiles(coordinate_for_filter=coordinate_for_filter)
+ n = len(unconstrained_empirical_quantiles)
+ assert n == 61
+ standard_gumbel_quantiles = one_fold_fit.standard_gumbel_quantiles(n=n)
+ ax.plot(standard_gumbel_quantiles, unconstrained_empirical_quantiles, linestyle='None',
+ label='{} m'.format(altitude), marker='o')
- unconstrained_empirical_quantiles = one_fold_fit.best_estimator.sorted_empirical_standard_gumbel_quantiles()
- all_quantiles = standard_gumbel_quantiles + unconstrained_empirical_quantiles
- ax.plot(standard_gumbel_quantiles, unconstrained_empirical_quantiles, linestyle='None',
- label=label + '\n(selected model is ${}$)'.format(model_name), marker='o')
+ all_quantiles.extend(standard_gumbel_quantiles)
+ all_quantiles.extend(unconstrained_empirical_quantiles)
size_label = 20
ax.set_xlabel("Theoretical quantile", fontsize=size_label)
@@ -538,7 +538,10 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
ticks = [i for i in range(ceil(ax_lim[0]), floor(ax_lim[1]) + 1)]
ax.set_xticks(ticks)
ax.set_yticks(ticks)
- ax.tick_params(labelsize=15)
+ labelsize = 15
+ ax.tick_params(labelsize=labelsize)
plot_name = 'qqplot/{}'.format(massif_name_corrected)
+ handles, labels = ax.get_legend_handles_labels()
+ ax.legend(handles[::-1], labels[::-1], prop={'size': labelsize})
self.studies.show_or_save_to_file(plot_name=plot_name, show=self.show, no_title=True)
plt.close()
diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py
index 75fcfd1b3fba77caba456cd07f7f7bc678973429..162e555928d27639c33d5d65d4f91fc91ad38284 100644
--- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py
+++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py
@@ -330,9 +330,10 @@ class OneFoldFit(object):
goodness_of_fit_anderson_test = goodness_of_fit_anderson(quantiles, self.SIGNIFICANCE_LEVEL)
return goodness_of_fit_anderson_test
- def standard_gumbel_quantiles(self):
+ def standard_gumbel_quantiles(self, n=None):
standard_gumbel_distribution = GevParams(loc=0, scale=1, shape=0)
- n = len(self.dataset.coordinates)
+ if n is None:
+ n = len(self.dataset.coordinates)
standard_gumbel_quantiles = [standard_gumbel_distribution.quantile(i / (n + 1)) for i in range(1, n + 1)]
return standard_gumbel_quantiles
diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/plot_total_aic.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/plot_total_aic.py
index 6ebcbd94da8c0a415c5261e94278a751abba583d..a8e3dffb24fb3e549525d4fbc4ad99ae92eaedeb 100644
--- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/plot_total_aic.py
+++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/plot_total_aic.py
@@ -15,8 +15,8 @@ from projects.exceeding_snow_loads.utils import dpi_paper1_figure
def plots(visualizer: AltitudesStudiesVisualizerForNonStationaryModels):
# visualizer.plot_shape_map()
- visualizer.plot_moments()
- # visualizer.plot_qqplots()
+ # visualizer.plot_moments()
+ visualizer.plot_qqplots()
# for std in [True, False]:
# visualizer.studies.plot_mean_maxima_against_altitude(std=std)
diff --git a/projects/altitude_spatial_model/altitudes_fit/plots/plot_coherence_curves.py b/projects/altitude_spatial_model/altitudes_fit/plots/plot_coherence_curves.py
index 9b7c74c7915bdacf1285bacb46df1024446e1148..d8b762225beb5cf457c55ede3f64b9d85e596f49 100644
--- a/projects/altitude_spatial_model/altitudes_fit/plots/plot_coherence_curves.py
+++ b/projects/altitude_spatial_model/altitudes_fit/plots/plot_coherence_curves.py
@@ -2,6 +2,8 @@ from typing import List
import matplotlib.pyplot as plt
import numpy as np
+from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_extract_eurocode_return_level import \
+ AbstractExtractEurocodeReturnLevel
from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitudes_studies_visualizer_for_non_stationary_models import \
AltitudesStudiesVisualizerForNonStationaryModels
from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.one_fold_fit import OneFoldFit
@@ -103,9 +105,9 @@ def plot_coherence_curve(ax, i, x_all_list, values_all_list, all_bound_list,
lower_bound, upper_bound = bounds
f = ax.fill_between if elevation_as_xaxis else ax.fill_betweenx
if legend and not legend_line:
- print('here2')
model_name = 'piecewise elevational-temporal models in 2019' if is_altitudinal else 'pointwise distributions'
- fill_label = "95\% confidence interval for the {}".format(model_name) if j == 0 else None
+ percentage = AbstractExtractEurocodeReturnLevel.percentage_confidence_interval
+ fill_label = "{}\% confidence interval for the {}".format(percentage, model_name) if j == 0 else None
f(x_list, lower_bound, upper_bound, color=color, alpha=0.2, label=fill_label)
else:
f(x_list, lower_bound, upper_bound, color=color, alpha=0.2)