From 1a6c3690856924390c60360fac32487e9dcaf475 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Mon, 16 Dec 2019 12:28:20 +0100
Subject: [PATCH] [PAPER 1] add gumbel type for mle fit. add
study_visualizer_for_shape_repartition both for gumbel VS stationary and
stationary VS non-stationary
---
.../study_visualization/study_visualizer.py | 7 +-
.../main_shape_repartition.py | 30 ++++++++
.../study_visualizer_for_shape_repartition.py | 61 ++++++++++++++++
...main_comparison_with_eurocode_examples.py} | 0
.../main_comparison_with_eurocode_global.py | 70 +++++++++++++++++++
.../paper_past_snow_loads/paper_main_utils.py | 8 ++-
...dy_visualizer_for_non_stationary_trends.py | 16 +++--
.../abstract_gev_trend_test.py | 6 +-
.../gev_trend_test_one_parameter.py | 44 ++++++++----
.../gev_trend_test_two_parameters.py | 4 +-
.../abstract_temporal_linear_margin_model.py | 10 ++-
.../temporal_linear_margin_models.py | 14 +++-
.../result_from_mle_extremes.py | 8 ++-
.../model/result_from_model_fit/utils.py | 10 ++-
.../test_gev_temporal_extremes_gumbel.py | 51 ++++++++++++++
15 files changed, 302 insertions(+), 37 deletions(-)
create mode 100644 experiment/paper_past_snow_loads/check_mle_convergence_for_trends/main_shape_repartition.py
create mode 100644 experiment/paper_past_snow_loads/check_mle_convergence_for_trends/study_visualizer_for_shape_repartition.py
rename experiment/paper_past_snow_loads/discussion_data_comparison_with_eurocode/{main_comparison_with_eurocode.py => main_comparison_with_eurocode_examples.py} (100%)
create mode 100644 experiment/paper_past_snow_loads/discussion_data_comparison_with_eurocode/main_comparison_with_eurocode_global.py
create mode 100644 test/test_extreme_fit/test_distribution/test_gev/test_gev_temporal_extremes_gumbel.py
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 132f14bb..d6da42c2 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
@@ -632,7 +632,10 @@ class StudyVisualizer(VisualizationParameters):
# Display the graph of the max on top
if ax is None:
ax = plt.gca()
- x, y = self.smooth_maxima_x_y(massif_names.index(massif_name))
+ if massif_name is None:
+ x, y = self.study.ordered_years, self.study.observations_annual_maxima.df_maxima_gev.mean(axis=0)
+ else:
+ x, y = self.smooth_maxima_x_y(massif_names.index(massif_name))
ax.plot(x, y, color=color, linewidth=5, label=label, linestyle=linestyle)
# ax.set_ylabel('{} (in {})'.format(snow_abbreviation, self.study.variable_unit), color=color, fontsize=15)
@@ -642,7 +645,7 @@ class StudyVisualizer(VisualizationParameters):
self.plot_name = plot_name
ax.set_ylabel('{} ({})'.format(snow_abbreviation, self.study.variable_unit), fontsize=15)
ax.set_xlabel('years', fontsize=15)
- if label is None:
+ if label is None and massif_name is not None:
ax.set_title('{} at {} m'.format(massif_name, altitude))
if last_plot:
ax.legend()
diff --git a/experiment/paper_past_snow_loads/check_mle_convergence_for_trends/main_shape_repartition.py b/experiment/paper_past_snow_loads/check_mle_convergence_for_trends/main_shape_repartition.py
new file mode 100644
index 00000000..26988557
--- /dev/null
+++ b/experiment/paper_past_snow_loads/check_mle_convergence_for_trends/main_shape_repartition.py
@@ -0,0 +1,30 @@
+from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal
+from experiment.paper_past_snow_loads.check_mle_convergence_for_trends.study_visualizer_for_shape_repartition import \
+ StudyVisualizerForShape, StudyVisualizerGumbel
+from experiment.paper_past_snow_loads.paper_main_utils import load_altitude_to_visualizer
+
+
+def main_shape_repartition(altitudes, massif_names=None,
+ non_stationary_uncertainty=None, uncertainty_methods=None,
+ study_class=CrocusSnowLoadTotal,
+ study_visualizer_class=StudyVisualizerForShape, save_to_file=True):
+ # Load altitude to visualizer
+ altitude_to_visualizer = load_altitude_to_visualizer(altitudes, massif_names, non_stationary_uncertainty,
+ study_class, uncertainty_methods,
+ study_visualizer_class=study_visualizer_class,
+ save_to_file=save_to_file)
+ altitudes_for_plot_trend = altitudes #[900, 1800, 2700]
+ visualizers_for_altitudes = [visualizer
+ for altitude, visualizer in altitude_to_visualizer.items()
+ if altitude in altitudes_for_plot_trend]
+ max_abs_tdrl = max([visualizer.max_abs_change for visualizer in visualizers_for_altitudes])
+ for visualizer in visualizers_for_altitudes:
+ # visualizer.plot_trends(max_abs_tdrl, add_colorbar=visualizer.study.altitude == 2700)
+ visualizer.plot_trends(max_abs_tdrl, add_colorbar=True)
+
+
+if __name__ == '__main__':
+ # main_shape_repartition([900], save_to_file=False)
+ # main_shape_repartition([900, 1800, 2700])
+ # main_shape_repartition([300, 600, 900, 1200, 1500, 1800, 2700])
+ main_shape_repartition([900], study_visualizer_class=StudyVisualizerGumbel)
diff --git a/experiment/paper_past_snow_loads/check_mle_convergence_for_trends/study_visualizer_for_shape_repartition.py b/experiment/paper_past_snow_loads/check_mle_convergence_for_trends/study_visualizer_for_shape_repartition.py
new file mode 100644
index 00000000..d7b62f8e
--- /dev/null
+++ b/experiment/paper_past_snow_loads/check_mle_convergence_for_trends/study_visualizer_for_shape_repartition.py
@@ -0,0 +1,61 @@
+from cached_property import cached_property
+
+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
+from experiment.trend_analysis.univariate_test.gev_trend_test_one_parameter import GevLocationTrendTest, \
+ GevStationaryVersusGumbel
+
+
+class StudyVisualizerForShape(StudyVisualizerForNonStationaryTrends):
+
+ 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, uncertainty_methods=None, non_stationary_contexts=None,
+ uncertainty_massif_names=None, effective_temporal_covariate=2017, relative_change_trend_plot=True):
+ 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, uncertainty_methods,
+ non_stationary_contexts, uncertainty_massif_names, effective_temporal_covariate,
+ relative_change_trend_plot)
+
+ @cached_property
+ def massif_name_to_unconstrained_shape_parameter(self):
+ return {m: t.unconstrained_estimator_gev_params.shape
+ for m, t in self.massif_name_to_minimized_aic_non_stationary_trend_test.items()}
+
+ @cached_property
+ def massif_name_to_change_value(self):
+ print(self.massif_name_to_unconstrained_shape_parameter)
+ return self.massif_name_to_unconstrained_shape_parameter
+
+ @property
+ def label(self):
+ return 'Shape parameter value'
+
+ @property
+ def graduation(self):
+ return 0.1
+
+
+class StudyVisualizerGumbel(StudyVisualizerForNonStationaryTrends):
+
+ 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, uncertainty_methods=None, non_stationary_contexts=None,
+ uncertainty_massif_names=None, effective_temporal_covariate=2017, relative_change_trend_plot=True):
+ 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, uncertainty_methods,
+ non_stationary_contexts, uncertainty_massif_names, effective_temporal_covariate,
+ relative_change_trend_plot)
+ # Assign default argument for the non stationary trends
+ self.non_stationary_trend_test = [GevStationaryVersusGumbel]
+ self.non_stationary_trend_test_to_marker = dict(zip(self.non_stationary_trend_test, ["o"]))
diff --git a/experiment/paper_past_snow_loads/discussion_data_comparison_with_eurocode/main_comparison_with_eurocode.py b/experiment/paper_past_snow_loads/discussion_data_comparison_with_eurocode/main_comparison_with_eurocode_examples.py
similarity index 100%
rename from experiment/paper_past_snow_loads/discussion_data_comparison_with_eurocode/main_comparison_with_eurocode.py
rename to experiment/paper_past_snow_loads/discussion_data_comparison_with_eurocode/main_comparison_with_eurocode_examples.py
diff --git a/experiment/paper_past_snow_loads/discussion_data_comparison_with_eurocode/main_comparison_with_eurocode_global.py b/experiment/paper_past_snow_loads/discussion_data_comparison_with_eurocode/main_comparison_with_eurocode_global.py
new file mode 100644
index 00000000..22b5593b
--- /dev/null
+++ b/experiment/paper_past_snow_loads/discussion_data_comparison_with_eurocode/main_comparison_with_eurocode_global.py
@@ -0,0 +1,70 @@
+from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusDepth
+from experiment.meteo_france_data.scm_models_data.crocus.crocus_variables import CrocusDepthVariable
+from experiment.meteo_france_data.scm_models_data.visualization.study_visualization.main_study_visualizer import \
+ study_iterator_global, SCM_STUDY_CLASS_TO_ABBREVIATION, snow_density_str
+from experiment.meteo_france_data.scm_models_data.visualization.study_visualization.study_visualizer import \
+ StudyVisualizer
+import matplotlib.pyplot as plt
+
+from experiment.paper_past_snow_loads.discussion_data_comparison_with_eurocode.crocus_study_comparison_with_eurocode import \
+ CrocusDifferenceSnowLoad, \
+ CrocusSnowDensityAtMaxofSwe, CrocusDifferenceSnowLoadRescaledAndEurocodeToSeeSynchronization, \
+ CrocusSnowDepthAtMaxofSwe, CrocusSnowDepthDifference
+from experiment.paper_past_snow_loads.paper_utils import dpi_paper1_figure
+
+
+def max_graph_annual_maxima_comparison():
+ """
+ We choose these massif because each represents a different eurocode region
+ we also choose them because they belong to a different climatic area
+ :return:
+ """
+ save_to_file = False
+ study_classes = [
+ CrocusSnowDensityAtMaxofSwe,
+ CrocusDifferenceSnowLoad,
+ CrocusSnowDepthDifference,
+ ][:]
+ study_class_to_ylim_and_yticks = {
+ CrocusSnowDensityAtMaxofSwe: ([150, 500], [50 * i for i in range(3, 11)]),
+ CrocusDifferenceSnowLoad: ([0, 10], [2 * i for i in range(0, 6)]),
+ CrocusSnowDepthDifference: ([0, 0.6], [0.2 * i for i in range(0, 4)]),
+ }
+ for study_class in study_classes:
+ ylim, yticks = study_class_to_ylim_and_yticks[study_class]
+
+ marker_altitude = [
+ ('deepskyblue', 900),
+ ('dodgerblue', 1800),
+ ('blue', 2700),
+ ][:]
+ ax = plt.gca()
+ for color, altitude in marker_altitude:
+ for study in study_iterator_global([study_class], altitudes=[altitude]):
+ study_visualizer = StudyVisualizer(study, save_to_file=save_to_file,
+ verbose=True,
+ multiprocessing=True)
+ snow_abbreviation = SCM_STUDY_CLASS_TO_ABBREVIATION[study_class]
+ nb_massifs = len(study_visualizer.study.study_massif_names)
+ label = 'Mean value at {}m (out of {} massifs)'.format(altitude, nb_massifs)
+ study_visualizer.visualize_max_graphs_poster(None, altitude, snow_abbreviation, color, label,
+ False, ax)
+
+ last_plot = altitude == 2700
+ if last_plot:
+ # if study_class == CrocusSnowDensityAtMaxofSwe:
+ # label = '{} for French standards'.format(snow_density_str)
+ # snow_density_eurocode = [150 for _ in study.ordered_years]
+ # ax.plot(study.ordered_years, snow_density_eurocode, color='k', label=label)
+ ax.legend()
+ tight_pad = {'h_pad': 0.2}
+ ax.set_ylim(ylim)
+ ax.set_xlim([1957, 2018])
+ ax.yaxis.set_ticks(yticks)
+ study_visualizer.show_or_save_to_file(no_title=True, tight_layout=True,
+ tight_pad=tight_pad, dpi=dpi_paper1_figure)
+ ax.clear()
+
+
+if __name__ == '__main__':
+ max_graph_annual_maxima_comparison()
diff --git a/experiment/paper_past_snow_loads/paper_main_utils.py b/experiment/paper_past_snow_loads/paper_main_utils.py
index 94493f2f..510dcdec 100644
--- a/experiment/paper_past_snow_loads/paper_main_utils.py
+++ b/experiment/paper_past_snow_loads/paper_main_utils.py
@@ -4,11 +4,13 @@ from experiment.paper_past_snow_loads.study_visualizer_for_non_stationary_trends
StudyVisualizerForNonStationaryTrends
-def load_altitude_to_visualizer(altitudes, massif_names, non_stationary_uncertainty, study_class, uncertainty_methods):
+def load_altitude_to_visualizer(altitudes, massif_names, non_stationary_uncertainty, study_class, uncertainty_methods,
+ study_visualizer_class=StudyVisualizerForNonStationaryTrends,
+ save_to_file=True):
altitude_to_visualizer = OrderedDict()
for altitude in altitudes:
- altitude_to_visualizer[altitude] = StudyVisualizerForNonStationaryTrends(
- study=study_class(altitude=altitude), multiprocessing=True, save_to_file=True,
+ altitude_to_visualizer[altitude] = study_visualizer_class(
+ study=study_class(altitude=altitude), multiprocessing=True, save_to_file=save_to_file,
uncertainty_massif_names=massif_names, uncertainty_methods=uncertainty_methods,
non_stationary_contexts=non_stationary_uncertainty)
return altitude_to_visualizer
diff --git a/experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py b/experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py
index 4fa5978d..976efe49 100644
--- a/experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py
+++ b/experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py
@@ -133,7 +133,6 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
ax.get_xaxis().set_visible(True)
ax.set_xticks([])
ax.set_xlabel('Altitude = {}m'.format(self.study.altitude), fontsize=15)
-
self.plot_name = 'tdlr_trends'
self.show_or_save_to_file(add_classic_title=False, tight_layout=True, no_title=True,
dpi=500)
@@ -157,16 +156,20 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
@property
def ticks_values_and_labels(self):
- graduation = 10 if self.relative_change_trend_plot else 0.2
positive_ticks = []
- tick = graduation
+ tick = self.graduation
while tick < self._max_abs_change:
positive_ticks.append(round(tick, 1))
- tick += graduation
+ tick += self.graduation
all_ticks_labels = [-t for t in positive_ticks] + [0] + positive_ticks
ticks_values = [((t / self._max_abs_change) + 1) / 2 for t in all_ticks_labels]
return ticks_values, all_ticks_labels
+ @property
+ def graduation(self):
+ graduation = 10 if self.relative_change_trend_plot else 0.2
+ return graduation
+
@cached_property
def massif_name_to_tdrl_value(self):
return {m: t.time_derivative_of_return_level for m, t in
@@ -291,3 +294,8 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
uncertainty.confidence_interval[1] > eurocode)
for eurocode, uncertainty in eurocode_and_uncertainties])
return 100 * np.mean(a, axis=0)
+
+
+
+
+
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 71ff1725..8d3b47b0 100644
--- a/experiment/trend_analysis/univariate_test/abstract_gev_trend_test.py
+++ b/experiment/trend_analysis/univariate_test/abstract_gev_trend_test.py
@@ -132,13 +132,13 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
AbstractCoordinates.COORDINATE_T)
@cached_property
- def non_stationary_constant_gev_params(self) -> GevParams:
+ def unconstrained_estimator_gev_params(self) -> GevParams:
# Constant parameters correspond to the gev params in 1958
return self.unconstrained_estimator.margin_function_from_fit.get_gev_params(coordinate=np.array([1958]),
is_transformed=False)
@cached_property
- def stationary_constant_gev_params(self) -> GevParams:
+ def constrained_estimator_gev_params(self) -> GevParams:
# Constant parameters correspond to any gev params
return self.constrained_estimator.margin_function_from_fit.get_gev_params(coordinate=np.array([1958]),
is_transformed=False)
@@ -192,4 +192,4 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
if self.crashed:
return 0.0
else:
- return self.non_stationary_constant_gev_params.quantile(p=self.quantile_level)
+ return self.unconstrained_estimator_gev_params.quantile(p=self.quantile_level)
diff --git a/experiment/trend_analysis/univariate_test/gev_trend_test_one_parameter.py b/experiment/trend_analysis/univariate_test/gev_trend_test_one_parameter.py
index 8af03890..5a89e207 100644
--- a/experiment/trend_analysis/univariate_test/gev_trend_test_one_parameter.py
+++ b/experiment/trend_analysis/univariate_test/gev_trend_test_one_parameter.py
@@ -1,28 +1,35 @@
from experiment.eurocode_data.utils import EUROCODE_QUANTILE
from experiment.trend_analysis.univariate_test.abstract_gev_trend_test import AbstractGevTrendTest
from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import \
- NonStationaryLocationTemporalModel, NonStationaryScaleTemporalModel, NonStationaryShapeTemporalModel
+ NonStationaryLocationTemporalModel, NonStationaryScaleTemporalModel, NonStationaryShapeTemporalModel, \
+ StationaryTemporalModel, GumbelTemporalModel
from extreme_fit.distribution.gev.gev_params import GevParams
class GevTrendTestOneParameter(AbstractGevTrendTest):
- def __init__(self, years, maxima, starting_year, unconstrained_model_class, gev_param_name, quantile_level=EUROCODE_QUANTILE):
- super().__init__(years, maxima, starting_year,
- unconstrained_model_class=unconstrained_model_class,
- quantile_level=quantile_level)
- self.gev_param_name = gev_param_name
-
@property
def degree_freedom_chi2(self) -> int:
return 1
+
+class GevTrendTestOneParameterAgainstStationary(GevTrendTestOneParameter):
+
+ def __init__(self, years, maxima, starting_year, unconstrained_model_class, gev_param_name,
+ quantile_level=EUROCODE_QUANTILE,
+ constrained_model_class=StationaryTemporalModel):
+ super().__init__(years, maxima, starting_year,
+ unconstrained_model_class=unconstrained_model_class,
+ quantile_level=quantile_level,
+ constrained_model_class=constrained_model_class)
+ self.gev_param_name = gev_param_name
+
@property
def non_stationary_linear_coef(self):
return self.get_non_stationary_linear_coef(gev_param_name=self.gev_param_name)
-class GevLocationTrendTest(GevTrendTestOneParameter):
+class GevLocationTrendTest(GevTrendTestOneParameterAgainstStationary):
def __init__(self, years, maxima, starting_year, quantile_level=EUROCODE_QUANTILE):
super().__init__(years, maxima, starting_year,
@@ -31,12 +38,12 @@ class GevLocationTrendTest(GevTrendTestOneParameter):
quantile_level=quantile_level)
def _slope_strength(self):
- return self.non_stationary_constant_gev_params.time_derivative_of_return_level(p=self.quantile_level,
+ return self.unconstrained_estimator_gev_params.time_derivative_of_return_level(p=self.quantile_level,
mu1=self.non_stationary_linear_coef)
@property
def mean_difference_same_sign_as_slope_strenght(self) -> bool:
- zeta0 = self.non_stationary_constant_gev_params.shape
+ zeta0 = self.unconstrained_estimator_gev_params.shape
mean_difference = self.mean_difference(zeta0=zeta0, mu1=self.non_stationary_linear_coef)
return self.same_sign(mean_difference, self._slope_strength())
@@ -45,7 +52,7 @@ class GevLocationTrendTest(GevTrendTestOneParameter):
return False
-class GevScaleTrendTest(GevTrendTestOneParameter):
+class GevScaleTrendTest(GevTrendTestOneParameterAgainstStationary):
def __init__(self, years, maxima, starting_year, quantile_level=EUROCODE_QUANTILE):
super().__init__(years, maxima, starting_year,
@@ -54,13 +61,13 @@ class GevScaleTrendTest(GevTrendTestOneParameter):
quantile_level=quantile_level)
def _slope_strength(self):
- return self.non_stationary_constant_gev_params.time_derivative_of_return_level(
+ return self.unconstrained_estimator_gev_params.time_derivative_of_return_level(
p=self.quantile_level,
sigma1=self.non_stationary_linear_coef)
@property
def mean_difference_same_sign_as_slope_strenght(self) -> bool:
- zeta0 = self.non_stationary_constant_gev_params.shape
+ zeta0 = self.unconstrained_estimator_gev_params.shape
mean_difference = self.mean_difference(zeta0=zeta0, sigma1=self.non_stationary_linear_coef)
return self.same_sign(mean_difference, self._slope_strength())
@@ -70,10 +77,19 @@ class GevScaleTrendTest(GevTrendTestOneParameter):
return self.same_sign(sigma1, self._slope_strength())
-class GevShapeTrendTest(GevTrendTestOneParameter):
+class GevShapeTrendTest(GevTrendTestOneParameterAgainstStationary):
def __init__(self, years, maxima, starting_year, quantile_level=EUROCODE_QUANTILE):
super().__init__(years, maxima, starting_year,
unconstrained_model_class=NonStationaryShapeTemporalModel,
gev_param_name=GevParams.SHAPE,
quantile_level=quantile_level)
+
+
+class GevStationaryVersusGumbel(GevTrendTestOneParameter):
+
+ def __init__(self, years, maxima, starting_year, quantile_level=EUROCODE_QUANTILE):
+ super().__init__(years, maxima, starting_year,
+ unconstrained_model_class=StationaryTemporalModel,
+ constrained_model_class=GumbelTemporalModel,
+ quantile_level=quantile_level)
diff --git a/experiment/trend_analysis/univariate_test/gev_trend_test_two_parameters.py b/experiment/trend_analysis/univariate_test/gev_trend_test_two_parameters.py
index 95525524..2bd8a106 100644
--- a/experiment/trend_analysis/univariate_test/gev_trend_test_two_parameters.py
+++ b/experiment/trend_analysis/univariate_test/gev_trend_test_two_parameters.py
@@ -29,13 +29,13 @@ class GevLocationAndScaleTrendTest(GevTrendTestTwoParameters):
return self.get_non_stationary_linear_coef(gev_param_name=GevParams.SCALE)
def _slope_strength(self):
- return self.non_stationary_constant_gev_params.time_derivative_of_return_level(p=self.quantile_level,
+ return self.unconstrained_estimator_gev_params.time_derivative_of_return_level(p=self.quantile_level,
mu1=self.mu1,
sigma1=self.sigma1)
@property
def mean_difference_same_sign_as_slope_strenght(self) -> bool:
- zeta0 = self.non_stationary_constant_gev_params.shape
+ zeta0 = self.unconstrained_estimator_gev_params.shape
mean_difference = self.mean_difference(zeta0=zeta0, mu1=self.mu1, sigma1=self.sigma1)
return self.same_sign(mean_difference, self._slope_strength())
diff --git a/extreme_fit/model/margin_model/linear_margin_model/abstract_temporal_linear_margin_model.py b/extreme_fit/model/margin_model/linear_margin_model/abstract_temporal_linear_margin_model.py
index 3dd02216..08700dfa 100644
--- a/extreme_fit/model/margin_model/linear_margin_model/abstract_temporal_linear_margin_model.py
+++ b/extreme_fit/model/margin_model/linear_margin_model/abstract_temporal_linear_margin_model.py
@@ -25,8 +25,10 @@ class AbstractTemporalLinearMarginModel(LinearMarginModel):
def __init__(self, coordinates: AbstractCoordinates, use_start_value=False, params_start_fit=None,
params_sample=None, starting_point=None, fit_method=TemporalMarginFitMethod.is_mev_gev_fit,
- nb_iterations_for_bayesian_fit=5000, params_start_fit_bayesian=None):
+ nb_iterations_for_bayesian_fit=5000, params_start_fit_bayesian=None,
+ type_for_MLE="GEV"):
super().__init__(coordinates, use_start_value, params_start_fit, params_sample, starting_point)
+ self.type_for_mle = type_for_MLE
self.params_start_fit_bayesian = params_start_fit_bayesian
self.nb_iterations_for_bayesian_fit = nb_iterations_for_bayesian_fit
assert isinstance(fit_method, TemporalMarginFitMethod)
@@ -59,10 +61,12 @@ class AbstractTemporalLinearMarginModel(LinearMarginModel):
res = safe_run_r_estimator(function=r('fevd_fixed'),
x=x,
data=y,
- method='MLE',
+ type=self.type_for_mle,
+ method="MLE",
**r_type_argument_kwargs
)
- return ResultFromMleExtremes(res, self.margin_function_start_fit.gev_param_name_to_dims)
+ return ResultFromMleExtremes(res, self.margin_function_start_fit.gev_param_name_to_dims,
+ type_for_mle=self.type_for_mle)
def extremes_fevd_bayesian_fit(self, x, df_coordinates_temp) -> AbstractResultFromExtremes:
r_type_argument_kwargs, y = self.extreme_arguments(df_coordinates_temp)
diff --git a/extreme_fit/model/margin_model/linear_margin_model/temporal_linear_margin_models.py b/extreme_fit/model/margin_model/linear_margin_model/temporal_linear_margin_models.py
index 33d267bc..e6a54751 100644
--- a/extreme_fit/model/margin_model/linear_margin_model/temporal_linear_margin_models.py
+++ b/extreme_fit/model/margin_model/linear_margin_model/temporal_linear_margin_models.py
@@ -1,7 +1,8 @@
from extreme_fit.model.margin_model.linear_margin_model.abstract_temporal_linear_margin_model import \
- AbstractTemporalLinearMarginModel
+ AbstractTemporalLinearMarginModel, TemporalMarginFitMethod
from extreme_fit.model.utils import r
from extreme_fit.distribution.gev.gev_params import GevParams
+from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
class StationaryTemporalModel(AbstractTemporalLinearMarginModel):
@@ -10,6 +11,15 @@ class StationaryTemporalModel(AbstractTemporalLinearMarginModel):
super().load_margin_functions({})
+class GumbelTemporalModel(StationaryTemporalModel):
+
+ def __init__(self, coordinates: AbstractCoordinates, use_start_value=False, params_start_fit=None,
+ params_sample=None, starting_point=None, fit_method=TemporalMarginFitMethod.is_mev_gev_fit,
+ nb_iterations_for_bayesian_fit=5000, params_start_fit_bayesian=None):
+ super().__init__(coordinates, use_start_value, params_start_fit, params_sample, starting_point, fit_method,
+ nb_iterations_for_bayesian_fit, params_start_fit_bayesian, type_for_MLE="Gumbel")
+
+
class NonStationaryLocationTemporalModel(AbstractTemporalLinearMarginModel):
def load_margin_functions(self, gev_param_name_to_dims=None):
@@ -62,4 +72,4 @@ class NonStationaryLocationAndScaleTemporalModel(AbstractTemporalLinearMarginMod
@property
def sigl(self):
- return 1
\ No newline at end of file
+ return 1
diff --git a/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_mle_extremes.py b/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_mle_extremes.py
index 56a167bf..9fa9b18f 100644
--- a/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_mle_extremes.py
+++ b/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_mle_extremes.py
@@ -1,5 +1,6 @@
import numpy as np
import rpy2
+from rpy2 import robjects
from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_result_from_extremes import \
AbstractResultFromExtremes
@@ -11,12 +12,17 @@ from extreme_fit.model.utils import r
class ResultFromMleExtremes(AbstractResultFromExtremes):
+ def __init__(self, result_from_fit: robjects.ListVector, gev_param_name_to_dim=None,
+ type_for_mle="GEV") -> None:
+ super().__init__(result_from_fit, gev_param_name_to_dim)
+ self.type_for_mle = type_for_mle
+
@property
def margin_coef_ordered_dict(self):
values = self.name_to_value['results']
d = self.get_python_dictionary(values)
values = {i: param for i, param in enumerate(np.array(d['par']))}
- return get_margin_coef_ordered_dict(self.gev_param_name_to_dim, values)
+ return get_margin_coef_ordered_dict(self.gev_param_name_to_dim, values, self.type_for_mle)
def _confidence_interval_method(self, common_kwargs, ci_method, return_period):
method_name = ci_method_to_method_name[ci_method]
diff --git a/extreme_fit/model/result_from_model_fit/utils.py b/extreme_fit/model/result_from_model_fit/utils.py
index 7dabcc92..03d95a07 100644
--- a/extreme_fit/model/result_from_model_fit/utils.py
+++ b/extreme_fit/model/result_from_model_fit/utils.py
@@ -11,15 +11,19 @@ def convertFloatVector_to_float(f):
return np.array(f)[0]
-def get_margin_coef_ordered_dict(gev_param_name_to_dim, mle_values):
+def get_margin_coef_ordered_dict(gev_param_name_to_dim, mle_values, type_for_mle="GEV"):
assert gev_param_name_to_dim is not None
# Build the Coeff dict from gev_param_name_to_dim
coef_dict = OrderedDict()
i = 0
for gev_param_name in GevParams.PARAM_NAMES:
- # Add intercept
+ # Add intercept (i.e. stationary parameter)
intercept_coef_name = LinearCoef.coef_template_str(gev_param_name, LinearCoef.INTERCEPT_NAME).format(1)
- coef_dict[intercept_coef_name] = mle_values[i]
+ if type_for_mle == "Gumbel" and i == 2:
+ coef_value = 0
+ else:
+ coef_value = mle_values[i]
+ coef_dict[intercept_coef_name] = coef_value
i += 1
# Add a potential linear temporal trend
if gev_param_name in gev_param_name_to_dim:
diff --git a/test/test_extreme_fit/test_distribution/test_gev/test_gev_temporal_extremes_gumbel.py b/test/test_extreme_fit/test_distribution/test_gev/test_gev_temporal_extremes_gumbel.py
new file mode 100644
index 00000000..2996bb21
--- /dev/null
+++ b/test/test_extreme_fit/test_distribution/test_gev/test_gev_temporal_extremes_gumbel.py
@@ -0,0 +1,51 @@
+import unittest
+
+import numpy as np
+import pandas as pd
+
+from experiment.trend_analysis.univariate_test.abstract_gev_trend_test import fitted_linear_margin_estimator
+from extreme_fit.model.margin_model.linear_margin_model.abstract_temporal_linear_margin_model import \
+ TemporalMarginFitMethod
+from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import StationaryTemporalModel, \
+ NonStationaryLocationTemporalModel, NonStationaryLocationAndScaleTemporalModel, GumbelTemporalModel
+from extreme_fit.model.utils import r, set_seed_r
+from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
+from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_coordinates import \
+ AbstractTemporalCoordinates
+from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
+from spatio_temporal_dataset.spatio_temporal_observations.abstract_spatio_temporal_observations import \
+ AbstractSpatioTemporalObservations
+
+
+class TestGevTemporalExtremesGumbel(unittest.TestCase):
+
+ def setUp(self) -> None:
+ set_seed_r()
+ r("""
+ N <- 50
+ loc = 0; scale = 1; shape <- 0
+ x_gev <- rgev(N, loc = loc, scale = scale, shape = shape)
+ start_loc = 0; start_scale = 1; start_shape = 1
+ """)
+ # Compute the stationary temporal margin with isMev
+ self.start_year = 0
+ df = pd.DataFrame({AbstractCoordinates.COORDINATE_T: range(self.start_year, self.start_year + 50)})
+ self.coordinates = AbstractTemporalCoordinates.from_df(df)
+ df2 = pd.DataFrame(data=np.array(r['x_gev']), index=df.index)
+ observations = AbstractSpatioTemporalObservations(df_maxima_gev=df2)
+ self.dataset = AbstractDataset(observations=observations, coordinates=self.coordinates)
+
+ def test_gev_temporal_margin_fit(self):
+ # Create estimator
+ estimator = fitted_linear_margin_estimator(GumbelTemporalModel,
+ self.coordinates, self.dataset,
+ starting_year=0,
+ fit_method=TemporalMarginFitMethod.extremes_fevd_mle)
+ ref = {'loc': -0.0862185692806497, 'scale': 1.0818465357627252, 'shape': 0}
+ for year in range(1, 3):
+ mle_params_estimated = estimator.margin_function_from_fit.get_gev_params(np.array([year])).to_dict()
+ for key in ref.keys():
+ self.assertAlmostEqual(ref[key], mle_params_estimated[key], places=3)
+
+if __name__ == '__main__':
+ unittest.main()
--
GitLab