From f23cb091c354ebde0aa5c6fb79b057d20b54f28c Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Tue, 17 Dec 2019 12:06:53 +0100
Subject: [PATCH] [PAPER 1] add comparison with gumbel model.
---
.../scm_models_data/abstract_study.py | 3 +-
...dy_visualizer_for_non_stationary_trends.py | 25 +++++++-----
.../abstract_gev_trend_test.py | 2 +-
.../gev_trend_test_one_parameter.py | 8 ++--
.../gumbel_trend_test_one_parameter.py | 24 +++++++++++
....py => gev_trend_test_three_parameters.py} | 6 ++-
.../gev_trend_test_two_parameters.py | 40 ++++++++++++++++++-
.../gumbel_test_two_parameters.py | 8 ++++
extreme_fit/distribution/gev/gev_params.py | 7 +++-
9 files changed, 103 insertions(+), 20 deletions(-)
rename experiment/trend_analysis/univariate_test/extreme_trend_test/trend_test_three_parameters/{gumbel_trend_test_three_parameters.py => gev_trend_test_three_parameters.py} (94%)
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 d511884e..cc7fe98e 100644
--- a/experiment/meteo_france_data/scm_models_data/abstract_study.py
+++ b/experiment/meteo_france_data/scm_models_data/abstract_study.py
@@ -404,7 +404,8 @@ class AbstractStudy(object):
# Add legend for the marker
if massif_name_to_marker_style is not None:
legend_elements = cls.get_legend_for_model_symbol(marker_style_to_label_name, markersize=8)
- ax.legend(handles=legend_elements, bbox_to_anchor=(0.01, 0.03), loc='lower left')
+ ax.legend(handles=legend_elements[:], loc='upper right')
+ # ax.legend(handles=legend_elements[4:], bbox_to_anchor=(0.01, 0.03), loc='lower left')
ax.annotate("Filled symbol = significant trend ", xy=(0.05, 0.015), xycoords='axes fraction', fontsize=7)
if show:
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 ba766f4c..acb422a3 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
@@ -1,3 +1,5 @@
+from collections import OrderedDict
+
import matplotlib.pyplot as plt
from multiprocessing.pool import Pool
from typing import Dict
@@ -19,10 +21,10 @@ from experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_one
GevLocationTrendTest, GevScaleTrendTest
from experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_one_parameter.gumbel_trend_test_one_parameter import \
GumbelLocationTrendTest, GevStationaryVersusGumbel, GumbelScaleTrendTest, GumbelVersusGumbel
-from experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_three_parameters.gumbel_trend_test_three_parameters import \
+from experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_three_parameters.gev_trend_test_three_parameters import \
GevLocationAndScaleTrendTestAgainstGumbel
from experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_two_parameters.gev_trend_test_two_parameters import \
- GevLocationAndScaleTrendTest
+ GevLocationAndScaleTrendTest, GevLocationAgainstGumbel, GevScaleAgainstGumbel
from experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_two_parameters.gumbel_test_two_parameters import \
GumbelLocationAndScaleTrendTest
from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import StationaryTemporalModel
@@ -67,17 +69,16 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
self.uncertainty_massif_names = self.study.study_massif_names
if self.non_stationary_trend_test_to_marker is None:
# Assign default argument for the non stationary trends
- # self.non_stationary_trend_test = [GumbelVersusGumbel,
- # GumbelLocationTrendTest, GumbelScaleTrendTest, GumbelLocationAndScaleTrendTest,
- # GevStationaryVersusGumbel,
- # GevLocationTrendTest, GevScaleTrendTest, GevLocationAndScaleTrendTest,
- # ]
- self.non_stationary_trend_test = [GumbelVersusGumbel, GevLocationAndScaleTrendTestAgainstGumbel]
+ self.non_stationary_trend_test = [GumbelVersusGumbel,
+ GumbelLocationTrendTest, GumbelScaleTrendTest, GumbelLocationAndScaleTrendTest,
+ GevStationaryVersusGumbel,
+ GevLocationAgainstGumbel, GevScaleAgainstGumbel, GevLocationAndScaleTrendTestAgainstGumbel
+ ]
self.non_stationary_trend_test_to_marker = {t: t.marker for t in self.non_stationary_trend_test}
# ["v", "^", "D", "X", "x", 7, 6, "d"]))
else:
self.non_stationary_trend_test = list(self.non_stationary_trend_test_to_marker.keys())
- self.marker_to_label = {t.marker: t.label for t in self.non_stationary_trend_test}
+ self.marker_to_label = OrderedDict([(t.marker, t.label) for t in self.non_stationary_trend_test])
self.global_max_abs_change = None
# Utils
@@ -136,12 +137,14 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
def plot_trends(self, max_abs_tdrl=None, add_colorbar=True):
if max_abs_tdrl is not None:
self.global_max_abs_change = max_abs_tdrl
+ marker_style_selected = set([d['marker'] for d in self.massif_name_to_marker_style.values()])
+ marker_style_to_label_name = {m: l for m, l in self.marker_to_label.items() if m in marker_style_selected}
ax = self.study.visualize_study(massif_name_to_value=self.massif_name_to_change_value,
replace_blue_by_white=False,
axis_off=False, show_label=False,
add_colorbar=add_colorbar,
massif_name_to_marker_style=self.massif_name_to_marker_style,
- marker_style_to_label_name=self.marker_to_label,
+ marker_style_to_label_name=marker_style_to_label_name,
massif_name_to_color=self.massif_name_to_color,
cmap=self.cmap,
show=False,
@@ -227,7 +230,7 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
for m, t in self.massif_name_to_minimized_aic_non_stationary_trend_test.items():
d[m] = {'marker': self.non_stationary_trend_test_to_marker[type(t)],
'color': 'k',
- 'markersize': 5,
+ 'markersize': 7,
'fillstyle': 'full' if t.is_significant else 'none'}
return d
diff --git a/experiment/trend_analysis/univariate_test/extreme_trend_test/abstract_gev_trend_test.py b/experiment/trend_analysis/univariate_test/extreme_trend_test/abstract_gev_trend_test.py
index 907716d9..5411a6ed 100644
--- a/experiment/trend_analysis/univariate_test/extreme_trend_test/abstract_gev_trend_test.py
+++ b/experiment/trend_analysis/univariate_test/extreme_trend_test/abstract_gev_trend_test.py
@@ -88,7 +88,7 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
@property
def total_number_of_parameters_for_unconstrained_model(self) -> int:
- return self.degree_freedom_chi2 + 3
+ raise NotImplementedError
@property
def aic(self):
diff --git a/experiment/trend_analysis/univariate_test/extreme_trend_test/trend_test_one_parameter/gev_trend_test_one_parameter.py b/experiment/trend_analysis/univariate_test/extreme_trend_test/trend_test_one_parameter/gev_trend_test_one_parameter.py
index d09c5a3b..b34c2bc7 100644
--- a/experiment/trend_analysis/univariate_test/extreme_trend_test/trend_test_one_parameter/gev_trend_test_one_parameter.py
+++ b/experiment/trend_analysis/univariate_test/extreme_trend_test/trend_test_one_parameter/gev_trend_test_one_parameter.py
@@ -4,6 +4,7 @@ from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_m
NonStationaryLocationTemporalModel, NonStationaryScaleTemporalModel, NonStationaryShapeTemporalModel, \
StationaryTemporalModel
from extreme_fit.distribution.gev.gev_params import GevParams
+from root_utils import classproperty
class GevTrendTestOneParameter(AbstractGevTrendTest):
@@ -31,9 +32,10 @@ class GevTrendTestOneParameterAgainstStationary(GevTrendTestOneParameter):
class GevLocationTrendTest(GevTrendTestOneParameterAgainstStationary):
- def __init__(self, years, maxima, starting_year, quantile_level=EUROCODE_QUANTILE):
+ def __init__(self, years, maxima, starting_year, quantile_level=EUROCODE_QUANTILE, constrained_model_class=StationaryTemporalModel):
super().__init__(years, maxima, starting_year,
unconstrained_model_class=NonStationaryLocationTemporalModel,
+ constrained_model_class=constrained_model_class,
gev_param_name=GevParams.LOC,
quantile_level=quantile_level)
@@ -51,12 +53,12 @@ class GevLocationTrendTest(GevTrendTestOneParameterAgainstStationary):
def variance_difference_same_sign_as_slope_strenght(self) -> bool:
return False
-
class GevScaleTrendTest(GevTrendTestOneParameterAgainstStationary):
- def __init__(self, years, maxima, starting_year, quantile_level=EUROCODE_QUANTILE):
+ def __init__(self, years, maxima, starting_year, quantile_level=EUROCODE_QUANTILE, constrained_model_class=StationaryTemporalModel):
super().__init__(years, maxima, starting_year,
unconstrained_model_class=NonStationaryScaleTemporalModel,
+ constrained_model_class=constrained_model_class,
gev_param_name=GevParams.SCALE,
quantile_level=quantile_level)
diff --git a/experiment/trend_analysis/univariate_test/extreme_trend_test/trend_test_one_parameter/gumbel_trend_test_one_parameter.py b/experiment/trend_analysis/univariate_test/extreme_trend_test/trend_test_one_parameter/gumbel_trend_test_one_parameter.py
index 71df5043..698a7063 100644
--- a/experiment/trend_analysis/univariate_test/extreme_trend_test/trend_test_one_parameter/gumbel_trend_test_one_parameter.py
+++ b/experiment/trend_analysis/univariate_test/extreme_trend_test/trend_test_one_parameter/gumbel_trend_test_one_parameter.py
@@ -55,6 +55,14 @@ class GevStationaryVersusGumbel(GevTrendTestOneParameter):
def _slope_strength(self):
return 0.0
+ @classproperty
+ def label(self):
+ return super().label % '\\zeta_0'
+
+ @classproperty
+ def marker(self):
+ return 'X'
+
class GumbelLocationTrendTest(GevTrendTestOneParameterAgainstStationary):
def __init__(self, years, maxima, starting_year, quantile_level=EUROCODE_QUANTILE):
@@ -72,6 +80,14 @@ class GumbelLocationTrendTest(GevTrendTestOneParameterAgainstStationary):
return self.unconstrained_estimator_gev_params.time_derivative_of_return_level(p=self.quantile_level,
mu1=self.non_stationary_linear_coef)
+ @classproperty
+ def label(self):
+ return super().label % '\\mu_1'
+
+ @classproperty
+ def marker(self):
+ return '.'
+
class GumbelScaleTrendTest(GevTrendTestOneParameterAgainstStationary):
@@ -90,3 +106,11 @@ class GumbelScaleTrendTest(GevTrendTestOneParameterAgainstStationary):
@property
def total_number_of_parameters_for_unconstrained_model(self) -> int:
return 3
+
+ @classproperty
+ def label(self):
+ return super().label % '\\sigma_1'
+
+ @classproperty
+ def marker(self):
+ return 11
\ No newline at end of file
diff --git a/experiment/trend_analysis/univariate_test/extreme_trend_test/trend_test_three_parameters/gumbel_trend_test_three_parameters.py b/experiment/trend_analysis/univariate_test/extreme_trend_test/trend_test_three_parameters/gev_trend_test_three_parameters.py
similarity index 94%
rename from experiment/trend_analysis/univariate_test/extreme_trend_test/trend_test_three_parameters/gumbel_trend_test_three_parameters.py
rename to experiment/trend_analysis/univariate_test/extreme_trend_test/trend_test_three_parameters/gev_trend_test_three_parameters.py
index b5f6cbcc..cb951aad 100644
--- a/experiment/trend_analysis/univariate_test/extreme_trend_test/trend_test_three_parameters/gumbel_trend_test_three_parameters.py
+++ b/experiment/trend_analysis/univariate_test/extreme_trend_test/trend_test_three_parameters/gev_trend_test_three_parameters.py
@@ -45,4 +45,8 @@ class GevLocationAndScaleTrendTestAgainstGumbel(GevTrendTestThreeParameters):
@classproperty
def marker(self):
- return 'd'
+ return 'D'
+
+ @property
+ def total_number_of_parameters_for_unconstrained_model(self) -> int:
+ return 5
\ No newline at end of file
diff --git a/experiment/trend_analysis/univariate_test/extreme_trend_test/trend_test_two_parameters/gev_trend_test_two_parameters.py b/experiment/trend_analysis/univariate_test/extreme_trend_test/trend_test_two_parameters/gev_trend_test_two_parameters.py
index 89cd18b3..7ef18928 100644
--- a/experiment/trend_analysis/univariate_test/extreme_trend_test/trend_test_two_parameters/gev_trend_test_two_parameters.py
+++ b/experiment/trend_analysis/univariate_test/extreme_trend_test/trend_test_two_parameters/gev_trend_test_two_parameters.py
@@ -1,9 +1,12 @@
from experiment.eurocode_data.utils import EUROCODE_QUANTILE
from experiment.trend_analysis.univariate_test.extreme_trend_test.abstract_gev_trend_test import AbstractGevTrendTest
+from experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_one_parameter.gev_trend_test_one_parameter import \
+ GevLocationTrendTest, GevScaleTrendTest
from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import \
NonStationaryLocationAndScaleTemporalModel, StationaryTemporalModel, NonStationaryLocationAndScaleGumbelModel, \
GumbelTemporalModel
from extreme_fit.distribution.gev.gev_params import GevParams
+from root_utils import classproperty
class GevTrendTestTwoParameters(AbstractGevTrendTest):
@@ -15,7 +18,8 @@ class GevTrendTestTwoParameters(AbstractGevTrendTest):
class GevLocationAndScaleTrendTest(GevTrendTestTwoParameters):
- def __init__(self, years, maxima, starting_year, constrained_model_class=StationaryTemporalModel, quantile_level=EUROCODE_QUANTILE):
+ def __init__(self, years, maxima, starting_year, constrained_model_class=StationaryTemporalModel,
+ quantile_level=EUROCODE_QUANTILE):
super().__init__(years, maxima, starting_year,
unconstrained_model_class=NonStationaryLocationAndScaleTemporalModel,
constrained_model_class=constrained_model_class,
@@ -45,3 +49,37 @@ class GevLocationAndScaleTrendTest(GevTrendTestTwoParameters):
return self.same_sign(self.sigma1, self._slope_strength())
+class GevLocationAgainstGumbel(GevTrendTestTwoParameters, GevLocationTrendTest):
+
+ def __init__(self, years, maxima, starting_year, quantile_level=EUROCODE_QUANTILE):
+ super().__init__(years, maxima, starting_year, quantile_level, GumbelTemporalModel)
+
+ @classproperty
+ def label(self):
+ return super().label % '\\zeta_0, \\mu_1'
+
+ @classproperty
+ def marker(self):
+ return 'o'
+
+ @property
+ def total_number_of_parameters_for_unconstrained_model(self) -> int:
+ return 4
+
+
+class GevScaleAgainstGumbel(GevTrendTestTwoParameters, GevScaleTrendTest):
+
+ def __init__(self, years, maxima, starting_year, quantile_level=EUROCODE_QUANTILE):
+ super().__init__(years, maxima, starting_year, quantile_level, GumbelTemporalModel)
+
+ @classproperty
+ def label(self):
+ return super().label % '\\zeta_0, \\sigma_1'
+
+ @classproperty
+ def marker(self):
+ return '^'
+
+ @property
+ def total_number_of_parameters_for_unconstrained_model(self) -> int:
+ return 4
diff --git a/experiment/trend_analysis/univariate_test/extreme_trend_test/trend_test_two_parameters/gumbel_test_two_parameters.py b/experiment/trend_analysis/univariate_test/extreme_trend_test/trend_test_two_parameters/gumbel_test_two_parameters.py
index ea7648ac..8eced3b2 100644
--- a/experiment/trend_analysis/univariate_test/extreme_trend_test/trend_test_two_parameters/gumbel_test_two_parameters.py
+++ b/experiment/trend_analysis/univariate_test/extreme_trend_test/trend_test_two_parameters/gumbel_test_two_parameters.py
@@ -4,6 +4,7 @@ from experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_two
from extreme_fit.distribution.gev.gev_params import GevParams
from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import \
NonStationaryLocationAndScaleGumbelModel, GumbelTemporalModel
+from root_utils import classproperty
class GumbelLocationAndScaleTrendTest(GevTrendTestTwoParameters):
@@ -30,3 +31,10 @@ class GumbelLocationAndScaleTrendTest(GevTrendTestTwoParameters):
return self.unconstrained_estimator_gev_params.time_derivative_of_return_level(p=self.quantile_level,
mu1=self.mu1,
sigma1=self.sigma1)
+ @classproperty
+ def label(self):
+ return super().label % '\\mu_1, \\sigma_1'
+
+ @classproperty
+ def marker(self):
+ return 'd'
diff --git a/extreme_fit/distribution/gev/gev_params.py b/extreme_fit/distribution/gev/gev_params.py
index 55f56115..5e957005 100644
--- a/extreme_fit/distribution/gev/gev_params.py
+++ b/extreme_fit/distribution/gev/gev_params.py
@@ -61,8 +61,11 @@ class GevParams(AbstractParams):
"""
quantile_annual_variation = mu1
if sigma1 != 0:
- power = np.float_power(- np.log(p), -self.shape)
- quantile_annual_variation -= (sigma1 / self.shape) * (1 - power)
+ if self.shape == 0:
+ quantile_annual_variation -= sigma1 * np.log(- np.log(p))
+ else:
+ power = np.float_power(- np.log(p), -self.shape)
+ quantile_annual_variation -= (sigma1 / self.shape) * (1 - power)
return quantile_annual_variation
# Compute some indicators (such as the mean and the variance)
--
GitLab