From 80137755aba3d4db8691dba84d28006ba048a6d1 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Sat, 23 Nov 2019 14:48:07 +0100
Subject: [PATCH] [TREND ANALYSIS] put MLE extreme package method as default
mle method for trend analysis. add aic score in trend analysis. start
study_visualizer for non stationary trends.
---
.../scm_models_data/abstract_study.py | 2 +-
.../study_visualization/study_visualizer.py | 3 -
.../method/main_result.py | 15 ++++
...dy_visualizer_for_non_stationary_trends.py | 90 +++++++++++++++++++
.../abstract_gev_trend_test.py | 45 +++++++---
.../abstract_univariate_test.py | 2 +-
.../univariate_test_results.py | 2 +-
.../abstract_result_from_model_fit.py | 2 +-
.../abstract_result_from_extremes.py | 8 ++
.../result_from_ismev.py | 4 -
10 files changed, 152 insertions(+), 21 deletions(-)
create mode 100644 experiment/paper_past_snow_loads/method/main_result.py
diff --git a/experiment/meteo_france_data/scm_models_data/abstract_study.py b/experiment/meteo_france_data/scm_models_data/abstract_study.py
index 752d7093..ac9208b6 100644
--- a/experiment/meteo_france_data/scm_models_data/abstract_study.py
+++ b/experiment/meteo_france_data/scm_models_data/abstract_study.py
@@ -299,7 +299,7 @@ class AbstractStudy(object):
default_color_for_nan_values='w',
massif_name_to_color=None,
show_label=True,
- scaled=False,
+ scaled=True,
fontsize=7,
axis_off=False,
massif_name_to_hatch_boolean_list=None,
diff --git a/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/study_visualizer.py b/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/study_visualizer.py
index ba731820..2829bf6e 100644
--- a/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/study_visualizer.py
+++ b/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/study_visualizer.py
@@ -711,9 +711,6 @@ class StudyVisualizer(VisualizationParameters):
self.massif_id_to_smooth_maxima[massif_id] = (x, y)
return self.massif_id_to_smooth_maxima[massif_id]
- def visualize_brown_resnick_fit(self):
- pass
-
def visualize_linear_margin_fit(self, only_first_max_stable=False):
margin_class = LinearNonStationaryLocationMarginModel if self.temporal_non_stationarity else LinearStationaryMarginModel
plot_name = 'Full Likelihood with Linear marginals and max stable dependency structure'
diff --git a/experiment/paper_past_snow_loads/method/main_result.py b/experiment/paper_past_snow_loads/method/main_result.py
new file mode 100644
index 00000000..085564f9
--- /dev/null
+++ b/experiment/paper_past_snow_loads/method/main_result.py
@@ -0,0 +1,15 @@
+from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal
+from experiment.paper_past_snow_loads.method.study_visualizer_for_non_stationary_trends import \
+ StudyVisualizerForNonStationaryTrends
+
+ALTITUDE_TRENDS = [1800]
+
+
+def draw_snow_load_map(altitude):
+ study = CrocusSnowLoadTotal(altitude=altitude)
+ visualizer = StudyVisualizerForNonStationaryTrends(study, multiprocessing=True)
+ visualizer.plot_trends()
+
+if __name__ == '__main__':
+ draw_snow_load_map(altitude=1800)
+
diff --git a/experiment/paper_past_snow_loads/method/study_visualizer_for_non_stationary_trends.py b/experiment/paper_past_snow_loads/method/study_visualizer_for_non_stationary_trends.py
index e69de29b..be8f34ad 100644
--- a/experiment/paper_past_snow_loads/method/study_visualizer_for_non_stationary_trends.py
+++ b/experiment/paper_past_snow_loads/method/study_visualizer_for_non_stationary_trends.py
@@ -0,0 +1,90 @@
+from typing import Dict
+
+import numpy as np
+from cached_property import cached_property
+
+from experiment.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
+from experiment.meteo_france_data.scm_models_data.visualization.study_visualization.study_visualizer import \
+ StudyVisualizer
+from experiment.trend_analysis.abstract_score import MeanScore
+from experiment.trend_analysis.univariate_test.abstract_gev_trend_test import AbstractGevTrendTest
+from experiment.trend_analysis.univariate_test.gev_trend_test_one_parameter import GevScaleTrendTest, \
+ GevLocationTrendTest
+from experiment.trend_analysis.univariate_test.gev_trend_test_two_parameters import GevLocationAndScaleTrendTest
+
+
+class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
+
+ def __init__(self, study: AbstractStudy, show=True, save_to_file=False, only_one_graph=False, only_first_row=False,
+ vertical_kde_plot=False, year_for_kde_plot=None, plot_block_maxima_quantiles=False,
+ temporal_non_stationarity=False, transformation_class=None, verbose=False, multiprocessing=False,
+ complete_non_stationary_trend_analysis=False, normalization_under_one_observations=True,
+ score_class=MeanScore):
+ super().__init__(study, show, save_to_file, only_one_graph, only_first_row, vertical_kde_plot,
+ year_for_kde_plot, plot_block_maxima_quantiles, temporal_non_stationarity,
+ transformation_class, verbose, multiprocessing, complete_non_stationary_trend_analysis,
+ normalization_under_one_observations, score_class)
+ self.non_stationary_trend_test = [GevLocationTrendTest, GevScaleTrendTest, GevLocationAndScaleTrendTest]
+
+ # Utils
+
+ @cached_property
+ def massif_name_to_years_and_maxima(self):
+ d = {}
+ df_maxima = self.study.observations_annual_maxima.df_maxima_gev
+ years = np.array(df_maxima.columns)
+ for massif_name, s_maxima in df_maxima.iterrows():
+ d[massif_name] = (years, np.array(s_maxima))
+ return d
+
+ @cached_property
+ def massif_name_to_psnow(self):
+ return {m: np.count_nonzero(maxima) / len(maxima) for m, (_, maxima) in
+ self.massif_name_to_years_and_maxima.items()}
+
+ @cached_property
+ def massif_name_to_non_null_years_and_maxima(self):
+ d = {}
+ for m, (years, maxima) in self.massif_name_to_years_and_maxima.items():
+ mask = np.nonzero(maxima)
+ d[m] = (years[mask], maxima[mask])
+ return d
+
+ @cached_property
+ def massif_name_to_minimized_aic_non_stationary_trend_test(self) -> Dict[str, AbstractGevTrendTest]:
+ starting_year = 1958
+ massif_name_to_trend_test_that_minimized_aic = {}
+ for massif_name, (x, y) in self.massif_name_to_non_null_years_and_maxima.items():
+ non_stationary_trend_test = [t(x, y, starting_year) for t in self.non_stationary_trend_test]
+ trend_test_that_minimized_aic = sorted(non_stationary_trend_test, key=lambda t: t.aic)[0]
+ massif_name_to_trend_test_that_minimized_aic[massif_name] = trend_test_that_minimized_aic
+ return massif_name_to_trend_test_that_minimized_aic
+
+ # Part 1 - Trends
+
+ def plot_trends(self):
+ v = max(abs(min(self.massif_name_to_tdrl_value.values())), max(self.massif_name_to_tdrl_value.values()))
+ vmin, vmax = -v, v
+ self.study.visualize_study(massif_name_to_value=self.massif_name_to_tdrl_value, vmin=vmin, vmax=vmax,
+ replace_blue_by_white=False, axis_off=True, show_label=False,
+ add_colorbar=True)
+
+ @cached_property
+ def massif_name_to_tdrl_value(self):
+ return {m: t.time_derivative_of_return_level for m, t in self.massif_name_to_minimized_aic_non_stationary_trend_test.items()}
+
+ @property
+ def massif_name_to_minimized_aic_model_class(self):
+ return {m: t.unconstrained_model_class for m, t in self.massif_name_to_minimized_aic_non_stationary_trend_test.items()}
+
+ @property
+ def massif_name_to_model_significance_symbol(self):
+ return {}
+
+
+ # Part 1 - Uncertainty return level plot
+
+
+
+ def massif_name_to_uncertainty(self):
+ pass
diff --git a/experiment/trend_analysis/univariate_test/abstract_gev_trend_test.py b/experiment/trend_analysis/univariate_test/abstract_gev_trend_test.py
index 96096bbd..bfe20929 100644
--- a/experiment/trend_analysis/univariate_test/abstract_gev_trend_test.py
+++ b/experiment/trend_analysis/univariate_test/abstract_gev_trend_test.py
@@ -27,15 +27,29 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
constrained_model_class=StationaryTemporalModel,
):
super().__init__(years, maxima, starting_year)
- self.fit_method = TemporalMarginFitMethod.is_mev_gev_fit
+ self.unconstrained_model_class = unconstrained_model_class
+ self.constrained_model_class = constrained_model_class
+ self.fit_method = TemporalMarginFitMethod.extremes_fevd_mle
# Load observations, coordinates and datasets
self.coordinates, self.dataset = load_temporal_coordinates_and_dataset(maxima, years)
+ # By default crashed boolean is False
+ self.crashed = False
try:
- # Fit constrained model
- self.constrained_estimator = fitted_linear_margin_estimator(constrained_model_class, self.coordinates, self.dataset, self.starting_year, self.fit_method)
- # Fit unconstrained model
- self.unconstrained_estimator = fitted_linear_margin_estimator(unconstrained_model_class, self.coordinates, self.dataset, self.starting_year, self.fit_method)
- self.crashed = False
+ pass
+ except SafeRunException:
+ self.crashed = True
+
+ @cached_property
+ def constrained_estimator(self):
+ try:
+ return fitted_linear_margin_estimator(self.constrained_model_class, self.coordinates, self.dataset, self.starting_year, self.fit_method)
+ except SafeRunException:
+ self.crashed = True
+
+ @cached_property
+ def unconstrained_estimator(self):
+ try:
+ return fitted_linear_margin_estimator(self.unconstrained_model_class, self.coordinates, self.dataset, self.starting_year, self.fit_method)
except SafeRunException:
self.crashed = True
@@ -69,6 +83,15 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
def degree_freedom_chi2(self) -> int:
raise NotImplementedError
+ @property
+ def total_number_of_parameters_for_unconstrained_model(self) -> int:
+ return self.degree_freedom_chi2 + 3
+
+ @property
+ def aic(self):
+ # deviance = - 2 * nllh
+ return 2 * self.total_number_of_parameters_for_unconstrained_model - self.unconstrained_model_deviance
+
@property
def likelihood_ratio(self):
return self.unconstrained_model_deviance - self.constrained_model_deviance
@@ -82,23 +105,25 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
@property
def unconstrained_model_deviance(self):
+ unconstrained_estimator = self.unconstrained_estimator
if self.crashed:
return np.nan
else:
- return self.unconstrained_estimator.result_from_model_fit.deviance
+ return unconstrained_estimator.result_from_model_fit.deviance
@property
def unconstained_nllh(self):
+ unconstrained_estimator = self.unconstrained_estimator
if self.crashed:
return np.nan
else:
- return self.unconstrained_estimator.result_from_model_fit.nllh
+ return unconstrained_estimator.result_from_model_fit.nllh
# Evolution of the GEV parameters and corresponding quantiles
@property
def test_sign(self) -> int:
- return np.sign(self.test_trend_slope_strength)
+ return np.sign(self.time_derivative_of_return_level)
def get_non_stationary_linear_coef(self, gev_param_name: str):
return self.unconstrained_estimator.margin_function_from_fit.get_coef(gev_param_name,
@@ -117,7 +142,7 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
is_transformed=False)
@property
- def test_trend_slope_strength(self):
+ def time_derivative_of_return_level(self):
if self.crashed:
return 0.0
else:
diff --git a/experiment/trend_analysis/univariate_test/abstract_univariate_test.py b/experiment/trend_analysis/univariate_test/abstract_univariate_test.py
index a039b9d9..cf3f9b5c 100644
--- a/experiment/trend_analysis/univariate_test/abstract_univariate_test.py
+++ b/experiment/trend_analysis/univariate_test/abstract_univariate_test.py
@@ -105,7 +105,7 @@ class AbstractUnivariateTest(object):
return len(self.years)
@property
- def test_trend_slope_strength(self):
+ def time_derivative_of_return_level(self):
return 0.0
@property
diff --git a/experiment/trend_analysis/univariate_test/univariate_test_results.py b/experiment/trend_analysis/univariate_test/univariate_test_results.py
index 3b3c0695..009da783 100644
--- a/experiment/trend_analysis/univariate_test/univariate_test_results.py
+++ b/experiment/trend_analysis/univariate_test/univariate_test_results.py
@@ -12,7 +12,7 @@ def compute_gev_change_point_test_result(smooth_maxima, starting_year, trend_tes
trend_test = trend_test_class(years, smooth_maxima, starting_year) # type: AbstractGevTrendTest
assert isinstance(trend_test, AbstractGevTrendTest)
return trend_test.test_trend_type, \
- trend_test.test_trend_slope_strength, \
+ trend_test.time_derivative_of_return_level, \
trend_test.unconstained_nllh, \
trend_test.test_trend_constant_quantile, \
trend_test.mean_difference_same_sign_as_slope_strenght, \
diff --git a/extreme_fit/model/result_from_model_fit/abstract_result_from_model_fit.py b/extreme_fit/model/result_from_model_fit/abstract_result_from_model_fit.py
index be61be0d..c82898e0 100644
--- a/extreme_fit/model/result_from_model_fit/abstract_result_from_model_fit.py
+++ b/extreme_fit/model/result_from_model_fit/abstract_result_from_model_fit.py
@@ -42,7 +42,7 @@ class AbstractResultFromModelFit(object):
@property
def deviance(self):
- raise NotImplementedError
+ return - 2 * self.nllh
@property
def convergence(self) -> str:
diff --git a/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_result_from_extremes.py b/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_result_from_extremes.py
index cbd16d9a..acf56fa5 100644
--- a/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_result_from_extremes.py
+++ b/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_result_from_extremes.py
@@ -14,6 +14,14 @@ class AbstractResultFromExtremes(AbstractResultFromModelFit):
super().__init__(result_from_fit)
self.gev_param_name_to_dim = gev_param_name_to_dim
+ @property
+ def results(self):
+ return self.get_python_dictionary(self.name_to_value['results'])
+
+ @property
+ def nllh(self):
+ return np.array(self.results['value'])[0]
+
@property
def is_non_stationary(self):
return len(self.gev_param_name_to_dim) > 0
diff --git a/extreme_fit/model/result_from_model_fit/result_from_ismev.py b/extreme_fit/model/result_from_model_fit/result_from_ismev.py
index 2332a660..22eef48d 100644
--- a/extreme_fit/model/result_from_model_fit/result_from_ismev.py
+++ b/extreme_fit/model/result_from_model_fit/result_from_ismev.py
@@ -26,10 +26,6 @@ class ResultFromIsmev(AbstractResultFromModelFit):
def nllh(self):
return convertFloatVector_to_float(self.name_to_value['nllh'])
- @property
- def deviance(self):
- return - 2 * self.nllh
-
@property
def convergence(self) -> str:
return convertFloatVector_to_float(self.name_to_value['conv']) == 0
--
GitLab