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 5e862f5527a7315a8005810b12d49b2d67f435e4..e9ec74bd310c4bc72e189d475fa05dd45f156364 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 @@ -3,6 +3,7 @@ from typing import Dict import matplotlib.pyplot as plt import numpy as np +import pandas as pd from matplotlib.ticker import PercentFormatter from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal @@ -157,6 +158,14 @@ def last_quantile(psnow): return standard_gumbel.quantile(last_proba) +def non_stationarity_psnow(altitude_to_visualizer): + all_relative_change_in_psnow = [] + for a, v in altitude_to_visualizer.items(): + all_relative_change_in_psnow.extend(list(v.massif_name_to_relative_change_in_psnow.values())) + s = pd.Series(all_relative_change_in_psnow) + print(s.describe()) + + if __name__ == '__main__': """ Worst examples: @@ -172,7 +181,7 @@ if __name__ == '__main__': 300 Haut_Var-Haut_Verdon 0.8446498197950775 """ - altitudes = [300, 600, 900, 1200, 1500, 1800][:-2] + altitudes = [300, 600, 900, 1200, 1500, 1800][:] # altitudes = ALL_ALTITUDES_WITHOUT_NAN # altitudes = [900, 1800, 2700] altitude_to_visualizer = {altitude: StudyVisualizerForNonStationaryTrends(CrocusSnowLoadTotal(altitude=altitude), @@ -183,7 +192,8 @@ if __name__ == '__main__': # plot_qqplot_wrt_standard_gumbel(altitude_to_visualizer) # plot_hist_psnow(altitude_to_visualizer) - plot_exceedance_psnow(altitude_to_visualizer) + # plot_exceedance_psnow(altitude_to_visualizer) + non_stationarity_psnow(altitude_to_visualizer) # plot_qqplot_for_time_series_examples(altitude_to_visualizer) # 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/study_visualizer_for_non_stationary_trends.py b/papers/exceeding_snow_loads/study_visualizer_for_non_stationary_trends.py index 7d71e1b7a33c54be571b892697ef45878f6de42d..2b82833c377eccfdf26bacf7621983c6ed5d71d2 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 @@ -474,3 +474,11 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer): percentage = 100 * np.array(eurocode_uncertainty.triplet) / eurocode_value percentages.append(percentage) return np.round(np.mean(percentages, axis=0)) + + @property + def massif_name_to_relative_change_in_psnow(self): + def compute_relative_change_in_psnow(maxima): + maxima_before, maxima_after = maxima[:30], maxima[30:] + psnow_before, psnow_after = [np.count_nonzero(s) / len(s) for s in [maxima_before, maxima_after]] + return 100 * (psnow_after - psnow_before) / psnow_before + return {m: compute_relative_change_in_psnow(self.massif_name_to_years_and_maxima[m][1]) for m in self.massifs_names_with_year_without_snow}