diff --git a/papers/exceeding_snow_loads/check_mle_convergence_for_trends/uncertainty_interval_size.py b/papers/exceeding_snow_loads/check_mle_convergence_for_trends/uncertainty_interval_size.py
new file mode 100644
index 0000000000000000000000000000000000000000..13767aa047818072bf42da928c967c39895463c5
--- /dev/null
+++ b/papers/exceeding_snow_loads/check_mle_convergence_for_trends/uncertainty_interval_size.py
@@ -0,0 +1,29 @@
+from typing import Dict
+
+import pandas as pd
+
+from experiment.eurocode_data.utils import EUROCODE_ALTITUDES
+from papers.exceeding_snow_loads.paper_utils import ModelSubsetForUncertainty
+from papers.exceeding_snow_loads.study_visualizer_for_non_stationary_trends import StudyVisualizerForNonStationaryTrends
+
+
+def uncertainty_interval_size(altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends]):
+ """ Plot one graph for each non-stationary context
+ :return:
+ """
+ altitude_to_visualizer = {a: v for a, v in altitude_to_visualizer.items() if a in EUROCODE_ALTITUDES}
+ visualizer = list(altitude_to_visualizer.values())[0]
+ for a, v in altitude_to_visualizer.items():
+ print(a)
+ interval_size(v)
+
+
+def interval_size(v: StudyVisualizerForNonStationaryTrends):
+ d = v.all_massif_name_to_eurocode_uncertainty_for_minimized_aic_model_class(
+ model_subset_for_uncertainty=ModelSubsetForUncertainty.stationary_gev)
+ # what we want is the confidence interval for the shape parameter
+ d = {m: [e.confidence_interval[0], e.confidence_interval[1], e.confidence_interval[1] - e.confidence_interval[0]]
+ for m, e in d.items()}
+ df = pd.DataFrame(d).transpose()
+ print((df.head()))
+ print(df.describe())
diff --git a/papers/exceeding_snow_loads/paper_utils.py b/papers/exceeding_snow_loads/paper_utils.py
index 801febc09d1e9187330458882cd994c55d235ff8..35e2c51921deab72ebb18c4d2cc303e409994f2e 100644
--- a/papers/exceeding_snow_loads/paper_utils.py
+++ b/papers/exceeding_snow_loads/paper_utils.py
@@ -32,3 +32,4 @@ class ModelSubsetForUncertainty(Enum):
stationary_gumbel_and_gev = 1
non_stationary_gumbel = 2
non_stationary_gumbel_and_gev = 3
+ stationary_gev = 4
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 86ed3b600e33b926f211e16bf1d331e71a08cd8a..c656915be054ddd2fe4feed59b59db790cf7ddbd 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
@@ -2,7 +2,10 @@ from multiprocessing.pool import Pool
import matplotlib as mpl
-from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal, CrocusSnowLoad3Days
+from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal, CrocusSnowLoad3Days, \
+ CrocusSnowLoad5Days, CrocusSnowLoad7Days
+from papers.exceeding_snow_loads.check_mle_convergence_for_trends.uncertainty_interval_size import \
+ uncertainty_interval_size
from papers.exceeding_snow_loads.paper_main_utils import load_altitude_to_visualizer
from papers.exceeding_snow_loads.paper_utils import paper_study_classes, paper_altitudes
from papers.exceeding_snow_loads.result_trends_and_return_levels.plot_diagnosis_risk import plot_diagnosis_risk
@@ -71,7 +74,7 @@ def intermediate_result(altitudes, massif_names=None,
# plot_uncertainty_massifs(altitude_to_visualizer)
plot_uncertainty_histogram(altitude_to_visualizer)
# plot_selection_curves(altitude_to_visualizer)
-
+ # uncertainty_interval_size(altitude_to_visualizer)
def major_result():
@@ -92,7 +95,7 @@ def major_result():
if __name__ == '__main__':
major_result()
- # intermediate_result(altitudes=[1800], massif_names=None,
+ # intermediate_result(altitudes=[300], massif_names=None,
# uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
# ConfidenceIntervalMethodFromExtremes.ci_mle][1:],
# multiprocessing=True)
diff --git a/papers/exceeding_snow_loads/result_trends_and_return_levels/plot_uncertainty_histogram.py b/papers/exceeding_snow_loads/result_trends_and_return_levels/plot_uncertainty_histogram.py
index dea54bfdb06f96477d8b213f75cc66b029de0b59..f81a1c81628d9bdb5029c29c833e2715ccab882a 100644
--- a/papers/exceeding_snow_loads/result_trends_and_return_levels/plot_uncertainty_histogram.py
+++ b/papers/exceeding_snow_loads/result_trends_and_return_levels/plot_uncertainty_histogram.py
@@ -64,7 +64,6 @@ def plot_histogram(altitude_to_visualizer, model_subset_for_uncertainty):
ax_twiny.set_xlim(ax.get_xlim())
ax_twiny.set_xticks(altitudes)
nb_massif_names = [len(v.massif_names_fitted) for v in altitude_to_visualizer.values()]
- print(nb_massif_names)
ax_twiny.set_xticklabels(nb_massif_names)
ax_twiny.set_xlabel('Total number of massifs at each altitude (for the percentage)', fontsize=fontsize_label)
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 c2d1a3c2d6dd020a03cefb5e581c16adf94b32c4..68b8b75c67c6255d376c8f0b10d3b0262b11ed00 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
@@ -30,7 +30,8 @@ from experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_two
GumbelLocationAndScaleTrendTest
from extreme_fit.model.margin_model.linear_margin_model.abstract_temporal_linear_margin_model import \
TemporalMarginFitMethod
-from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import GumbelTemporalModel
+from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import GumbelTemporalModel, \
+ StationaryTemporalModel
from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import \
ConfidenceIntervalMethodFromExtremes
from extreme_fit.model.result_from_model_fit.result_from_extremes.eurocode_return_level_uncertainties import \
@@ -322,6 +323,8 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
def massif_name_and_model_subset_to_model_class(self, massif_name, model_subset_for_uncertainty):
if model_subset_for_uncertainty is ModelSubsetForUncertainty.stationary_gumbel:
return GumbelTemporalModel
+ if model_subset_for_uncertainty is ModelSubsetForUncertainty.stationary_gev:
+ return StationaryTemporalModel
elif model_subset_for_uncertainty is ModelSubsetForUncertainty.stationary_gumbel_and_gev:
return self.massif_name_to_stationary_trend_test_that_minimized_aic[massif_name].unconstrained_model_class
elif model_subset_for_uncertainty is ModelSubsetForUncertainty.non_stationary_gumbel:
@@ -395,7 +398,7 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
triplet = [(massif_name_to_eurocode_region[massif_name],
self.massif_name_to_eurocode_values[massif_name],
self.triplet_to_eurocode_uncertainty[(ci_method, model_subset_for_uncertainty, massif_name)])
- for massif_name in self.uncertainty_massif_names]
+ for massif_name in self.massif_names_fitted]
# First array for histogram
a = 100 * np.array([(uncertainty.confidence_interval[0] > eurocode,
uncertainty.mean_estimate > eurocode,