diff --git a/extreme_data/meteo_france_data/scm_models_data/abstract_study.py b/extreme_data/meteo_france_data/scm_models_data/abstract_study.py
index 9376dc8c6f680dc22c0ed24aca238869726cef5e..feaef969992ea98f44d21b7286a1d04133a73e29 100644
--- a/extreme_data/meteo_france_data/scm_models_data/abstract_study.py
+++ b/extreme_data/meteo_france_data/scm_models_data/abstract_study.py
@@ -100,8 +100,8 @@ class AbstractStudy(object):
assert slope in SLOPES
self.orientation = orientation
self.slope = slope
- # Add some options
- self.fit_method = MarginFitMethod.is_mev_gev_fit
+ # Add some cache for computation
+ self._cache_for_pointwise_fit = {}
""" Time """
@@ -878,11 +878,17 @@ class AbstractStudy(object):
d, _ = self._massif_name_to_stationary_gev_params_and_confidence(quantile_level=None)
return d
- @cached_property
- def massif_name_to_stationary_gev_params_and_confidence_for_return_level_100(self):
- return self._massif_name_to_stationary_gev_params_and_confidence(quantile_level=0.99)
+ def massif_name_to_stationary_gev_params_and_confidence(self, quantile_level,
+ confidence_interval_based_on_delta_method):
+ t = (quantile_level, confidence_interval_based_on_delta_method)
+ if t in self._cache_for_pointwise_fit:
+ return self._cache_for_pointwise_fit[t]
+ else:
+ res = self._massif_name_to_stationary_gev_params_and_confidence(*t)
+ self._cache_for_pointwise_fit[t] = res
+ return res
- def _massif_name_to_stationary_gev_params_and_confidence(self, quantile_level=None):
+ def _massif_name_to_stationary_gev_params_and_confidence(self, quantile_level=None, confidence_interval_based_on_delta_method=True):
""" at least 90% of values must be above zero"""
massif_name_to_stationary_gev_params = {}
massif_name_to_confidence = {}
@@ -891,7 +897,8 @@ class AbstractStudy(object):
if percentage_of_non_zeros > 90:
gev_params, confidence = fitted_stationary_gev_with_uncertainty_interval(annual_maxima,
fit_method=MarginFitMethod.extremes_fevd_mle,
- quantile_level=quantile_level)
+ quantile_level=quantile_level,
+ confidence_interval_based_on_delta_method=confidence_interval_based_on_delta_method)
if -0.5 <= gev_params.shape <= 0.5:
massif_name_to_stationary_gev_params[massif_name] = gev_params
massif_name_to_confidence[massif_name] = confidence
diff --git a/extreme_data/meteo_france_data/scm_models_data/utils_function.py b/extreme_data/meteo_france_data/scm_models_data/utils_function.py
index ed545b2f6c1fc7351e05f62ee78b787d7046ae0f..58b2554094955c9dcfe979a1c3e362f9bd9276c5 100644
--- a/extreme_data/meteo_france_data/scm_models_data/utils_function.py
+++ b/extreme_data/meteo_france_data/scm_models_data/utils_function.py
@@ -1,22 +1,51 @@
+from multiprocessing import Pool
+
+import numpy as np
+from sklearn.utils import resample
+
from extreme_fit.estimator.margin_estimator.utils import _fitted_stationary_gev
from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import StationaryTemporalModel
from extreme_fit.model.margin_model.utils import MarginFitMethod
+from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_extract_eurocode_return_level import \
+ AbstractExtractEurocodeReturnLevel
from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import \
ConfidenceIntervalMethodFromExtremes
from extreme_fit.model.result_from_model_fit.result_from_extremes.eurocode_return_level_uncertainties import \
EurocodeConfidenceIntervalFromExtremes
+from root_utils import NB_CORES
def fitted_stationary_gev_with_uncertainty_interval(x_gev, fit_method=MarginFitMethod.is_mev_gev_fit,
model_class=StationaryTemporalModel,
starting_year=None,
- quantile_level=0.98):
+ quantile_level=0.98,
+ confidence_interval_based_on_delta_method=True):
estimator, gev_param = _fitted_stationary_gev(fit_method, model_class, starting_year, x_gev)
if quantile_level is not None:
EurocodeConfidenceIntervalFromExtremes.quantile_level = quantile_level
coordinate = estimator.dataset.coordinates.df_all_coordinates.iloc[0].values
- confidence_interval = EurocodeConfidenceIntervalFromExtremes.from_estimator_extremes(estimator, ConfidenceIntervalMethodFromExtremes.ci_mle,
- coordinate)
+ if confidence_interval_based_on_delta_method:
+ confidence_interval = EurocodeConfidenceIntervalFromExtremes.from_estimator_extremes(estimator, ConfidenceIntervalMethodFromExtremes.ci_mle,
+ coordinate)
+ else:
+ # Bootstrap method
+ nb_bootstrap = AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP
+ x_gev_list = [resample(x_gev) for _ in range(nb_bootstrap)]
+ arguments = [(fit_method, model_class, starting_year, x, quantile_level) for x in x_gev_list]
+ multiprocess = True
+ if multiprocess:
+ with Pool(NB_CORES) as p:
+ return_level_list = p.map(_compute_return_level, arguments)
+ else:
+ return_level_list = [_compute_return_level(argument) for argument in arguments]
+ confidence_interval = tuple([np.quantile(return_level_list, q)
+ for q in AbstractExtractEurocodeReturnLevel.bottom_and_upper_quantile])
+ mean_estimate = gev_param.quantile(quantile_level)
+ confidence_interval = EurocodeConfidenceIntervalFromExtremes(mean_estimate, confidence_interval)
else:
confidence_interval = None
return gev_param, confidence_interval
+
+def _compute_return_level(t):
+ fit_method, model_class, starting_year, x, quantile_level = t
+ return _fitted_stationary_gev(fit_method, model_class, starting_year, x)[1].quantile(quantile_level)
\ No newline at end of file
diff --git a/extreme_fit/estimator/margin_estimator/abstract_margin_estimator.py b/extreme_fit/estimator/margin_estimator/abstract_margin_estimator.py
index f4791e07d45dd8e23303b6241191fe88ce2aeefd..f233fa0ae67efe0cd7d16f5984298ccdc6c8c93c 100644
--- a/extreme_fit/estimator/margin_estimator/abstract_margin_estimator.py
+++ b/extreme_fit/estimator/margin_estimator/abstract_margin_estimator.py
@@ -88,6 +88,18 @@ class LinearMarginEstimator(AbstractMarginEstimator):
sorted_empirical_quantiles = sorted(sorted_empirical_quantiles)
return sorted_empirical_quantiles
+ def coordinate_values_to_maxima_from_standard_gumbel_quantiles(self, standard_gumbel_quantiles, split=Split.all):
+ coordinate_values_to_maxima = {}
+ coordinate_values = self.dataset.df_coordinates(split=split).values
+ assert len(standard_gumbel_quantiles) == len(coordinate_values)
+ for quantile, coordinate in zip(standard_gumbel_quantiles, coordinate_values):
+ gev_param = self.function_from_fit.get_params(
+ coordinate=coordinate,
+ is_transformed=False)
+ maximum = gev_param.gumbel_inverse_standardization(quantile)
+ coordinate_values_to_maxima[tuple(coordinate)] = np.array([maximum])
+ return coordinate_values_to_maxima
+
def deviance(self, split=Split.all):
return 2 * self.nllh(split=split)
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 aa9b8afb76d6f61c76623a4041a7a546b7b750a6..ac2c6d9d213698a7d890f57120677dce18c74b43 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,6 +16,13 @@ from root_utils import classproperty
class AbstractExtractEurocodeReturnLevel(object):
ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY = 0.05
+ NB_BOOTSTRAP = 1000
+
+ @classproperty
+ def bottom_and_upper_quantile(cls):
+ bottom_quantile = cls.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY / 2
+ bottom_and_upper_quantile = (bottom_quantile, 1 - bottom_quantile)
+ return bottom_and_upper_quantile
def __init__(self, estimator: LinearMarginEstimator, ci_method, temporal_covariate, quantile_level=EUROCODE_QUANTILE):
self.ci_method = ci_method
@@ -100,7 +107,5 @@ class ExtractEurocodeReturnLevelFromMyBayesianExtremes(AbstractExtractEurocodeRe
@property
def confidence_interval(self):
# Bottom and upper quantile correspond to the quantile
- bottom_quantile = self.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY / 2
- bottom_and_upper_quantile = (bottom_quantile, 1 - bottom_quantile)
return [np.quantile(self.posterior_eurocode_return_level_samples_for_temporal_covariate, q=q)
- for q in bottom_and_upper_quantile]
+ for q in self.bottom_and_upper_quantile]
diff --git a/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py b/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
index be76e1c31dab2037ac1a1fc18282c19284efe92b..f0b704dfcd5e009ea0bb96e7c8a36aa3fbf037c7 100644
--- a/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
+++ b/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
@@ -4,6 +4,7 @@ from typing import List
import matplotlib as mpl
from extreme_data.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal, CrocusSnowLoad3Days
+from extreme_fit.model.utils import set_seed_for_test
from projects.altitude_spatial_model.altitudes_fit.plots.plot_coherence_curves import plot_coherence_curves
from projects.altitude_spatial_model.altitudes_fit.plots.plot_histogram_altitude_studies import \
plot_histogram_all_models_against_altitudes, plot_histogram_all_trends_against_altitudes, \
@@ -26,13 +27,15 @@ def main():
study_classes = [SafranSnowfall1Day, SafranSnowfall3Days, SafranSnowfall5Days, SafranSnowfall7Days][:1]
seasons = [Season.annual, Season.winter, Season.spring, Season.automn][:1]
- fast = True
+ set_seed_for_test()
+
+ fast = None
if fast is None:
- massif_names = None
- altitudes_list = altitudes_for_groups[1:2]
- elif fast:
massif_names = ['Vanoise', 'Haute-Maurienne', 'Vercors'][:1]
altitudes_list = altitudes_for_groups[:]
+ elif fast:
+ massif_names = ['Vanoise', 'Haute-Maurienne', 'Vercors'][:1]
+ altitudes_list = altitudes_for_groups[3:]
else:
massif_names = None
altitudes_list = altitudes_for_groups[:]
@@ -46,7 +49,8 @@ def main_loop(altitudes_list, massif_names, seasons, study_classes):
for season in seasons:
for study_class in study_classes:
print('Inner loop', season, study_class)
- visualizer_list = load_visualizer_list(season, study_class, altitudes_list, massif_names)
+ visualizer_list = load_visualizer_list(season, study_class, altitudes_list, massif_names
+ )
plot_visualizers(massif_names, visualizer_list)
for visualizer in visualizer_list:
plot_visualizer(massif_names, visualizer)
diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
index 8b7f57f6e0fec896c9a6e71fe4ed4bb8484c536d..f41da81eb03b1d8c973f7ad56e5a75b626ef4372 100644
--- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
+++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
@@ -40,16 +40,18 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
fit_method=MarginFitMethod.extremes_fevd_mle,
temporal_covariate_for_fit=None,
display_only_model_that_pass_anderson_test=True,
+ confidence_interval_based_on_delta_method=False
):
super().__init__(studies.study, show=show, save_to_file=not show)
self.studies = studies
self.non_stationary_models = model_classes
self.fit_method = fit_method
self.temporal_covariate_for_fit = temporal_covariate_for_fit
- self.display_only_model_that_pass_anderson_test = display_only_model_that_pass_anderson_test
+ self.display_only_model_that_pass_test = display_only_model_that_pass_anderson_test
self.massif_names = massif_names if massif_names is not None else self.study.all_massif_names()
self.massif_name_to_massif_id = {m: i for i, m in enumerate(self.massif_names)}
self.altitude_group = get_altitude_group_from_altitudes(self.studies.altitudes)
+ self.confidence_interval_based_on_delta_method = confidence_interval_based_on_delta_method
# Load one fold fit
self.massif_name_to_massif_altitudes = {}
self._massif_name_to_one_fold_fit = {}
@@ -64,7 +66,8 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
old_fold_fit = OneFoldFit(massif_name, dataset, model_classes, self.fit_method,
self.temporal_covariate_for_fit,
type(self.altitude_group),
- self.display_only_model_that_pass_anderson_test)
+ self.display_only_model_that_pass_test,
+ self.confidence_interval_based_on_delta_method)
self._massif_name_to_one_fold_fit[massif_name] = old_fold_fit
# Print number of massif without any validated fit
massifs_without_any_validated_fit = [massif_name
@@ -103,7 +106,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
@property
def massif_name_to_one_fold_fit(self) -> Dict[str, OneFoldFit]:
return {massif_name: old_fold_fit for massif_name, old_fold_fit in self._massif_name_to_one_fold_fit.items()
- if not self.display_only_model_that_pass_anderson_test
+ if not self.display_only_model_that_pass_test
or old_fold_fit.has_at_least_one_valid_model}
def plot_moments(self):
diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py
index 35c892f12500108b5ba88c723ec94becc3a40ebb..82a1d73830059d01c6463d621dd6e982fdc588d8 100644
--- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py
+++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py
@@ -1,8 +1,11 @@
+from multiprocessing import Pool
+
import numpy.testing as npt
import numpy as np
import rpy2
from cached_property import cached_property
from scipy.stats import chi2
+from sklearn.utils import resample
from extreme_fit.distribution.gev.gev_params import GevParams
from extreme_fit.distribution.gumbel.gumbel_gof import goodness_of_fit_anderson
@@ -18,14 +21,18 @@ from extreme_fit.model.margin_model.polynomial_margin_model.gumbel_altitudinal_m
from extreme_fit.model.margin_model.polynomial_margin_model.models_based_on_pariwise_analysis.gev_with_linear_shape_wrt_altitude import \
AltitudinalShapeLinearTimeStationary
from extreme_fit.model.margin_model.utils import MarginFitMethod
+from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_extract_eurocode_return_level import \
+ AbstractExtractEurocodeReturnLevel
from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import \
ConfidenceIntervalMethodFromExtremes
from extreme_fit.model.result_from_model_fit.result_from_extremes.eurocode_return_level_uncertainties import \
EurocodeConfidenceIntervalFromExtremes
from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitude_group import AbstractAltitudeGroup, \
DefaultAltitudeGroup
-from root_utils import classproperty
+from root_utils import classproperty, NB_CORES
from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
+from spatio_temporal_dataset.slicer.split import Split
+from spatio_temporal_dataset.spatio_temporal_observations.annual_maxima_observations import AnnualMaxima
class OneFoldFit(object):
@@ -36,11 +43,14 @@ class OneFoldFit(object):
nb_years = 60
def __init__(self, massif_name: str, dataset: AbstractDataset, models_classes,
- fit_method=MarginFitMethod.extremes_fevd_mle, temporal_covariate_for_fit=None,
+ fit_method=MarginFitMethod.extremes_fevd_mle,
+ temporal_covariate_for_fit=None,
altitude_class=DefaultAltitudeGroup,
- only_models_that_pass_anderson_test=True,
+ only_models_that_pass_goodness_of_fit_test=True,
+ confidence_interval_based_on_delta_method=False,
):
- self.only_models_that_pass_anderson_test = only_models_that_pass_anderson_test
+ self.confidence_interval_based_on_delta_method = confidence_interval_based_on_delta_method
+ self.only_models_that_pass_goodness_of_fit_test = only_models_that_pass_goodness_of_fit_test
self.altitude_group = altitude_class()
self.massif_name = massif_name
self.dataset = dataset
@@ -51,16 +61,19 @@ class OneFoldFit(object):
# Fit Estimators
self.model_class_to_estimator = {}
for model_class in models_classes:
- self.model_class_to_estimator[model_class] = fitted_linear_margin_estimator_short(model_class=model_class,
- dataset=self.dataset,
- fit_method=self.fit_method,
- temporal_covariate_for_fit=self.temporal_covariate_for_fit)
+ self.model_class_to_estimator[model_class] = self.fitted_linear_margin_estimator(model_class, self.dataset)
# Best estimator definition
self.best_estimator_class_for_total_aic = None
# Cached object
self._folder_to_goodness = {}
+ def fitted_linear_margin_estimator(self, model_class, dataset):
+ return fitted_linear_margin_estimator_short(model_class=model_class,
+ dataset=dataset,
+ fit_method=self.fit_method,
+ temporal_covariate_for_fit=self.temporal_covariate_for_fit)
+
@classproperty
def folder_for_plots(cls):
return 'Total aic/' if cls.best_estimator_minimizes_total_aic else ''
@@ -129,7 +142,7 @@ class OneFoldFit(object):
@cached_property
def sorted_estimators_with_stationary(self):
- if self.only_models_that_pass_anderson_test:
+ if self.only_models_that_pass_goodness_of_fit_test:
return [e for e in self.sorted_estimators
if self.goodness_of_fit_test(e)
# and self.sensitivity_of_fit_test_top_maxima(e)
@@ -148,7 +161,7 @@ class OneFoldFit(object):
@cached_property
def sorted_estimators_without_stationary(self):
- if self.only_models_that_pass_anderson_test:
+ if self.only_models_that_pass_goodness_of_fit_test:
return [e for e in self._sorted_estimators_without_stationary if self.goodness_of_fit_test(e)]
else:
return self._sorted_estimators_without_stationary
@@ -240,13 +253,22 @@ class OneFoldFit(object):
@property
def is_significant(self) -> bool:
- stationary_model_classes = [StationaryAltitudinal, StationaryGumbelAltitudinal,
- AltitudinalShapeLinearTimeStationary]
- if any([isinstance(self.best_estimator.margin_model, c)
- for c in stationary_model_classes]):
- return False
+ if self.confidence_interval_based_on_delta_method:
+ stationary_model_classes = [StationaryAltitudinal, StationaryGumbelAltitudinal,
+ AltitudinalShapeLinearTimeStationary]
+ if any([isinstance(self.best_estimator.margin_model, c)
+ for c in stationary_model_classes]):
+ return False
+ else:
+ return self.likelihood_ratio > chi2.ppf(q=1 - self.SIGNIFICANCE_LEVEL, df=self.degree_freedom_chi2)
else:
- return self.likelihood_ratio > chi2.ppf(q=1 - self.SIGNIFICANCE_LEVEL, df=self.degree_freedom_chi2)
+ # Bootstrap based significance
+ print('new significance is used')
+ sign_of_changes = [self.sign_of_change(f) for f in self.bootstrap_fitted_functions_from_fit]
+ if self.sign_of_change(self.best_function_from_fit) > 0:
+ return np.quantile(sign_of_changes, self.SIGNIFICANCE_LEVEL) > 0
+ else:
+ return np.quantile(sign_of_changes, 1 - self.SIGNIFICANCE_LEVEL) < 0
# @property
# def goodness_of_fit_anderson_test(self):
@@ -272,7 +294,7 @@ class OneFoldFit(object):
fit_method=self.fit_method,
temporal_covariate_for_fit=self.temporal_covariate_for_fit)
# Compare sign of change
- has_not_opposite_sign = self.sign_of_change(estimator) * self.sign_of_change(new_estimator) >= 0
+ has_not_opposite_sign = self.sign_of_change(estimator.function_from_fit) * self.sign_of_change(new_estimator.function_from_fit) >= 0
return has_not_opposite_sign
def sensitivity_of_fit_test_last_years(self, estimator: LinearMarginEstimator):
@@ -287,16 +309,16 @@ class OneFoldFit(object):
fit_method=self.fit_method,
temporal_covariate_for_fit=self.temporal_covariate_for_fit)
# Compare sign of change
- has_not_opposite_sign = self.sign_of_change(estimator) * self.sign_of_change(new_estimator) >= 0
+ has_not_opposite_sign = self.sign_of_change(estimator.function_from_fit) * self.sign_of_change(new_estimator.function_from_fit) >= 0
# if not has_not_opposite_sign:
# print('Last years', self.massif_name, model_class, self.sign_of_change(estimator), self.sign_of_change(new_estimator))
return has_not_opposite_sign
- def sign_of_change(self, estimator):
+ def sign_of_change(self, function_from_fit):
return_levels = []
for year in [2019 - self.nb_years, 2019]:
coordinate = np.array([self.altitude_plot, year])
- return_level = estimator.function_from_fit.get_params(
+ return_level = function_from_fit.get_params(
coordinate=coordinate,
is_transformed=False).return_level(return_period=self.return_period)
return_levels.append(return_level)
@@ -314,7 +336,54 @@ class OneFoldFit(object):
return standard_gumbel_quantiles
def best_confidence_interval(self, altitude, year) -> EurocodeConfidenceIntervalFromExtremes:
- EurocodeConfidenceIntervalFromExtremes.quantile_level = self.quantile_level
- return EurocodeConfidenceIntervalFromExtremes.from_estimator_extremes(estimator_extremes=self.best_estimator,
- ci_method=ConfidenceIntervalMethodFromExtremes.ci_mle,
- coordinate=np.array([altitude, year]))
+ coordinate = np.array([altitude, year])
+ if self.confidence_interval_based_on_delta_method:
+ EurocodeConfidenceIntervalFromExtremes.quantile_level = self.quantile_level
+ return EurocodeConfidenceIntervalFromExtremes.from_estimator_extremes(
+ estimator_extremes=self.best_estimator,
+ ci_method=ConfidenceIntervalMethodFromExtremes.ci_mle,
+ coordinate=coordinate)
+ else:
+ print('nb of bootstrap for confidence interval=', AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP)
+ mean_estimate = self.get_return_level(self.best_function_from_fit, coordinate)
+ bootstrap_return_levels = [self.get_return_level(f, coordinate) for f in self.bootstrap_fitted_functions_from_fit]
+ print(mean_estimate, bootstrap_return_levels)
+ confidence_interval = tuple([np.quantile(bootstrap_return_levels, q)
+ for q in AbstractExtractEurocodeReturnLevel.bottom_and_upper_quantile])
+ print(mean_estimate, confidence_interval)
+ return EurocodeConfidenceIntervalFromExtremes(mean_estimate, confidence_interval)
+
+ def get_return_level(self, function_from_fit, coordinate):
+ return function_from_fit.get_params(coordinate).return_level(self.return_period)
+
+ @property
+ def bootstrap_data(self):
+ bootstrap = []
+ for _ in range(AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP):
+ residuals = self.best_estimator.sorted_empirical_standard_gumbel_quantiles(split=Split.all)
+ resample_residuals = resample(residuals)
+ coordinate_values_to_maxima = self.best_estimator. \
+ coordinate_values_to_maxima_from_standard_gumbel_quantiles(standard_gumbel_quantiles=resample_residuals)
+ bootstrap.append(coordinate_values_to_maxima)
+
+ return bootstrap
+
+ @cached_property
+ def bootstrap_fitted_functions_from_fit(self):
+ multiprocess = True
+ if multiprocess:
+ with Pool(NB_CORES) as p:
+ functions_from_fit = p.map(self.fit_one_bootstrap_estimator, self.bootstrap_data)
+ else:
+ functions_from_fit = []
+ for coordinate_values_to_maxima in self.bootstrap_data:
+ estimator = self.fit_one_bootstrap_estimator(coordinate_values_to_maxima)
+ functions_from_fit.append(estimator)
+ return functions_from_fit
+
+ def fit_one_bootstrap_estimator(self, coordinate_values_to_maxima):
+ coordinates = self.dataset.coordinates
+ observations = AnnualMaxima.from_coordinates(coordinates, coordinate_values_to_maxima)
+ dataset = AbstractDataset(observations=observations, coordinates=coordinates)
+ model_class = type(self.best_margin_model)
+ return self.fitted_linear_margin_estimator(model_class, dataset).function_from_fit
diff --git a/projects/altitude_spatial_model/altitudes_fit/plots/plot_coherence_curves.py b/projects/altitude_spatial_model/altitudes_fit/plots/plot_coherence_curves.py
index 7e47ac6edc0622b8bd379672c8ec6957ae0b6f27..9b7c74c7915bdacf1285bacb46df1024446e1148 100644
--- a/projects/altitude_spatial_model/altitudes_fit/plots/plot_coherence_curves.py
+++ b/projects/altitude_spatial_model/altitudes_fit/plots/plot_coherence_curves.py
@@ -126,23 +126,32 @@ def load_all_list(massif_name, visualizer_list, altitudinal_model=True, year=201
if massif_name in visualizer.massif_name_to_one_fold_fit:
if altitudinal_model:
+ # Piecewise altitudinal temporal model
min_altitude, *_, max_altitude = visualizer.massif_name_to_massif_altitudes[massif_name]
one_fold_fit = visualizer.massif_name_to_one_fold_fit[massif_name]
- altitudes_list = list(range(min_altitude, max_altitude, 10))
+ step = 300
+ altitudes_list = list(range(min_altitude, max_altitude + step, step))
gev_params_list = [one_fold_fit.get_gev_params(altitude, year) for altitude in altitudes_list]
confidence_interval_values = [one_fold_fit.best_confidence_interval(altitude, year) for altitude in
altitudes_list]
else:
+ # Pointwise distribution
assert OneFoldFit.return_period == 100, 'change the call below'
altitudes_list, study_list_valid = zip(*[(a, s) for a, s in visualizer.studies.altitude_to_study.items()
if massif_name in
- s.massif_name_to_stationary_gev_params_and_confidence_for_return_level_100[
+ s.massif_name_to_stationary_gev_params_and_confidence
+ (OneFoldFit.quantile_level,
+ visualizer.confidence_interval_based_on_delta_method)[
0]])
gev_params_list = [
- study.massif_name_to_stationary_gev_params_and_confidence_for_return_level_100[0][massif_name]
+ study.massif_name_to_stationary_gev_params_and_confidence
+ (OneFoldFit.quantile_level,
+ visualizer.confidence_interval_based_on_delta_method)[0][massif_name]
for study in study_list_valid]
confidence_interval_values = [
- study.massif_name_to_stationary_gev_params_and_confidence_for_return_level_100[1][massif_name]
+ study.massif_name_to_stationary_gev_params_and_confidence
+ (OneFoldFit.quantile_level,
+ visualizer.confidence_interval_based_on_delta_method)[1][massif_name]
for study in study_list_valid]
# Checks
diff --git a/projects/altitude_spatial_model/altitudes_fit/utils_altitude_studies_visualizer.py b/projects/altitude_spatial_model/altitudes_fit/utils_altitude_studies_visualizer.py
index 3eaf91fece20944a2fda1dd21d8fb637bbc8eb65..9ce6d3825a4c8f521016fcbbabe578c3bc835f53 100644
--- a/projects/altitude_spatial_model/altitudes_fit/utils_altitude_studies_visualizer.py
+++ b/projects/altitude_spatial_model/altitudes_fit/utils_altitude_studies_visualizer.py
@@ -16,6 +16,7 @@ def load_visualizer_list(season, study_class, altitudes_list, massif_names, **kw
massif_names=massif_names,
show=False,
temporal_covariate_for_fit=None,
+ confidence_interval_based_on_delta_method=False
)
visualizer_list.append(visualizer)
compute_and_assign_max_abs(visualizer_list)