From f10e815593cc9538db5607d0d141310f6908f963 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Mon, 25 Nov 2019 13:34:55 +0100
Subject: [PATCH] [PAPER 1] taking into account psnow when computing TDRL and
also when computing return level uncertainties
---
.../altitude_hypercube_visualizer.py | 2 +-
.../eurocode_visualizer.py | 45 +++++++++++--------
.../main_result_trends_and_return_levels.py | 14 +++---
...dy_visualizer_for_non_stationary_trends.py | 10 +++--
...bstract_comparison_non_stationary_model.py | 11 +++--
.../abstract_gev_trend_test.py | 12 ++---
.../gev_trend_test_one_parameter.py | 33 +++++++++-----
.../gev_trend_test_two_parameters.py | 13 +++---
extreme_fit/distribution/gev/gev_params.py | 3 +-
.../abstract_extract_eurocode_return_level.py | 12 ++---
.../abstract_result_from_extremes.py | 4 +-
.../eurocode_return_level_uncertainties.py | 17 ++++---
.../result_from_bayesian_extremes.py | 4 +-
.../result_from_mle_extremes.py | 4 +-
14 files changed, 107 insertions(+), 77 deletions(-)
diff --git a/experiment/meteo_france_data/scm_models_data/visualization/hypercube_visualization/altitude_hypercube_visualizer.py b/experiment/meteo_france_data/scm_models_data/visualization/hypercube_visualization/altitude_hypercube_visualizer.py
index d918abe9..ae0e5389 100644
--- a/experiment/meteo_france_data/scm_models_data/visualization/hypercube_visualization/altitude_hypercube_visualizer.py
+++ b/experiment/meteo_france_data/scm_models_data/visualization/hypercube_visualization/altitude_hypercube_visualizer.py
@@ -396,7 +396,7 @@ class AltitudeHypercubeVisualizer(AbstractHypercubeVisualizer):
title += ' until starting_year={}'.format(self.last_starting_year)
title += ' with {} test'.format(get_display_name_from_object_type(self.trend_test_class))
if first_title:
- title += '\nEvolution of the quantile {} every {} years'.format(AbstractGevTrendTest.quantile_for_strength,
+ title += '\nEvolution of the Eurocode quantile every {} years'.format(
AbstractGevTrendTest.nb_years_for_quantile_evolution)
else:
title += '\nStarting years'
diff --git a/experiment/paper_past_snow_loads/result_trends_and_return_levels/eurocode_visualizer.py b/experiment/paper_past_snow_loads/result_trends_and_return_levels/eurocode_visualizer.py
index 78bd644d..fcf66c6f 100644
--- a/experiment/paper_past_snow_loads/result_trends_and_return_levels/eurocode_visualizer.py
+++ b/experiment/paper_past_snow_loads/result_trends_and_return_levels/eurocode_visualizer.py
@@ -108,25 +108,10 @@ def plot_single_uncertainty_massif_and_non_stationary_context(ax, massif_name, n
if j == 0:
eurocode_region.plot_eurocode_snow_load_on_ground_characteristic_value_variable_action(ax,
altitudes=altitudes)
- # Compute ordered_return_level_uncertaines for a given massif_name, uncertainty methods, and non stationary context
- ordered_return_level_uncertaines = []
- for visualizer in altitude_to_visualizer.values():
- u = visualizer.triplet_to_eurocode_uncertainty[(uncertainty_method, non_stationary_context, massif_name)]
- ordered_return_level_uncertaines.append(u)
- # Display
- mean = [r.mean_estimate for r in ordered_return_level_uncertaines]
- # Filter and keep only non nan values
- not_nan_index = [not np.isnan(m) for m in mean]
- mean = list(np.array(mean)[not_nan_index])
- altitudes = list(np.array(altitudes)[not_nan_index])
- ordered_return_level_uncertaines = list(np.array(ordered_return_level_uncertaines)[not_nan_index])
-
- ci_method_name = str(uncertainty_method).split('.')[1].replace('_', ' ')
- ax.plot(altitudes, mean, linestyle='--', marker='o', color=color,
- label=get_label_name(non_stationary_context, ci_method_name))
- lower_bound = [r.confidence_interval[0] for r in ordered_return_level_uncertaines]
- upper_bound = [r.confidence_interval[1] for r in ordered_return_level_uncertaines]
- ax.fill_between(altitudes, lower_bound, upper_bound, color=color, alpha=alpha)
+ # Plot uncertainties
+ plot_valid_return_level_uncertainties(alpha, altitude_to_visualizer, altitudes, ax, color, massif_name,
+ non_stationary_context, uncertainty_method)
+
ax.legend(loc=2)
ax.set_ylim([0.0, 16])
massif_name_str = massif_name.replace('_', ' ')
@@ -144,3 +129,25 @@ def plot_single_uncertainty_massif_and_non_stationary_context(ax, massif_name, n
ax.set_ylabel('50-year return level of SL (kN $m^-2$)')
ax.set_xlabel('Altitude (m)')
ax.grid()
+
+
+def plot_valid_return_level_uncertainties(alpha, altitude_to_visualizer, altitudes, ax, color, massif_name,
+ non_stationary_context, uncertainty_method):
+ # Compute ordered_return_level_uncertaines for a given massif_name, uncertainty methods, and non stationary context
+ ordered_return_level_uncertainties = []
+ for visualizer in altitude_to_visualizer.values():
+ u = visualizer.triplet_to_eurocode_uncertainty[(uncertainty_method, non_stationary_context, massif_name)]
+ ordered_return_level_uncertainties.append(u)
+ # Display
+ mean = [r.mean_estimate for r in ordered_return_level_uncertainties]
+ # Filter and keep only non nan values
+ not_nan_index = [not np.isnan(m) for m in mean]
+ mean = list(np.array(mean)[not_nan_index])
+ valid_altitudes = list(np.array(altitudes)[not_nan_index])
+ ordered_return_level_uncertainties = list(np.array(ordered_return_level_uncertainties)[not_nan_index])
+ ci_method_name = str(uncertainty_method).split('.')[1].replace('_', ' ')
+ ax.plot(valid_altitudes, mean, linestyle='--', marker='o', color=color,
+ label=get_label_name(non_stationary_context, ci_method_name))
+ lower_bound = [r.confidence_interval[0] for r in ordered_return_level_uncertainties]
+ upper_bound = [r.confidence_interval[1] for r in ordered_return_level_uncertainties]
+ ax.fill_between(valid_altitudes, lower_bound, upper_bound, color=color, alpha=alpha)
diff --git a/experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py b/experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py
index 8554af5c..6c867dfb 100644
--- a/experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py
+++ b/experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py
@@ -67,16 +67,16 @@ def major_result():
if __name__ == '__main__':
# minor_result(altitude=1800)
- # intermediate_result(altitudes=[1500, 1800], massif_names=['Chartreuse'],
- # uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.ci_mle],
- # non_stationary_uncertainty=[False])
+ intermediate_result(altitudes=[1500, 1800], massif_names=['Chartreuse'],
+ uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.ci_mle],
+ non_stationary_uncertainty=[False])
# intermediate_result(altitudes=[1500, 1800], massif_names=None,
# uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.ci_mle],
# non_stationary_uncertainty=[False])
# intermediate_result(altitudes=[300, 600, 900, 1200, 1500, 1800], massif_names=None,
# uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.ci_mle],
# non_stationary_uncertainty=[False])
- intermediate_result(altitudes=[300, 600, 900, 1200, 1500, 1800], massif_names=None,
- uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.ci_mle,
- ConfidenceIntervalMethodFromExtremes.ci_bayes],
- non_stationary_uncertainty=[False, True])
+ # intermediate_result(altitudes=[300, 600, 900, 1200, 1500, 1800], massif_names=['Vercors', 'Oisans', 'Devoluy'],
+ # uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.ci_mle,
+ # ConfidenceIntervalMethodFromExtremes.ci_bayes],
+ # non_stationary_uncertainty=[False, True])
diff --git a/experiment/paper_past_snow_loads/result_trends_and_return_levels/study_visualizer_for_non_stationary_trends.py b/experiment/paper_past_snow_loads/result_trends_and_return_levels/study_visualizer_for_non_stationary_trends.py
index b80f46ac..85272623 100644
--- a/experiment/paper_past_snow_loads/result_trends_and_return_levels/study_visualizer_for_non_stationary_trends.py
+++ b/experiment/paper_past_snow_loads/result_trends_and_return_levels/study_visualizer_for_non_stationary_trends.py
@@ -72,7 +72,7 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
self.massif_name_to_years_and_maxima.items()}
@cached_property
- def massif_name_to_eurocode_quantile_in_practice(self):
+ def massif_name_to_eurocode_quantile_level_in_practice(self):
"""Due to missing data, the the eurocode quantile which 0.98 if we have all the data
correspond in practice to the quantile psnow x 0.98 of the data where there is snow"""
return {m: p * EUROCODE_QUANTILE for m, p in self.massif_name_to_psnow.items()}
@@ -90,7 +90,9 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
starting_year = None
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]
+ quantile_level = self.massif_name_to_eurocode_quantile_level_in_practice[massif_name]
+ non_stationary_trend_test = [t(years=x, maxima=y, starting_year=starting_year, quantile_level=quantile_level)
+ 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
@@ -151,7 +153,9 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
arguments = [
[self.massif_name_to_non_null_years_and_maxima[m],
self.massif_name_to_model_class(m, non_stationary_model),
- ci_method, self.effective_temporal_covariate] for m in self.uncertainty_massif_names]
+ ci_method, self.effective_temporal_covariate,
+ self.massif_name_to_eurocode_quantile_level_in_practice[m]
+ ] for m in self.uncertainty_massif_names]
if self.multiprocessing:
with Pool(NB_CORES) as p:
res = p.starmap(compute_eurocode_confidence_interval, arguments)
diff --git a/experiment/trend_analysis/univariate_test/abstract_comparison_non_stationary_model.py b/experiment/trend_analysis/univariate_test/abstract_comparison_non_stationary_model.py
index 51f0cc77..cfd2dc5e 100644
--- a/experiment/trend_analysis/univariate_test/abstract_comparison_non_stationary_model.py
+++ b/experiment/trend_analysis/univariate_test/abstract_comparison_non_stationary_model.py
@@ -1,3 +1,4 @@
+from experiment.eurocode_data.utils import EUROCODE_QUANTILE
from experiment.trend_analysis.univariate_test.abstract_gev_trend_test import AbstractGevTrendTest
from experiment.trend_analysis.univariate_test.gev_trend_test_two_parameters import GevLocationAndScaleTrendTest
from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import \
@@ -27,11 +28,13 @@ class AbstractComparisonNonStationaryModelOneParameter(AbstractComparisonNonStat
class ComparisonAgainstMu(AbstractComparisonNonStationaryModelOneParameter, GevLocationAndScaleTrendTest):
- def __init__(self, years, maxima, starting_year):
- super().__init__(years, maxima, starting_year, constrained_model_class=NonStationaryLocationTemporalModel)
+ def __init__(self, years, maxima, starting_year, quantile_level=EUROCODE_QUANTILE):
+ super().__init__(years, maxima, starting_year, constrained_model_class=NonStationaryLocationTemporalModel,
+ quantile_level=quantile_level)
class ComparisonAgainstSigma(AbstractComparisonNonStationaryModelOneParameter, GevLocationAndScaleTrendTest):
- def __init__(self, years, maxima, starting_year):
- super().__init__(years, maxima, starting_year, constrained_model_class=NonStationaryScaleTemporalModel)
+ def __init__(self, years, maxima, starting_year, quantile_level=EUROCODE_QUANTILE):
+ super().__init__(years, maxima, starting_year, constrained_model_class=NonStationaryScaleTemporalModel,
+ quantile_level=quantile_level)
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 bfe20929..59184411 100644
--- a/experiment/trend_analysis/univariate_test/abstract_gev_trend_test.py
+++ b/experiment/trend_analysis/univariate_test/abstract_gev_trend_test.py
@@ -19,14 +19,14 @@ from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoo
class AbstractGevTrendTest(AbstractUnivariateTest):
RRunTimeError_TREND = 'R RunTimeError trend'
- # I should use the quantile from the Eurocode for the buildings
- quantile_for_strength = EUROCODE_QUANTILE
nb_years_for_quantile_evolution = 10
def __init__(self, years, maxima, starting_year, unconstrained_model_class,
constrained_model_class=StationaryTemporalModel,
+ quantile_level=EUROCODE_QUANTILE,
):
super().__init__(years, maxima, starting_year)
+ self.quantile_level = quantile_level
self.unconstrained_model_class = unconstrained_model_class
self.constrained_model_class = constrained_model_class
self.fit_method = TemporalMarginFitMethod.extremes_fevd_mle
@@ -42,14 +42,16 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
@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)
+ 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)
+ return fitted_linear_margin_estimator(self.unconstrained_model_class, self.coordinates, self.dataset,
+ self.starting_year, self.fit_method)
except SafeRunException:
self.crashed = True
@@ -176,4 +178,4 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
if self.crashed:
return 0.0
else:
- return self.non_stationary_constant_gev_params.quantile(p=self.quantile_for_strength)
+ return self.non_stationary_constant_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 9883777a..8af03890 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,3 +1,4 @@
+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
@@ -6,8 +7,10 @@ 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):
- super().__init__(years, maxima, starting_year, unconstrained_model_class)
+ 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
@@ -21,13 +24,15 @@ class GevTrendTestOneParameter(AbstractGevTrendTest):
class GevLocationTrendTest(GevTrendTestOneParameter):
- def __init__(self, years, maxima, starting_year):
+ def __init__(self, years, maxima, starting_year, quantile_level=EUROCODE_QUANTILE):
super().__init__(years, maxima, starting_year,
- NonStationaryLocationTemporalModel, GevParams.LOC)
+ unconstrained_model_class=NonStationaryLocationTemporalModel,
+ gev_param_name=GevParams.LOC,
+ quantile_level=quantile_level)
def _slope_strength(self):
- return self.non_stationary_constant_gev_params.quantile_strength_evolution(p=self.quantile_for_strength,
- mu1=self.non_stationary_linear_coef)
+ return self.non_stationary_constant_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:
@@ -42,13 +47,15 @@ class GevLocationTrendTest(GevTrendTestOneParameter):
class GevScaleTrendTest(GevTrendTestOneParameter):
- def __init__(self, years, maxima, starting_year):
+ def __init__(self, years, maxima, starting_year, quantile_level=EUROCODE_QUANTILE):
super().__init__(years, maxima, starting_year,
- NonStationaryScaleTemporalModel, GevParams.SCALE)
+ unconstrained_model_class=NonStationaryScaleTemporalModel,
+ gev_param_name=GevParams.SCALE,
+ quantile_level=quantile_level)
def _slope_strength(self):
- return self.non_stationary_constant_gev_params.quantile_strength_evolution(
- p=self.quantile_for_strength,
+ return self.non_stationary_constant_gev_params.time_derivative_of_return_level(
+ p=self.quantile_level,
sigma1=self.non_stationary_linear_coef)
@property
@@ -65,6 +72,8 @@ class GevScaleTrendTest(GevTrendTestOneParameter):
class GevShapeTrendTest(GevTrendTestOneParameter):
- def __init__(self, years, maxima, starting_year):
+ def __init__(self, years, maxima, starting_year, quantile_level=EUROCODE_QUANTILE):
super().__init__(years, maxima, starting_year,
- NonStationaryShapeTemporalModel, GevParams.SHAPE)
+ unconstrained_model_class=NonStationaryShapeTemporalModel,
+ gev_param_name=GevParams.SHAPE,
+ 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 85456427..95525524 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
@@ -1,3 +1,4 @@
+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 \
NonStationaryLocationAndScaleTemporalModel, StationaryTemporalModel
@@ -13,9 +14,11 @@ class GevTrendTestTwoParameters(AbstractGevTrendTest):
class GevLocationAndScaleTrendTest(GevTrendTestTwoParameters):
- def __init__(self, years, maxima, starting_year, constrained_model_class=StationaryTemporalModel):
+ def __init__(self, years, maxima, starting_year, constrained_model_class=StationaryTemporalModel, quantile_level=EUROCODE_QUANTILE):
super().__init__(years, maxima, starting_year,
- NonStationaryLocationAndScaleTemporalModel, constrained_model_class=constrained_model_class)
+ unconstrained_model_class=NonStationaryLocationAndScaleTemporalModel,
+ constrained_model_class=constrained_model_class,
+ quantile_level=quantile_level)
@property
def mu1(self):
@@ -26,9 +29,9 @@ 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.quantile_strength_evolution(p=self.quantile_for_strength,
- mu1=self.mu1,
- sigma1=self.sigma1)
+ return self.non_stationary_constant_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:
diff --git a/extreme_fit/distribution/gev/gev_params.py b/extreme_fit/distribution/gev/gev_params.py
index 9811f7dd..55f56115 100644
--- a/extreme_fit/distribution/gev/gev_params.py
+++ b/extreme_fit/distribution/gev/gev_params.py
@@ -49,8 +49,7 @@ class GevParams(AbstractParams):
def __str__(self):
return self.to_dict().__str__()
- def quantile_strength_evolution(self, p=0.99, mu1=0.0, sigma1=0.0):
- # todo: chagne the name to time derivative
+ def time_derivative_of_return_level(self, p=0.99, mu1=0.0, sigma1=0.0):
"""
Compute the variation of some quantile with respect to time.
(when mu1 and sigma1 can be modified with time)
diff --git a/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_extract_eurocode_return_level.py b/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_extract_eurocode_return_level.py
index c0cd76d2..d02ef69d 100644
--- a/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_extract_eurocode_return_level.py
+++ b/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_extract_eurocode_return_level.py
@@ -16,13 +16,13 @@ from extreme_fit.model.result_from_model_fit.result_from_extremes.result_from_ml
class AbstractExtractEurocodeReturnLevel(object):
ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY = 0.05
- def __init__(self, estimator: LinearMarginEstimator, ci_method, temporal_covariate):
+ def __init__(self, estimator: LinearMarginEstimator, ci_method, temporal_covariate, quantile_level=EUROCODE_QUANTILE):
self.ci_method = ci_method
self.estimator = estimator
self.result_from_fit = self.estimator.result_from_model_fit
self.temporal_covariate = temporal_covariate
# Fixed Parameters
- self.eurocode_quantile = EUROCODE_QUANTILE
+ self.eurocode_quantile_level = quantile_level
self.alpha_for_confidence_interval = self.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY
@property
@@ -43,7 +43,7 @@ class ExtractEurocodeReturnLevelFromCiMethod(AbstractExtractEurocodeReturnLevel)
@cached_property
def confidence_interval_method(self):
- return self.result_from_fit.confidence_interval_method(self.eurocode_quantile,
+ return self.result_from_fit.confidence_interval_method(self.eurocode_quantile_level,
self.alpha_for_confidence_interval,
self.transformed_temporal_covariate,
self.ci_method)
@@ -60,8 +60,8 @@ class ExtractEurocodeReturnLevelFromCiMethod(AbstractExtractEurocodeReturnLevel)
class ExtractEurocodeReturnLevelFromMyBayesianExtremes(AbstractExtractEurocodeReturnLevel):
result_from_fit: ResultFromBayesianExtremes
- def __init__(self, estimator: LinearMarginEstimator, ci_method, temporal_covariate):
- super().__init__(estimator, ci_method, temporal_covariate)
+ def __init__(self, estimator: LinearMarginEstimator, ci_method, temporal_covariate, quantile_level=EUROCODE_QUANTILE):
+ super().__init__(estimator, ci_method, temporal_covariate, quantile_level)
assert isinstance(self.result_from_fit, ResultFromBayesianExtremes)
@property
@@ -81,7 +81,7 @@ class ExtractEurocodeReturnLevelFromMyBayesianExtremes(AbstractExtractEurocodeRe
@cached_property
def posterior_eurocode_return_level_samples_for_temporal_covariate(self) -> np.ndarray:
return np.array(
- [p.quantile(self.eurocode_quantile) for p in self.gev_params_from_fit_for_temporal_covariate])
+ [p.quantile(self.eurocode_quantile_level) for p in self.gev_params_from_fit_for_temporal_covariate])
@property
def mean_estimate(self) -> np.ndarray:
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 acf56fa5..240d6535 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
@@ -48,8 +48,8 @@ class AbstractResultFromExtremes(AbstractResultFromModelFit):
qcov = r("make.qcov")(self.result_from_fit,
**kwargs)
common_kwargs['qcov'] = qcov
- mean_estimate, confidence_interval = self._confidence_interval_method(common_kwargs, ci_method)
+ mean_estimate, confidence_interval = self._confidence_interval_method(common_kwargs, ci_method, return_period)
return mean_estimate, confidence_interval
- def _confidence_interval_method(self, common_kwargs, ci_method):
+ def _confidence_interval_method(self, common_kwargs, ci_method, return_period):
raise NotImplementedError
diff --git a/extreme_fit/model/result_from_model_fit/result_from_extremes/eurocode_return_level_uncertainties.py b/extreme_fit/model/result_from_model_fit/result_from_extremes/eurocode_return_level_uncertainties.py
index 71d7e183..7b25ef08 100644
--- a/extreme_fit/model/result_from_model_fit/result_from_extremes/eurocode_return_level_uncertainties.py
+++ b/extreme_fit/model/result_from_model_fit/result_from_extremes/eurocode_return_level_uncertainties.py
@@ -1,5 +1,6 @@
from enum import Enum
+from experiment.eurocode_data.utils import EUROCODE_QUANTILE
from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_extract_eurocode_return_level import \
ExtractEurocodeReturnLevelFromMyBayesianExtremes, ExtractEurocodeReturnLevelFromCiMethod
from experiment.trend_analysis.univariate_test.utils import load_temporal_coordinates_and_dataset, \
@@ -11,9 +12,9 @@ from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_int
ConfidenceIntervalMethodFromExtremes
-def compute_eurocode_confidence_interval(smooth_maxima_x_y, model_class, ci_method, temporal_covariate):
+def compute_eurocode_confidence_interval(smooth_maxima_x_y, model_class, ci_method, temporal_covariate, quantile_level):
years, smooth_maxima = smooth_maxima_x_y
- return EurocodeConfidenceIntervalFromExtremes.from_maxima_years_model_class(smooth_maxima, years, model_class, temporal_covariate, ci_method)
+ return EurocodeConfidenceIntervalFromExtremes.from_maxima_years_model_class(smooth_maxima, years, model_class, temporal_covariate, ci_method, quantile_level)
class EurocodeConfidenceIntervalFromExtremes(object):
@@ -29,17 +30,19 @@ class EurocodeConfidenceIntervalFromExtremes(object):
@classmethod
def from_estimator_extremes(cls, estimator_extremes: LinearMarginEstimator,
ci_method: ConfidenceIntervalMethodFromExtremes,
- temporal_covariate: int):
+ temporal_covariate: int,
+ quantile_level=EUROCODE_QUANTILE):
if ci_method == ConfidenceIntervalMethodFromExtremes.my_bayes:
- extractor = ExtractEurocodeReturnLevelFromMyBayesianExtremes(estimator_extremes, ci_method, temporal_covariate)
+ extractor = ExtractEurocodeReturnLevelFromMyBayesianExtremes(estimator_extremes, ci_method, temporal_covariate, quantile_level)
else:
- extractor = ExtractEurocodeReturnLevelFromCiMethod(estimator_extremes, ci_method, temporal_covariate)
+ extractor = ExtractEurocodeReturnLevelFromCiMethod(estimator_extremes, ci_method, temporal_covariate, quantile_level)
return cls(extractor.mean_estimate, extractor.confidence_interval)
@classmethod
def from_maxima_years_model_class(cls, maxima, years, model_class,
temporal_covariate,
- ci_method=ConfidenceIntervalMethodFromExtremes.ci_bayes):
+ ci_method=ConfidenceIntervalMethodFromExtremes.ci_bayes,
+ quantile_level=EUROCODE_QUANTILE):
# Load coordinates and dataset
coordinates, dataset = load_temporal_coordinates_and_dataset(maxima, years)
# Select fit method depending on the ci_method
@@ -52,7 +55,7 @@ class EurocodeConfidenceIntervalFromExtremes(object):
fitted_estimator = fitted_linear_margin_estimator(model_class, coordinates, dataset, starting_year=None,
fit_method=fit_method)
# Load object from result from extremes
- return cls.from_estimator_extremes(fitted_estimator, ci_method, temporal_covariate)
+ return cls.from_estimator_extremes(fitted_estimator, ci_method, temporal_covariate, quantile_level)
diff --git a/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_bayesian_extremes.py b/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_bayesian_extremes.py
index b83915d8..acb087e6 100644
--- a/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_bayesian_extremes.py
+++ b/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_bayesian_extremes.py
@@ -47,7 +47,7 @@ class ResultFromBayesianExtremes(AbstractResultFromExtremes):
mean_posterior_parameters = self.df_posterior_samples.iloc[:, :-2].mean(axis=0)
return self.get_coef_dict_from_posterior_sample(mean_posterior_parameters)
- def _confidence_interval_method(self, common_kwargs, ci_method):
+ def _confidence_interval_method(self, common_kwargs, ci_method, return_period):
bayesian_ci_parameters = {
'burn.in': self.burn_in_nb,
'FUN': "mean",
@@ -58,7 +58,7 @@ class ResultFromBayesianExtremes(AbstractResultFromExtremes):
lower, mean_estimate, upper = a
else:
d = self.get_python_dictionary(res)
- keys = ['Posterior Mean 50-year level', '95% lower CI', '95% upper CI']
+ keys = ['Posterior Mean {}-year level'.format(return_period), '95% lower CI', '95% upper CI']
mean_estimate, lower, upper = [np.array(d[k])[0] for k in keys]
return mean_estimate, (lower, upper)
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 68708932..243eb2e9 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
@@ -18,7 +18,7 @@ class ResultFromMleExtremes(AbstractResultFromExtremes):
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)
- def _confidence_interval_method(self, common_kwargs, ci_method):
+ def _confidence_interval_method(self, common_kwargs, ci_method, return_period):
method_name = ci_method_to_method_name[ci_method]
mle_ci_parameters = {
'method': method_name,
@@ -33,7 +33,7 @@ class ResultFromMleExtremes(AbstractResultFromExtremes):
lower, mean_estimate, upper, _ = a
else:
d = self.get_python_dictionary(res)
- keys = ['50-year return level', '95% lower CI', '95% upper CI']
+ keys = ['{}-year return level'.format(return_period), '95% lower CI', '95% upper CI']
mean_estimate, lower, upper = [np.array(d[k])[0] for k in keys]
return mean_estimate, (lower, upper)
--
GitLab