From b081c3ec7e7f89f27899659b12335f094f001542 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Mon, 13 Jan 2020 17:56:33 +0100
Subject: [PATCH] [paper 1] add upper and lower bound into the gev param
object. Add script to run fit wihout maximum, and then compare the upper
bound with the maximum values that was removed. The test is ok for all
altitudes, the upper bound found was always higher than the maximum values
---
.../without_maximum/__init__.py | 0
.../main_fit_without_maximum.py | 20 ++++++
...study_visualizer_for_fit_witout_maximum.py | 61 +++++++++++++++++++
extreme_fit/distribution/gev/gev_params.py | 20 +++++-
4 files changed, 100 insertions(+), 1 deletion(-)
create mode 100644 experiment/paper_past_snow_loads/check_mle_convergence_for_trends/without_maximum/__init__.py
create mode 100644 experiment/paper_past_snow_loads/check_mle_convergence_for_trends/without_maximum/main_fit_without_maximum.py
create mode 100644 experiment/paper_past_snow_loads/check_mle_convergence_for_trends/without_maximum/study_visualizer_for_fit_witout_maximum.py
diff --git a/experiment/paper_past_snow_loads/check_mle_convergence_for_trends/without_maximum/__init__.py b/experiment/paper_past_snow_loads/check_mle_convergence_for_trends/without_maximum/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/experiment/paper_past_snow_loads/check_mle_convergence_for_trends/without_maximum/main_fit_without_maximum.py b/experiment/paper_past_snow_loads/check_mle_convergence_for_trends/without_maximum/main_fit_without_maximum.py
new file mode 100644
index 00000000..7184ee57
--- /dev/null
+++ b/experiment/paper_past_snow_loads/check_mle_convergence_for_trends/without_maximum/main_fit_without_maximum.py
@@ -0,0 +1,20 @@
+from typing import Dict
+
+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 experiment.paper_past_snow_loads.check_mle_convergence_for_trends.without_maximum.study_visualizer_for_fit_witout_maximum import \
+ StudyVisualizerForFitWithoutMaximum
+
+
+def fit_without_maximum_value(altitude_to_visualizer: Dict[int, StudyVisualizerForFitWithoutMaximum]):
+ for v in altitude_to_visualizer.values():
+ v.maximum_value_test()
+
+
+if __name__ == '__main__':
+ altitudes = ALL_ALTITUDES_WITHOUT_NAN[:]
+ altitude_to_visualizer = {altitude: StudyVisualizerForFitWithoutMaximum(CrocusSnowLoadTotal(altitude=altitude),
+ multiprocessing=True)
+ for altitude in altitudes}
+ fit_without_maximum_value(altitude_to_visualizer)
diff --git a/experiment/paper_past_snow_loads/check_mle_convergence_for_trends/without_maximum/study_visualizer_for_fit_witout_maximum.py b/experiment/paper_past_snow_loads/check_mle_convergence_for_trends/without_maximum/study_visualizer_for_fit_witout_maximum.py
new file mode 100644
index 00000000..898d5840
--- /dev/null
+++ b/experiment/paper_past_snow_loads/check_mle_convergence_for_trends/without_maximum/study_visualizer_for_fit_witout_maximum.py
@@ -0,0 +1,61 @@
+from typing import Dict, Tuple
+
+import matplotlib
+import numpy as np
+from cached_property import cached_property
+
+from experiment.meteo_france_data.plot.create_shifted_cmap import get_shifted_map
+from experiment.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
+from experiment.paper_past_snow_loads.study_visualizer_for_non_stationary_trends import \
+ StudyVisualizerForNonStationaryTrends
+from experiment.trend_analysis.abstract_score import MeanScore
+
+
+class StudyVisualizerForFitWithoutMaximum(StudyVisualizerForNonStationaryTrends):
+
+ def __init__(self, study: AbstractStudy, **kwargs):
+ super().__init__(study, **kwargs)
+
+ @cached_property
+ def massif_name_to_maximum_index_for_non_null_values(self) -> Tuple[Dict, Dict]:
+ d = super().massif_name_to_non_null_years_and_maxima
+ d_without_maximum = {}
+ d_maximum = {}
+ for m, (years, maxima) in d.items():
+ idx = np.argmax(maxima)
+ maximum = maxima[idx]
+ maxima = np.delete(maxima, idx)
+ years = np.delete(years, idx)
+ d_without_maximum[m] = (years, maxima)
+ d_maximum[m] = maximum
+ return d_without_maximum, d_maximum
+
+ @property
+ def massif_name_to_maximum(self) -> Dict:
+ return self.massif_name_to_maximum_index_for_non_null_values[1]
+
+ @cached_property
+ def massif_name_to_non_null_years_and_maxima(self):
+ return self.massif_name_to_maximum_index_for_non_null_values[0]
+
+ def maximum_value_test(self):
+ diff = []
+ for massif_name, maximum in self.massif_name_to_maximum.items():
+ t = self.massif_name_to_trend_test_that_minimized_aic[massif_name]
+ msg = '{} {}m'.format(massif_name, self.study.altitude)
+ upper_bound = t.unconstrained_estimator_gev_params.density_upper_bound
+ if not np.isinf(upper_bound):
+ diff.append(upper_bound - maximum)
+ assert maximum <= upper_bound, msg
+ if len(diff) > 1:
+ print('{} mean difference={}'.format(self.study.altitude, min(diff)))
+
+
+
+
+
+
+
+
+
+
diff --git a/extreme_fit/distribution/gev/gev_params.py b/extreme_fit/distribution/gev/gev_params.py
index f19b54d2..eb957c35 100644
--- a/extreme_fit/distribution/gev/gev_params.py
+++ b/extreme_fit/distribution/gev/gev_params.py
@@ -130,4 +130,22 @@ class GevParams(AbstractParams):
if self.shape == 0:
return x
else:
- return np.log(1 + self.shape * x) / self.shape
\ No newline at end of file
+ return np.log(1 + self.shape * x) / self.shape
+
+ @property
+ def bound(self):
+ return self.location - (self.scale / self.shape)
+
+ @property
+ def density_upper_bound(self):
+ if self.shape >= 0:
+ return np.inf
+ else:
+ return self.bound
+
+ @property
+ def density_lower_bound(self):
+ if self.shape <= 0:
+ return np.inf
+ else:
+ return self.bound
\ No newline at end of file
--
GitLab