Skip to content
GitLab
Commit 8be09426 authored by Le Roux Erwan's avatar Le Roux Erwan
Browse files

[paper 1] add file to see the impact of mixed distribution

parent f869c2a4
master
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 115 additions and 6 deletions
+115 -6
papers/exceeding_snow_loads/check_mle_convergence_for_trends/mixed_distribution/mixed_distribution_impact.py 0 → 100644
+ 105
− 0
View file @ 8be09426
import pandas as pd
from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal
from experiment.meteo_france_data.scm_models_data.visualization.study_visualization.main_study_visualizer import \
ALL_ALTITUDES_WITHOUT_NAN
from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import \
ConfidenceIntervalMethodFromExtremes
from papers.exceeding_snow_loads.study_visualizer_for_non_stationary_trends import StudyVisualizerForNonStationaryTrends
def mix_dsitrbution_impact():
altitudes = [300, 600, 900]
altitude_to_couple = {altitude: [None, None] for altitude in altitudes}
for j, fit_gev_only_on_non_null_maxima in enumerate([True, False]):
for altitude in altitudes:
altitude_to_couple[altitude][j] = StudyVisualizerForNonStationaryTrends(
CrocusSnowLoadTotal(altitude=altitude),
select_only_acceptable_shape_parameter=True,
multiprocessing=True,
save_to_file=True,
fit_gev_only_on_non_null_maxima=fit_gev_only_on_non_null_maxima,
fit_only_time_series_with_ninety_percent_of_non_null_values=True,
show=False)
pd.options.display.width = 0
# Compare return level
for altitude, (viz1, viz2) in altitude_to_couple.items():
# for m in viz1.all_massif_name_to_eurocode_uncertainty_for_minimized_aic_model_class:
m_to_eurocode_return_level_1 = viz1.all_massif_name_to_eurocode_uncertainty_for_minimized_aic_model_class()
m_to_eurocode_return_level_2 = viz2.all_massif_name_to_eurocode_uncertainty_for_minimized_aic_model_class()
s_diff = []
massifs_names = set(viz1.massifs_names_with_year_without_snow).intersection(m_to_eurocode_return_level_1.keys())
for m in massifs_names:
return_level1 = m_to_eurocode_return_level_1[m].mean_estimate
return_level2 = m_to_eurocode_return_level_2[m].mean_estimate
difference = return_level1 - return_level2
abs_difference = abs(difference)
s_diff.append({'massif': m, 'abs(r1 - r2)': abs_difference, 'abs(r1 - r2)/r1': 100 * abs_difference/return_level1,
'abs(r1 - r2)/r2': 100 * abs_difference / return_level2,
'r1 (with mixed)': return_level1, 'r2 (without)': return_level2
})
df = pd.DataFrame(s_diff)
print(df)
print(altitude, '\n', df.describe())
"""
abs(r1 - r2) abs(r1 - r2)/r1 abs(r1 - r2)/r2 massif r1 (with mixed) r2 (without)
0 0.000000 0.000000 0.000000 Belledonne 0.370008 0.370008
1 0.000564 0.186687 0.186340 Bauges 0.302025 0.302589
2 0.000000 0.000000 0.000000 Mont-Blanc 0.407016 0.407016
3 0.017822 2.427636 2.370098 Aravis 0.734127 0.751949
4 0.025847 7.197779 6.714485 Chablais 0.359097 0.384944
5 0.060221 9.472780 8.653092 Vanoise 0.635731 0.695952
6 0.099971 21.315869 17.570553 Oisans 0.468997 0.568967
7 0.004269 1.617252 1.643837 Beaufortain 0.263942 0.259673
8 0.000731 0.206005 0.205581 Maurienne 0.354998 0.355729
300
abs(r1 - r2) abs(r1 - r2)/r1 abs(r1 - r2)/r2 r1 (with mixed) r2 (without)
count 9.000000 9.000000 9.000000 9.000000 9.000000
mean 0.023269 4.713779 4.149332 0.432882 0.455203
std 0.034915 7.110878 5.938525 0.156124 0.174963
min 0.000000 0.000000 0.000000 0.263942 0.259673
25% 0.000564 0.186687 0.186340 0.354998 0.355729
50% 0.004269 1.617252 1.643837 0.370008 0.384944
75% 0.025847 7.197779 6.714485 0.468997 0.568967
max 0.099971 21.315869 17.570553 0.734127 0.751949
abs(r1 - r2) abs(r1 - r2)/r1 abs(r1 - r2)/r2 massif r1 (with mixed) r2 (without)
0 0.0 0.0 0.0 Haut_Var-Haut_Verdon 1.977354 1.977354
1 0.0 0.0 0.0 Ubaye 1.607606 1.607606
2 0.0 0.0 0.0 Parpaillon 0.435603 0.435603
600
abs(r1 - r2) abs(r1 - r2)/r1 abs(r1 - r2)/r2 r1 (with mixed) r2 (without)
count 3.0 3.0 3.0 3.000000 3.000000
mean 0.0 0.0 0.0 1.340187 1.340187
std 0.0 0.0 0.0 0.804912 0.804912
min 0.0 0.0 0.0 0.435603 0.435603
25% 0.0 0.0 0.0 1.021604 1.021604
50% 0.0 0.0 0.0 1.607606 1.607606
75% 0.0 0.0 0.0 1.792480 1.792480
max 0.0 0.0 0.0 1.977354 1.977354
abs(r1 - r2) abs(r1 - r2)/r1 abs(r1 - r2)/r2 massif r1 (with mixed) r2 (without)
0 0.0 0.0 0.0 Mercantour 2.818572 2.818572
900
abs(r1 - r2) abs(r1 - r2)/r1 abs(r1 - r2)/r2 r1 (with mixed) r2 (without)
count 1.0 1.0 1.0 1.000000 1.000000
mean 0.0 0.0 0.0 2.818572 2.818572
std NaN NaN NaN NaN NaN
min 0.0 0.0 0.0 2.818572 2.818572
25% 0.0 0.0 0.0 2.818572 2.818572
50% 0.0 0.0 0.0 2.818572 2.818572
75% 0.0 0.0 0.0 2.818572 2.818572
max 0.0 0.0 0.0 2.818572 2.818572
Process finished with exit code 0
"""
"""
Au grand maximum: l ecart est inférieur à 0.1 kN.... entre si on utilise la mixed distribution ou non
"""
if __name__ == '__main__':
mix_dsitrbution_impact()
papers/exceeding_snow_loads/check_mle_convergence_for_trends/qqplot/plot_qqplot.py
+ 1
− 0
View file @ 8be09426
......@@ -193,6 +193,7 @@ if __name__ == '__main__':
fit_method=TemporalMarginFitMethod.extremes_fevd_mle,
multiprocessing=True,
fit_gev_only_on_non_null_maxima=True,
fit_only_time_series_with_ninety_percent_of_non_null_values=True,
save_to_file=True,
show=False)
for altitude in altitudes}
......
papers/exceeding_snow_loads/study_visualizer_for_non_stationary_trends.py
+ 9
− 6
View file @ 8be09426
......@@ -53,7 +53,7 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
non_stationary_trend_test_to_marker=None,
fit_method=TemporalMarginFitMethod.extremes_fevd_mle,
select_only_acceptable_shape_parameter=True,
fit_gev_only_on_non_null_maxima=True,
fit_gev_only_on_non_null_maxima=False,
fit_only_time_series_with_ninety_percent_of_non_null_values=True,
):
super().__init__(study, show, save_to_file, only_one_graph, only_first_row, vertical_kde_plot,
......@@ -126,8 +126,9 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
d[m] = (years[mask], maxima[mask])
else:
d = self.massif_name_to_years_and_maxima
if self.fit_only_time_series_with_ninety_percent_of_non_null_values:
d = {m: v for m, v in d.items() if self.massif_name_to_psnow[m] >= 0.90}
# In both cases, we remove any massif with psnow < 0.9
if self.fit_only_time_series_with_ninety_percent_of_non_null_values:
d = {m: v for m, v in d.items() if self.massif_name_to_psnow[m] >= 0.9}
return d
@property
......@@ -324,16 +325,18 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
else:
raise ValueError(model_subset_for_uncertainty)
def all_massif_name_to_eurocode_uncertainty_for_minimized_aic_model_class(self, ci_method,
model_subset_for_uncertainty) \
def all_massif_name_to_eurocode_uncertainty_for_minimized_aic_model_class(self, ci_method=ConfidenceIntervalMethodFromExtremes.ci_mle,
model_subset_for_uncertainty=ModelSubsetForUncertainty.non_stationary_gumbel_and_gev) \
-> Dict[str, EurocodeConfidenceIntervalFromExtremes]:
# Compute for the uncertainty massif names
massifs_names = set(self.massif_name_to_years_and_maxima_for_model_fitting.keys()).\
intersection(self.uncertainty_massif_names)
arguments = [
[self.massif_name_to_years_and_maxima_for_model_fitting[m],
self.massif_name_and_model_subset_to_model_class(m, model_subset_for_uncertainty),
ci_method, self.effective_temporal_covariate,
self.massif_name_to_eurocode_quantile_level_in_practice[m]
] for m in self.uncertainty_massif_names]
] for m in massifs_names]
if self.multiprocessing:
with Pool(NB_CORES) as p:
res = p.starmap(compute_eurocode_confidence_interval, arguments)
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment