diff --git a/extreme_trend/abstract_gev_trend_test.py b/extreme_trend/abstract_gev_trend_test.py
index ce17e3c3a9f5e34931dde65bfc9af9e0c5d17f62..57f8921bc81ffd34f4c11a7696f83b36faa9e927 100644
--- a/extreme_trend/abstract_gev_trend_test.py
+++ b/extreme_trend/abstract_gev_trend_test.py
@@ -3,7 +3,9 @@ from math import ceil, floor
import matplotlib.pyplot as plt
import numpy as np
from cached_property import cached_property
-from scipy.stats import chi2
+from scipy.stats import chi2, kstest, anderson
+from scipy.stats.morestats import AndersonResult
+from scipy.stats.stats import KstestResult
from extreme_data.eurocode_data.utils import EUROCODE_QUANTILE, YEAR_OF_INTEREST_FOR_RETURN_LEVEL
from extreme_data.meteo_france_data.scm_models_data.crocus.crocus_variables import AbstractSnowLoadVariable
@@ -54,6 +56,22 @@ class AbstractGevTrendTest(object):
def is_significant(self) -> bool:
return self.likelihood_ratio > chi2.ppf(q=1 - self.SIGNIFICANCE_LEVEL, df=self.degree_freedom_chi2)
+ @property
+ def goodness_of_fit_ks_test(self):
+ quantiles = self.compute_empirical_quantiles(estimator=self.unconstrained_estimator)
+ test_res = kstest(rvs=quantiles, cdf="gumbel_l") # type: KstestResult
+ return test_res.pvalue < self.SIGNIFICANCE_LEVEL
+
+ @property
+ def goodness_of_fit_anderson_test(self):
+ assert self.SIGNIFICANCE_LEVEL == 0.05
+ # significance_level=array([25. , 10. , 5. , 2.5, 1. ]))
+ index_for_significance_level_5_percent = 4
+ quantiles = self.compute_empirical_quantiles(estimator=self.unconstrained_estimator)
+ test_res = anderson(quantiles, dist='gumbel') # type: AndersonResult
+ print(test_res)
+ return test_res.statistic > test_res.critical_values[index_for_significance_level_5_percent]
+
@property
def degree_freedom_chi2(self) -> int:
raise NotImplementedError
diff --git a/extreme_trend/visualizers/study_visualizer_for_non_stationary_trends.py b/extreme_trend/visualizers/study_visualizer_for_non_stationary_trends.py
index 75188fd51f9772547325653e3575bc20da90ac41..419cca60583e5954ef8fc76357f8127dc736c4a5 100644
--- a/extreme_trend/visualizers/study_visualizer_for_non_stationary_trends.py
+++ b/extreme_trend/visualizers/study_visualizer_for_non_stationary_trends.py
@@ -305,6 +305,18 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
return Counter(
[type(t) for t in self.massif_name_to_trend_test_that_minimized_aic.values() if t.is_significant])
+ @cached_property
+ def selected_and_anderson_goodness_of_fit_trend_test_class_counter(self):
+ return Counter(
+ [type(t) for t in self.massif_name_to_trend_test_that_minimized_aic.values()
+ if t.goodness_of_fit_anderson_test])
+
+ @cached_property
+ def selected_and_kstest_goodness_of_fit_trend_test_class_counter(self):
+ return Counter(
+ [type(t) for t in self.massif_name_to_trend_test_that_minimized_aic.values()
+ if t.goodness_of_fit_ks_test])
+
@cached_property
def massif_name_to_marker_style(self):
d = {}
diff --git a/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/main_snowfall_article.py b/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/main_snowfall_article.py
index f7137c76b5694c3065eec31699a256ea81d507cd..282300f3e25efa01346f78767c08e10318302401 100644
--- a/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/main_snowfall_article.py
+++ b/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/main_snowfall_article.py
@@ -74,8 +74,8 @@ def intermediate_result(altitudes, massif_names=None,
# Plots
validation_plot(altitude_to_visualizer, order_derivative=0)
validation_plot(altitude_to_visualizer, order_derivative=1)
- # plot_snowfall_mean(altitude_to_visualizer)
- # plot_snowfall_time_derivative_mean(altitude_to_visualizer)
+ plot_snowfall_mean(altitude_to_visualizer)
+ plot_snowfall_time_derivative_mean(altitude_to_visualizer)
def major_result():
diff --git a/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/validation_plot.py b/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/validation_plot.py
index 4847c939f6b725798ec440679f146cd6044c4763..32c30b15aa6aa229cb3d67cc69a916ff7b0736e1 100644
--- a/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/validation_plot.py
+++ b/projects/contrasting_trends_in_snow_loads/snowfall_versus_time_and_altitude/validation_plot.py
@@ -24,20 +24,23 @@ def validation_plot(altitude_to_visualizer: Dict[int, StudyVisualizerForMeanValu
altitude_to_relative_differences[altitude] = plot_function(visualizer)
study_visualizer.show_or_save_to_file(add_classic_title=False, dpi=500)
# Shoe plot with respect to the altitude.
- plot_shoe_relative_differences_distribution(altitude_to_relative_differences, altitudes, study_visualizer, order_derivative)
+ plot_shoe_relative_differences_distribution(altitude_to_relative_differences, altitudes, study_visualizer,
+ order_derivative)
study_visualizer.show_or_save_to_file(add_classic_title=False, dpi=500)
plt.close()
-def plot_shoe_relative_differences_distribution(altitude_to_relative_differences, altitudes, visualizer, order_derivative):
+def plot_shoe_relative_differences_distribution(altitude_to_relative_differences, altitudes, visualizer,
+ order_derivative):
study = visualizer.study
ax = plt.gca()
width = 150
ax.boxplot([altitude_to_relative_differences[a] for a in altitudes], positions=altitudes, widths=width)
ax.set_xlim([min(altitudes) - width, max(altitudes) + width])
moment = '' if order_derivative == 0 else 'time derivative of '
- ylabel = 'Relative difference of the {} model mean \n' \
- 'w.r.t. the {}empirical mean of {} (\%)'.format(moment, moment, SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)])
+ ylabel = 'Global relative difference of the {} model mean \n' \
+ 'w.r.t. the {}empirical mean of {} (\%)'.format(moment, moment,
+ SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)])
visualizer.plot_trends = ylabel
ax.set_ylabel(ylabel)
ax.set_xlabel('Altitude (m)')
@@ -57,14 +60,17 @@ def plot_relative_difference_map_order_zero(visualizer: StudyVisualizerForMeanVa
'for the ' + label, graduation=1)
return list(visualizer.massif_name_to_relative_difference_for_mean.values())
+
def plot_relative_difference_map_order_one(visualizer: StudyVisualizerForMeanValues):
study = visualizer.study
- label = ' time derivative of mean annual maxima of {} ({})'.format(SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)], study.variable_unit)
+ label = ' time derivative of mean annual maxima of {} ({})'.format(SCM_STUDY_CLASS_TO_ABBREVIATION[type(study)],
+ study.variable_unit)
visualizer.plot_abstract_fast(massif_name_to_value=visualizer.massif_name_to_change_ratio_in_empirical_mean,
label='Empirical' + label, negative_and_positive_values=False, graduation=0.5)
visualizer.plot_abstract_fast(massif_name_to_value=visualizer.massif_name_to_change_ratio_in_model_mean,
label='Model' + label, negative_and_positive_values=False, graduation=0.5)
- visualizer.plot_abstract_fast(massif_name_to_value=visualizer.massif_name_to_relative_difference_for_change_ratio_in_mean,
- label='Relative difference of the model mean w.r.t. the empirical mean \n'
- 'for the ' + label, graduation=1)
- return list(visualizer.massif_name_to_relative_difference_for_mean.values())
\ No newline at end of file
+ visualizer.plot_abstract_fast(
+ massif_name_to_value=visualizer.massif_name_to_relative_difference_for_change_ratio_in_mean,
+ label='Relative difference of the model mean w.r.t. the empirical mean \n'
+ 'for the ' + label, graduation=5)
+ return list(visualizer.massif_name_to_relative_difference_for_change_ratio_in_mean.values())
diff --git a/projects/exceeding_snow_loads/section_results/main_result_trends_and_return_levels.py b/projects/exceeding_snow_loads/section_results/main_result_trends_and_return_levels.py
index ab830a32368616df25ea7d967be1cb745e5637cf..d777a99b6eb2ad17d9f912f96a17b312123bce81 100644
--- a/projects/exceeding_snow_loads/section_results/main_result_trends_and_return_levels.py
+++ b/projects/exceeding_snow_loads/section_results/main_result_trends_and_return_levels.py
@@ -2,6 +2,7 @@ from multiprocessing.pool import Pool
import matplotlib as mpl
+from projects.exceeding_snow_loads.section_results.plot_selection_curves import plot_selection_curves
from projects.exceeding_snow_loads.section_results.plot_trend_curves import plot_trend_map
mpl.rcParams['text.usetex'] = True
@@ -64,11 +65,11 @@ def intermediate_result(altitudes, massif_names=None,
_ = compute_minimized_aic(visualizer)
# Plots
- plot_trend_map(altitude_to_visualizer)
+ # plot_trend_map(altitude_to_visualizer)
# plot_trend_curves(altitude_to_visualizer={a: v for a, v in altitude_to_visualizer.items() if a >= 900})
# plot_uncertainty_massifs(altitude_to_visualizer)
# plot_uncertainty_histogram(altitude_to_visualizer)
- # plot_selection_curves(altitude_to_visualizer)
+ plot_selection_curves(altitude_to_visualizer)
# uncertainty_interval_size(altitude_to_visualizer)
@@ -84,9 +85,12 @@ def major_result():
# ModelSubsetForUncertainty.non_stationary_gumbel_and_gev]
model_subsets_for_uncertainty = None
# study_classes = [CrocusSnowLoad3Days, CrocusSnowLoad5Days, CrocusSnowLoad7Days][::-1]
+ altitudes = paper_altitudes
+ # altitudes = [900]
for study_class in study_classes:
- intermediate_result(paper_altitudes, massif_names, model_subsets_for_uncertainty,
- uncertainty_methods, study_class)
+ intermediate_result(altitudes, massif_names, model_subsets_for_uncertainty,
+ uncertainty_methods, study_class,
+ multiprocessing=True)
if __name__ == '__main__':
diff --git a/projects/exceeding_snow_loads/section_results/plot_selection_curves.py b/projects/exceeding_snow_loads/section_results/plot_selection_curves.py
index 2b4ced66b3fe27be9069b3e2df222ac1abc17175..da02e8ecc517ba0d7993fcaa4daa3c2a6e63cdb3 100644
--- a/projects/exceeding_snow_loads/section_results/plot_selection_curves.py
+++ b/projects/exceeding_snow_loads/section_results/plot_selection_curves.py
@@ -21,6 +21,8 @@ def plot_selection_curves(altitude_to_visualizer: Dict[int, StudyVisualizerForNo
selected_counter = merge_counter([v.selected_trend_test_class_counter for v in altitude_to_visualizer.values()])
selected_and_significative_counter = merge_counter([v.selected_and_significative_trend_test_class_counter for v in altitude_to_visualizer.values()])
+ # selected_and_significative_counter = merge_counter([v.selected_and_anderson_goodness_of_fit_trend_test_class_counter for v in altitude_to_visualizer.values()])
+ # selected_and_significative_counter = merge_counter([v.selected_and_kstest_goodness_of_fit_trend_test_class_counter for v in altitude_to_visualizer.values()])
total_of_selected_models = sum(selected_counter.values())
l = sorted(enumerate(visualizer.non_stationary_trend_test), key=lambda e: selected_counter[e[1]])
permutation = [i for i, v in l][::-1]