diff --git a/experiment/eurocode_data/eurocode_return_level_uncertainties.py b/experiment/eurocode_data/eurocode_return_level_uncertainties.py
deleted file mode 100644
index c81824a5b574e680afb55d17555d7fbf07fcb67f..0000000000000000000000000000000000000000
--- a/experiment/eurocode_data/eurocode_return_level_uncertainties.py
+++ /dev/null
@@ -1,116 +0,0 @@
-from enum import Enum
-from typing import List
-
-import numpy as np
-from cached_property import cached_property
-
-from experiment.eurocode_data.utils import EUROCODE_QUANTILE, YEAR_OF_INTEREST_FOR_RETURN_LEVEL
-from experiment.trend_analysis.univariate_test.utils import load_temporal_coordinates_and_dataset, \
- fitted_linear_margin_estimator
-from extreme_fit.distribution.gev.gev_params import GevParams
-from extreme_fit.estimator.abstract_estimator import AbstractEstimator
-from extreme_fit.estimator.margin_estimator.abstract_margin_estimator import LinearMarginEstimator
-from extreme_fit.estimator.utils import load_margin_function
-from extreme_fit.model.margin_model.linear_margin_model.abstract_temporal_linear_margin_model import \
- AbstractTemporalLinearMarginModel, TemporalMarginFitMethod
-from extreme_fit.model.margin_model.margin_function.linear_margin_function import LinearMarginFunction
-from extreme_fit.model.result_from_model_fit.result_from_extremes import ResultFromExtremes, ResultFromBayesianExtremes
-
-
-class ConfidenceIntervalMethodFromExtremes(Enum):
- # Confidence interval from the ci function
- bayes = 0
- normal = 1
- boot = 2
- proflik = 3
- # Confidence interval from my functions
- my_bayes = 4
-
-
-def compute_eurocode_level_uncertainty(last_year_for_the_data, smooth_maxima_x_y, model_class, ci_method):
- years, smooth_maxima = smooth_maxima_x_y
- idx = years.index(last_year_for_the_data) + 1
- years, smooth_maxima = years[:idx], smooth_maxima[:idx]
- return EurocodeLevelUncertaintyFromExtremes.from_maxima_years_model_class(smooth_maxima, years, model_class, ci_method)
-
-
-class EurocodeLevelUncertaintyFromExtremes(object):
- YEAR_OF_INTEREST = 2017
-
- def __init__(self, posterior_mean, poster_uncertainty_interval):
- self.posterior_mean = posterior_mean
- self.poster_uncertainty_interval = poster_uncertainty_interval
-
- @classmethod
- def from_estimator_extremes(cls, estimator_extremes: LinearMarginEstimator,
- ci_method: ConfidenceIntervalMethodFromExtremes):
- extractor = ExtractEurocodeReturnLevelFromExtremes(estimator_extremes, ci_method, cls.YEAR_OF_INTEREST)
- return cls(extractor.posterior_mean_eurocode_return_level_for_the_year_of_interest,
- extractor.posterior_eurocode_return_level_uncertainty_interval_for_the_year_of_interest)
-
- @classmethod
- def from_maxima_years_model_class(cls, maxima, years, model_class,
- ci_method=ConfidenceIntervalMethodFromExtremes.bayes):
- # Load coordinates and dataset
- coordinates, dataset = load_temporal_coordinates_and_dataset(maxima, years)
- # Select fit method depending on the ci_method
- if ci_method in [ConfidenceIntervalMethodFromExtremes.bayes,
- ConfidenceIntervalMethodFromExtremes.my_bayes]:
- fit_method = TemporalMarginFitMethod.extremes_fevd_bayesian
- else:
- fit_method = TemporalMarginFitMethod.extremes_fevd_mle
- # Fitted estimator
- fitted_estimator = fitted_linear_margin_estimator(model_class, coordinates, dataset, starting_year=1958,
- fit_method=fit_method)
- # Load object from result from extremes
- return cls.from_estimator_extremes(fitted_estimator, ci_method)
-
-
-class ExtractFromExtremes(object):
- pass
-
-
-class ExtractEurocodeReturnLevelFromExtremes(object):
- ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY = 0.05
-
- def __init__(self, estimator: LinearMarginEstimator,
- ci_method,
- year_of_interest: int = YEAR_OF_INTEREST_FOR_RETURN_LEVEL,
- alpha_for_confidence_interval: int = ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY,
- ):
- self.estimator = estimator
- self.result_from_fit = self.estimator.result_from_model_fit # type: ResultFromBayesianExtremes
- assert isinstance(self.result_from_fit, ResultFromBayesianExtremes)
- self.year_of_interest = year_of_interest
- self.alpha_for_confidence_interval = alpha_for_confidence_interval
-
- @property
- def margin_functions_from_fit(self) -> List[LinearMarginFunction]:
- margin_functions = []
- for _, s in self.result_from_fit.df_posterior_samples.iterrows():
- coef_dict = self.result_from_fit.get_coef_dict_from_posterior_sample(s)
- margin_function = load_margin_function(self.estimator, self.estimator.margin_model, coef_dict=coef_dict)
- margin_functions.append(margin_function)
- return margin_functions
-
- @property
- def gev_params_from_fit_for_year_of_interest(self) -> List[GevParams]:
- return [margin_function.get_gev_params(coordinate=np.array([self.year_of_interest]), is_transformed=False)
- for margin_function in self.margin_functions_from_fit]
-
- @cached_property
- def posterior_eurocode_return_level_samples_for_year_of_interest(self) -> np.ndarray:
- """We divide by 100 to transform the snow water equivalent into snow load"""
- return np.array([p.quantile(EUROCODE_QUANTILE) for p in self.gev_params_from_fit_for_year_of_interest]) / 100
-
- @property
- def posterior_mean_eurocode_return_level_for_the_year_of_interest(self) -> np.ndarray:
- return np.mean(self.posterior_eurocode_return_level_samples_for_year_of_interest)
-
- @property
- def posterior_eurocode_return_level_uncertainty_interval_for_the_year_of_interest(self):
- # Bottom and upper quantile correspond to the quantile
- bottom_quantile = self.alpha_for_confidence_interval / 2
- bottom_and_upper_quantile = (bottom_quantile, 1 - bottom_quantile)
- return [np.quantile(self.posterior_eurocode_return_level_samples_for_year_of_interest, q=q)
- for q in bottom_and_upper_quantile]
diff --git a/experiment/eurocode_data/eurocode_visualizer.py b/experiment/eurocode_data/eurocode_visualizer.py
index eb23bfdceccd5ab9708409f815611733b4a98a37..9b10db963eaf980e09794e5bbef337f8c505083c 100644
--- a/experiment/eurocode_data/eurocode_visualizer.py
+++ b/experiment/eurocode_data/eurocode_visualizer.py
@@ -2,9 +2,8 @@ from typing import Dict, List
import matplotlib.pyplot as plt
-from experiment.eurocode_data.eurocode_return_level_uncertainties import EurocodeLevelUncertaintyFromExtremes
-from experiment.eurocode_data.massif_name_to_departement import DEPARTEMENT_TYPES, massif_name_to_eurocode_region
-from experiment.eurocode_data.utils import EUROCODE_QUANTILE, EUROCODE_ALTITUDES
+from extreme_fit.model.result_from_model_fit.result_from_extremes.eurocode_return_level_uncertainties import EurocodeConfidenceIntervalFromExtremes
+from experiment.eurocode_data.massif_name_to_departement import massif_name_to_eurocode_region
from experiment.meteo_france_data.scm_models_data.visualization.utils import create_adjusted_axes
from root_utils import get_display_name_from_object_type
@@ -55,7 +54,7 @@ def plot_model_name_to_uncertainty_method_to_ordered_dict(d, massif_name, axes):
def plot_label_to_ordered_return_level_uncertainties(ax, massif_name, model_name,
label_to_ordered_return_level_uncertainties:
Dict[str, List[
- EurocodeLevelUncertaintyFromExtremes]]):
+ EurocodeConfidenceIntervalFromExtremes]]):
""" Generic function that might be used by many other more global functions"""
colors = ['tab:blue', 'tab:orange', 'tab:purple', 'tab:olive']
alpha = 0.2
diff --git a/experiment/eurocode_data/main_eurocode_drawing.py b/experiment/eurocode_data/main_eurocode_drawing.py
index eda7cf77180b1e4c6f089a511806d8c548a910ff..054c66be68fa34dcba18212c192b736f509a9308 100644
--- a/experiment/eurocode_data/main_eurocode_drawing.py
+++ b/experiment/eurocode_data/main_eurocode_drawing.py
@@ -1,10 +1,11 @@
import time
from collections import OrderedDict
-from experiment.eurocode_data.eurocode_return_level_uncertainties import ConfidenceIntervalMethodFromExtremes
+from extreme_fit.model.result_from_model_fit.result_from_extremes.eurocode_return_level_uncertainties import \
+ ConfidenceIntervalMethodFromExtremes
from experiment.eurocode_data.eurocode_visualizer import \
plot_massif_name_to_model_name_to_uncertainty_method_to_ordered_dict, get_model_name
-from experiment.eurocode_data.massif_name_to_departement import DEPARTEMENT_TYPES, MASSIF_NAMES_ALPS
+from experiment.eurocode_data.massif_name_to_departement import MASSIF_NAMES_ALPS
from experiment.eurocode_data.utils import EUROCODE_ALTITUDES, LAST_YEAR_FOR_EUROCODE
from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSwe3Days
from experiment.meteo_france_data.scm_models_data.visualization.hypercube_visualization.altitude_hypercube_visualizer import \
@@ -13,7 +14,6 @@ from experiment.meteo_france_data.scm_models_data.visualization.hypercube_visual
load_altitude_visualizer
from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import StationaryTemporalModel, \
NonStationaryLocationAndScaleTemporalModel
-from root_utils import get_display_name_from_object_type
# Model class
import matplotlib as mpl
@@ -22,7 +22,8 @@ mpl.rcParams['text.usetex'] = True
mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
-def massif_name_to_ordered_return_level_uncertainties(model_class, last_year_for_the_data, altitudes, massif_names, uncertainty_methods):
+def massif_name_to_ordered_return_level_uncertainties(model_class, last_year_for_the_data, altitudes, massif_names,
+ uncertainty_methods):
# Load model name
model_name = get_model_name(model_class)
# Load altitude visualizer
@@ -35,20 +36,20 @@ def massif_name_to_ordered_return_level_uncertainties(model_class, last_year_for
trend_test_class=None) # type: AltitudeHypercubeVisualizer
# Loop on the data
assert isinstance(altitude_visualizer.tuple_to_study_visualizer, OrderedDict)
- massif_name_to_ordered_eurocode_level_uncertainty = {massif_name: {ci_method: [] for ci_method in uncertainty_methods} for massif_name in massif_names}
+ massif_name_to_ordered_eurocode_level_uncertainty = {
+ massif_name: {ci_method: [] for ci_method in uncertainty_methods} for massif_name in massif_names}
for altitude, visualizer in altitude_visualizer.tuple_to_study_visualizer.items():
print('{} processing altitude = {} '.format(model_name, altitude))
for ci_method in uncertainty_methods:
- d = visualizer.massif_name_to_altitude_and_eurocode_level_uncertainty(model_class, last_year_for_the_data, massif_names, ci_method)
+ d = visualizer.massif_name_to_altitude_and_eurocode_level_uncertainty(model_class, last_year_for_the_data,
+ massif_names, ci_method)
# Append the altitude one by one
for massif_name, return_level_uncertainty in d.items():
- massif_name_to_ordered_eurocode_level_uncertainty[massif_name][ci_method].append(return_level_uncertainty)
+ massif_name_to_ordered_eurocode_level_uncertainty[massif_name][ci_method].append(
+ return_level_uncertainty)
return {model_name: massif_name_to_ordered_eurocode_level_uncertainty}
-
-
-
def main_drawing():
fast_plot = [True, False][0]
# Select parameters
@@ -60,8 +61,10 @@ def main_drawing():
][1:]
altitudes = EUROCODE_ALTITUDES[:]
uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.my_bayes, ConfidenceIntervalMethodFromExtremes.bayes]
+ show = True
if fast_plot:
+ show = True
model_class_and_last_year = model_class_and_last_year[:1]
altitudes = altitudes[2:4]
# altitudes = altitudes[:]
@@ -72,7 +75,8 @@ def main_drawing():
for model_class, last_year_for_the_data in model_class_and_last_year:
start = time.time()
model_name_to_massif_name_to_ordered_return_level.update(
- massif_name_to_ordered_return_level_uncertainties(model_class, last_year_for_the_data, altitudes, massif_names, uncertainty_methods))
+ massif_name_to_ordered_return_level_uncertainties(model_class, last_year_for_the_data, altitudes,
+ massif_names, uncertainty_methods))
duration = time.time() - start
print(model_class, duration)
# Transform the dictionary into the desired format
@@ -84,7 +88,7 @@ def main_drawing():
# Plot graph
plot_massif_name_to_model_name_to_uncertainty_method_to_ordered_dict(
massif_name_to_model_name_to_ordered_return_level_uncertainties, nb_massif_names=len(massif_names),
- nb_model_names=len(model_class_and_last_year), show=True)
+ nb_model_names=len(model_class_and_last_year), show=show)
if __name__ == '__main__':
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 5c996b9025df3c3ce7187441edce9a62f7897948..35e3b9eb6f9f4f7e0857ff11acb3c3a5e93ac327 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
@@ -11,9 +11,8 @@ import numpy as np
import pandas as pd
import seaborn as sns
-from experiment.eurocode_data.eurocode_return_level_uncertainties import ExtractEurocodeReturnLevelFromExtremes, \
- EurocodeLevelUncertaintyFromExtremes, compute_eurocode_level_uncertainty
-from experiment.eurocode_data.massif_name_to_departement import dep_class_to_massif_names
+from extreme_fit.model.result_from_model_fit.result_from_extremes.eurocode_return_level_uncertainties import \
+ EurocodeConfidenceIntervalFromExtremes, compute_eurocode_confidence_interval
from experiment.meteo_france_data.scm_models_data.abstract_extended_study import AbstractExtendedStudy
from experiment.trend_analysis.abstract_score import MeanScore, AbstractTrendScore
from experiment.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
@@ -355,14 +354,14 @@ class StudyVisualizer(VisualizationParameters):
start_year, stop_year = self.study.start_year_and_stop_year
return list(range(start_year, stop_year))
- def massif_name_to_altitude_and_eurocode_level_uncertainty(self, model_class, last_year_for_the_data, massif_names, ci_method) -> Dict[str, Tuple[int, EurocodeLevelUncertaintyFromExtremes]]:
+ def massif_name_to_altitude_and_eurocode_level_uncertainty(self, model_class, last_year_for_the_data, massif_names, ci_method) -> Dict[str, Tuple[int, EurocodeConfidenceIntervalFromExtremes]]:
massif_ids_and_names = [(massif_id, massif_name) for massif_id, massif_name in enumerate(self.study.study_massif_names) if massif_name in massif_names]
arguments = [[last_year_for_the_data, self.smooth_maxima_x_y(massif_id), model_class, ci_method] for massif_id, _ in massif_ids_and_names]
if self.multiprocessing:
with Pool(NB_CORES) as p:
- res = p.starmap(compute_eurocode_level_uncertainty, arguments)
+ res = p.starmap(compute_eurocode_confidence_interval, arguments)
else:
- res = [compute_eurocode_level_uncertainty(*argument) for argument in arguments]
+ res = [compute_eurocode_confidence_interval(*argument) for argument in arguments]
res_and_altitude = [(self.study.altitude, r) for r in res]
massif_name_to_eurocode_return_level_uncertainty = OrderedDict(zip([massif_name for _, massif_name in massif_ids_and_names], res_and_altitude))
return massif_name_to_eurocode_return_level_uncertainty
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 bba4e48b1e3fffa4a68abb0ad6f136edc7cd7acf..851b9f5477e6bb11ae8414be4330b17cf5c73029 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
@@ -6,9 +6,10 @@ import pandas as pd
from extreme_fit.distribution.gev.gev_params import GevParams
from extreme_fit.model.margin_model.linear_margin_model.linear_margin_model import LinearMarginModel
from extreme_fit.model.result_from_model_fit.abstract_result_from_model_fit import AbstractResultFromModelFit
-from extreme_fit.model.result_from_model_fit.result_from_extremes import ResultFromExtremes, ResultFromBayesianExtremes
+from extreme_fit.model.result_from_model_fit.result_from_extremes.result_from_bayesian_extremes import AbstractResultFromExtremes, ResultFromBayesianExtremes
+from extreme_fit.model.result_from_model_fit.result_from_extremes.result_from_mle_extremes import ResultFromMleExtremes
from extreme_fit.model.result_from_model_fit.result_from_ismev import ResultFromIsmev
-from extreme_fit.model.utils import r, ro, get_null, get_margin_formula_extremes, get_coord, get_coord_df
+from extreme_fit.model.utils import r, ro, get_null, get_margin_formula_extremes, get_coord_df
from extreme_fit.model.utils import safe_run_r_estimator
from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
@@ -21,8 +22,6 @@ class TemporalMarginFitMethod(Enum):
class AbstractTemporalLinearMarginModel(LinearMarginModel):
"""Linearity only with respect to the temporal coordinates"""
- ISMEV_GEV_FIT_METHOD_STR = 'isMev.gev.fit'
- EXTREMES_FEVD_BAYESIAN_FIT_METHOD_STR = 'extRemes.fevd.Bayesian'
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):
@@ -52,7 +51,7 @@ class AbstractTemporalLinearMarginModel(LinearMarginModel):
# Gev fit with extRemes package
- def extremes_fevd_mle_fit(self, x, df_coordinates_temp) -> ResultFromExtremes:
+ def extremes_fevd_mle_fit(self, x, df_coordinates_temp) -> AbstractResultFromExtremes:
r_type_argument_kwargs, y = self.extreme_arguments(df_coordinates_temp)
res = safe_run_r_estimator(function=r('fevd_fixed'),
x=x,
@@ -60,9 +59,9 @@ class AbstractTemporalLinearMarginModel(LinearMarginModel):
method='MLE',
**r_type_argument_kwargs
)
- return ResultFromExtremes(res, self.margin_function_start_fit.gev_param_name_to_dims)
+ return ResultFromMleExtremes(res, self.margin_function_start_fit.gev_param_name_to_dims)
- def extremes_fevd_bayesian_fit(self, x, df_coordinates_temp) -> ResultFromExtremes:
+ def extremes_fevd_bayesian_fit(self, x, df_coordinates_temp) -> AbstractResultFromExtremes:
r_type_argument_kwargs, y = self.extreme_arguments(df_coordinates_temp)
# Assert for any non-stationary model that the shape parameter is constant
# (because the prior function considers that the last parameter should be the shape)
diff --git a/extreme_fit/model/result_from_model_fit/abstract_result_from_model_fit.py b/extreme_fit/model/result_from_model_fit/abstract_result_from_model_fit.py
index 3a91b1dfa6497cf9be9d8ee1a476062784a0f04f..acd280fb1735c4c6b8b11e311da8a5c8cdcaa6c7 100644
--- a/extreme_fit/model/result_from_model_fit/abstract_result_from_model_fit.py
+++ b/extreme_fit/model/result_from_model_fit/abstract_result_from_model_fit.py
@@ -7,10 +7,14 @@ class AbstractResultFromModelFit(object):
def __init__(self, result_from_fit: robjects.ListVector) -> None:
if hasattr(result_from_fit, 'names'):
- self.name_to_value = {name: result_from_fit.rx2(name) for name in result_from_fit.names}
+ self.name_to_value = self.get_python_dictionary(result_from_fit)
else:
self.name_to_value = {}
+ @staticmethod
+ def get_python_dictionary(r_dictionary):
+ return {name: r_dictionary.rx2(name) for name in r_dictionary.names}
+
@property
def names(self):
return self.name_to_value.keys()
diff --git a/extreme_fit/model/result_from_model_fit/result_from_extremes/__init__.py b/extreme_fit/model/result_from_model_fit/result_from_extremes/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
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
new file mode 100644
index 0000000000000000000000000000000000000000..4c0255e3385fb8fc3250b8803c4a39dd3eb051c1
--- /dev/null
+++ b/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_extract_eurocode_return_level.py
@@ -0,0 +1,74 @@
+from typing import List
+
+import numpy as np
+from cached_property import cached_property
+
+from experiment.eurocode_data.utils import EUROCODE_QUANTILE, YEAR_OF_INTEREST_FOR_RETURN_LEVEL
+from extreme_fit.distribution.gev.gev_params import GevParams
+from extreme_fit.estimator.margin_estimator.abstract_margin_estimator import LinearMarginEstimator
+from extreme_fit.estimator.utils import load_margin_function
+from extreme_fit.model.margin_model.margin_function.linear_margin_function import LinearMarginFunction
+from extreme_fit.model.result_from_model_fit.result_from_extremes.result_from_bayesian_extremes import \
+ ResultFromBayesianExtremes
+
+
+class AbstractExtractEurocodeReturnLevel(object):
+ ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY = 0.05
+
+ def __init__(self, estimator: LinearMarginEstimator,
+ ci_method,
+ year_of_interest: int = YEAR_OF_INTEREST_FOR_RETURN_LEVEL,
+ ):
+ self.ci_method = ci_method
+ self.estimator = estimator
+ self.result_from_fit = self.estimator.result_from_model_fit
+
+ self.year_of_interest = year_of_interest
+
+ self.eurocode_quantile = EUROCODE_QUANTILE
+ self.alpha_for_confidence_interval = self.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY
+
+
+class ExtractEurocodeReturnLevelFromCiMethod(AbstractExtractEurocodeReturnLevel):
+ pass
+
+
+class ExtractEurocodeReturnLevelFromMyBayesianExtremes(AbstractExtractEurocodeReturnLevel):
+ result_from_fit: ResultFromBayesianExtremes
+
+ def __init__(self, estimator: LinearMarginEstimator, ci_method,
+ year_of_interest: int = YEAR_OF_INTEREST_FOR_RETURN_LEVEL):
+ super().__init__(estimator, ci_method, year_of_interest)
+ assert isinstance(self.result_from_fit, ResultFromBayesianExtremes)
+
+ @property
+ def margin_functions_from_fit(self) -> List[LinearMarginFunction]:
+ margin_functions = []
+ for _, s in self.result_from_fit.df_posterior_samples.iterrows():
+ coef_dict = self.result_from_fit.get_coef_dict_from_posterior_sample(s)
+ margin_function = load_margin_function(self.estimator, self.estimator.margin_model, coef_dict=coef_dict)
+ margin_functions.append(margin_function)
+ return margin_functions
+
+ @property
+ def gev_params_from_fit_for_year_of_interest(self) -> List[GevParams]:
+ return [margin_function.get_gev_params(coordinate=np.array([self.year_of_interest]), is_transformed=False)
+ for margin_function in self.margin_functions_from_fit]
+
+ @cached_property
+ def posterior_eurocode_return_level_samples_for_year_of_interest(self) -> np.ndarray:
+ """We divide by 100 to transform the snow water equivalent into snow load"""
+ return np.array(
+ [p.quantile(self.eurocode_quantile) for p in self.gev_params_from_fit_for_year_of_interest]) / 100
+
+ @property
+ def posterior_mean_eurocode_return_level_for_the_year_of_interest(self) -> np.ndarray:
+ return np.mean(self.posterior_eurocode_return_level_samples_for_year_of_interest)
+
+ @property
+ def posterior_eurocode_return_level_uncertainty_interval_for_the_year_of_interest(self):
+ # Bottom and upper quantile correspond to the quantile
+ bottom_quantile = self.alpha_for_confidence_interval / 2
+ bottom_and_upper_quantile = (bottom_quantile, 1 - bottom_quantile)
+ return [np.quantile(self.posterior_eurocode_return_level_samples_for_year_of_interest, q=q)
+ for q in bottom_and_upper_quantile]
\ No newline at end of file
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
new file mode 100644
index 0000000000000000000000000000000000000000..17522c74f1aa62745596d861da9f3d12f50fa228
--- /dev/null
+++ b/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_result_from_extremes.py
@@ -0,0 +1,18 @@
+import numpy as np
+import pandas as pd
+from rpy2 import robjects
+
+from extreme_fit.model.result_from_model_fit.abstract_result_from_model_fit import \
+ AbstractResultFromModelFit
+from extreme_fit.model.utils import r
+
+
+class AbstractResultFromExtremes(AbstractResultFromModelFit):
+
+ def __init__(self, result_from_fit: robjects.ListVector, gev_param_name_to_dim=None) -> None:
+ super().__init__(result_from_fit)
+ self.gev_param_name_to_dim = gev_param_name_to_dim
+
+ def load_dataframe_from_r_matrix(self, name):
+ r_matrix = self.name_to_value[name]
+ return pd.DataFrame(np.array(r_matrix), columns=r.colnames(r_matrix))
diff --git a/extreme_fit/model/result_from_model_fit/result_from_extremes/confidence_interval_method.py b/extreme_fit/model/result_from_model_fit/result_from_extremes/confidence_interval_method.py
new file mode 100644
index 0000000000000000000000000000000000000000..f2f03c3ef1902c58bebcc1226f750227a5e170ca
--- /dev/null
+++ b/extreme_fit/model/result_from_model_fit/result_from_extremes/confidence_interval_method.py
@@ -0,0 +1,11 @@
+from enum import Enum
+
+
+class ConfidenceIntervalMethodFromExtremes(Enum):
+ # Confidence interval from the ci function
+ bayes = 0
+ normal = 1
+ boot = 2
+ proflik = 3
+ # Confidence interval from my functions
+ my_bayes = 4
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
new file mode 100644
index 0000000000000000000000000000000000000000..427801c3160c571ea5dc9c4af342796f6f368f55
--- /dev/null
+++ b/extreme_fit/model/result_from_model_fit/result_from_extremes/eurocode_return_level_uncertainties.py
@@ -0,0 +1,56 @@
+from enum import Enum
+
+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, \
+ fitted_linear_margin_estimator
+from extreme_fit.estimator.margin_estimator.abstract_margin_estimator import LinearMarginEstimator
+from extreme_fit.model.margin_model.linear_margin_model.abstract_temporal_linear_margin_model import \
+ TemporalMarginFitMethod
+from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import \
+ ConfidenceIntervalMethodFromExtremes
+
+
+def compute_eurocode_confidence_interval(last_year_for_the_data, smooth_maxima_x_y, model_class, ci_method):
+ years, smooth_maxima = smooth_maxima_x_y
+ idx = years.index(last_year_for_the_data) + 1
+ years, smooth_maxima = years[:idx], smooth_maxima[:idx]
+ return EurocodeConfidenceIntervalFromExtremes.from_maxima_years_model_class(smooth_maxima, years, model_class, ci_method)
+
+
+class EurocodeConfidenceIntervalFromExtremes(object):
+ YEAR_OF_INTEREST = 2017
+
+ def __init__(self, posterior_mean, poster_uncertainty_interval):
+ self.posterior_mean = posterior_mean
+ self.poster_uncertainty_interval = poster_uncertainty_interval
+
+ @classmethod
+ def from_estimator_extremes(cls, estimator_extremes: LinearMarginEstimator,
+ ci_method: ConfidenceIntervalMethodFromExtremes):
+ if ci_method == ConfidenceIntervalMethodFromExtremes.my_bayes:
+ extractor = ExtractEurocodeReturnLevelFromMyBayesianExtremes(estimator_extremes, ci_method, cls.YEAR_OF_INTEREST)
+ else:
+ extractor = ExtractEurocodeReturnLevelFromCiMethod(estimator_extremes, ci_method, cls.YEAR_OF_INTEREST)
+ return cls(extractor.posterior_mean_eurocode_return_level_for_the_year_of_interest,
+ extractor.posterior_eurocode_return_level_uncertainty_interval_for_the_year_of_interest)
+
+ @classmethod
+ def from_maxima_years_model_class(cls, maxima, years, model_class,
+ ci_method=ConfidenceIntervalMethodFromExtremes.bayes):
+ # Load coordinates and dataset
+ coordinates, dataset = load_temporal_coordinates_and_dataset(maxima, years)
+ # Select fit method depending on the ci_method
+ if ci_method in [ConfidenceIntervalMethodFromExtremes.bayes,
+ ConfidenceIntervalMethodFromExtremes.my_bayes]:
+ fit_method = TemporalMarginFitMethod.extremes_fevd_bayesian
+ else:
+ fit_method = TemporalMarginFitMethod.extremes_fevd_mle
+ # Fitted estimator
+ fitted_estimator = fitted_linear_margin_estimator(model_class, coordinates, dataset, starting_year=1958,
+ fit_method=fit_method)
+ # Load object from result from extremes
+ return cls.from_estimator_extremes(fitted_estimator, ci_method)
+
+
+
diff --git a/extreme_fit/model/result_from_model_fit/result_from_extremes.py b/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_bayesian_extremes.py
similarity index 68%
rename from extreme_fit/model/result_from_model_fit/result_from_extremes.py
rename to extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_bayesian_extremes.py
index 194a3f60c34e8b962041f3ed86543ff8c53d4b32..296a7e4fbd6f2984769e588585b6f9b2337bb6e4 100644
--- a/extreme_fit/model/result_from_model_fit/result_from_extremes.py
+++ b/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_bayesian_extremes.py
@@ -1,26 +1,12 @@
-import numpy as np
import pandas as pd
from rpy2 import robjects
-from extreme_fit.distribution.gev.gev_params import GevParams
-from extreme_fit.model.result_from_model_fit.abstract_result_from_model_fit import \
- AbstractResultFromModelFit
+from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_result_from_extremes import \
+ AbstractResultFromExtremes
from extreme_fit.model.result_from_model_fit.utils import get_margin_coef_ordered_dict
-from extreme_fit.model.utils import r
-class ResultFromExtremes(AbstractResultFromModelFit):
-
- def __init__(self, result_from_fit: robjects.ListVector, gev_param_name_to_dim=None) -> None:
- super().__init__(result_from_fit)
- self.gev_param_name_to_dim = gev_param_name_to_dim
-
- def load_dataframe_from_r_matrix(self, name):
- r_matrix = self.name_to_value[name]
- return pd.DataFrame(np.array(r_matrix), columns=r.colnames(r_matrix))
-
-
-class ResultFromBayesianExtremes(ResultFromExtremes):
+class ResultFromBayesianExtremes(AbstractResultFromExtremes):
def __init__(self, result_from_fit: robjects.ListVector, gev_param_name_to_dim=None,
burn_in_percentage=0.5) -> None:
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
new file mode 100644
index 0000000000000000000000000000000000000000..079c4d63b529071b4ab22aa754c7f8f18af48f2f
--- /dev/null
+++ b/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_mle_extremes.py
@@ -0,0 +1,15 @@
+import numpy as np
+
+from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_result_from_extremes import \
+ AbstractResultFromExtremes
+from extreme_fit.model.result_from_model_fit.utils import get_margin_coef_ordered_dict
+
+
+class ResultFromMleExtremes(AbstractResultFromExtremes):
+
+ @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)
diff --git a/test/test_extreme_fit/test_distribution/test_gev/test_gev_temporal_bayesian.py b/test/test_extreme_fit/test_distribution/test_gev/test_gev_temporal_extremes_bayesian.py
similarity index 85%
rename from test/test_extreme_fit/test_distribution/test_gev/test_gev_temporal_bayesian.py
rename to test/test_extreme_fit/test_distribution/test_gev/test_gev_temporal_extremes_bayesian.py
index 7a620674f0040c6901165787a712fd0885eda2e5..868521460f7a94c890694fbc4f9b5b8e77a892cc 100644
--- a/test/test_extreme_fit/test_distribution/test_gev/test_gev_temporal_bayesian.py
+++ b/test/test_extreme_fit/test_distribution/test_gev/test_gev_temporal_extremes_bayesian.py
@@ -4,13 +4,10 @@ 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.distribution.gev.gev_params import GevParams
-from extreme_fit.estimator.margin_estimator.abstract_margin_estimator import LinearMarginEstimator
from extreme_fit.model.margin_model.linear_margin_model.abstract_temporal_linear_margin_model import \
- AbstractTemporalLinearMarginModel, TemporalMarginFitMethod
+ TemporalMarginFitMethod
from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import StationaryTemporalModel, \
NonStationaryLocationTemporalModel, NonStationaryLocationAndScaleTemporalModel
-from extreme_fit.model.result_from_model_fit.result_from_extremes import ResultFromExtremes
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 \
@@ -18,10 +15,9 @@ from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_
from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
from spatio_temporal_dataset.spatio_temporal_observations.abstract_spatio_temporal_observations import \
AbstractSpatioTemporalObservations
-from test.test_utils import load_non_stationary_temporal_margin_models
-class TestGevTemporalBayesian(unittest.TestCase):
+class TestGevTemporalExtremesBayesian(unittest.TestCase):
def setUp(self) -> None:
set_seed_r()
@@ -42,12 +38,12 @@ class TestGevTemporalBayesian(unittest.TestCase):
def test_gev_temporal_margin_fit_stationary(self):
# Create estimator
- estimator_fitted = fitted_linear_margin_estimator(StationaryTemporalModel, self.coordinates, self.dataset,
+ estimator = fitted_linear_margin_estimator(StationaryTemporalModel, self.coordinates, self.dataset,
starting_year=0,
fit_method=self.fit_method)
ref = {'loc': 0.34272436381693616, 'scale': 1.3222588712831973, 'shape': 0.30491484962825105}
for year in range(1, 3):
- mle_params_estimated = estimator_fitted.margin_function_from_fit.get_gev_params(np.array([year])).to_dict()
+ 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)
diff --git a/test/test_extreme_fit/test_distribution/test_gev/test_gev_temporal_extremes_mle.py b/test/test_extreme_fit/test_distribution/test_gev/test_gev_temporal_extremes_mle.py
new file mode 100644
index 0000000000000000000000000000000000000000..0d420386f879e805c32d2de566462b366df78c3b
--- /dev/null
+++ b/test/test_extreme_fit/test_distribution/test_gev/test_gev_temporal_extremes_mle.py
@@ -0,0 +1,73 @@
+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
+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 TestGevTemporalExtremesMle(unittest.TestCase):
+
+ def setUp(self) -> None:
+ set_seed_r()
+ r("""
+ N <- 50
+ loc = 0; scale = 1; shape <- 1
+ 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)
+ self.fit_method = TemporalMarginFitMethod.extremes_fevd_mle
+
+ def test_gev_temporal_margin_fit_stationary(self):
+ # Create estimator
+ estimator = fitted_linear_margin_estimator(StationaryTemporalModel, self.coordinates, self.dataset,
+ starting_year=0,
+ fit_method=self.fit_method)
+ ref = {'loc': 0.02191974259369493, 'scale': 1.0347946062900268, 'shape': 0.829052520147379}
+ 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)
+
+ def test_gev_temporal_margin_fit_non_stationary_location(self):
+ # Create estimator
+ estimator = fitted_linear_margin_estimator(NonStationaryLocationTemporalModel, self.coordinates, self.dataset,
+ starting_year=0,
+ fit_method=self.fit_method)
+ # Checks that parameters returned are indeed different
+ mle_params_estimated_year1 = estimator.margin_function_from_fit.get_gev_params(np.array([1])).to_dict()
+ mle_params_estimated_year3 = estimator.margin_function_from_fit.get_gev_params(np.array([3])).to_dict()
+ self.assertNotEqual(mle_params_estimated_year1, mle_params_estimated_year3)
+
+ def test_gev_temporal_margin_fit_non_stationary_location_and_scale(self):
+ # Create estimator
+ estimator = fitted_linear_margin_estimator(NonStationaryLocationAndScaleTemporalModel, self.coordinates,
+ self.dataset,
+ starting_year=0,
+ fit_method=self.fit_method)
+ # Checks that parameters returned are indeed different
+ mle_params_estimated_year1 = estimator.margin_function_from_fit.get_gev_params(np.array([1])).to_dict()
+ mle_params_estimated_year3 = estimator.margin_function_from_fit.get_gev_params(np.array([3])).to_dict()
+ self.assertNotEqual(mle_params_estimated_year1, mle_params_estimated_year3)
+
+
+if __name__ == '__main__':
+ unittest.main()
diff --git a/test/test_extreme_fit/test_distribution/test_gev/test_gev_temporal.py b/test/test_extreme_fit/test_distribution/test_gev/test_gev_temporal_is_mev.py
similarity index 87%
rename from test/test_extreme_fit/test_distribution/test_gev/test_gev_temporal.py
rename to test/test_extreme_fit/test_distribution/test_gev/test_gev_temporal_is_mev.py
index 116c041e60357d80307e66787879c6e06e228080..2e52087c5324cd49b23f4ab264dbb1b20a324fec 100644
--- a/test/test_extreme_fit/test_distribution/test_gev/test_gev_temporal.py
+++ b/test/test_extreme_fit/test_distribution/test_gev/test_gev_temporal_is_mev.py
@@ -3,7 +3,10 @@ import unittest
import numpy as np
import pandas as pd
+from experiment.trend_analysis.univariate_test.utils import fitted_linear_margin_estimator
from extreme_fit.estimator.margin_estimator.abstract_margin_estimator import LinearMarginEstimator
+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
from extreme_fit.model.utils import r, set_seed_r
@@ -33,12 +36,13 @@ class TestGevTemporal(unittest.TestCase):
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)
+ self.fit_method = TemporalMarginFitMethod.is_mev_gev_fit
def test_gev_temporal_margin_fit_stationary(self):
# Create estimator
- margin_model = StationaryTemporalModel(self.coordinates)
- estimator = LinearMarginEstimator(self.dataset, margin_model)
- estimator.fit()
+ estimator = fitted_linear_margin_estimator(StationaryTemporalModel, self.coordinates, self.dataset,
+ starting_year=0,
+ fit_method=self.fit_method)
ref = {'loc': 0.04309190816463247, 'scale': 2.0688696961628437, 'shape': 0.8291528207825063}
for year in range(1, 3):
mle_params_estimated = estimator.margin_function_from_fit.get_gev_params(np.array([year])).to_dict()
@@ -49,7 +53,6 @@ class TestGevTemporal(unittest.TestCase):
# Create estimator
margin_models = load_non_stationary_temporal_margin_models(self.coordinates)
for margin_model in margin_models:
- # margin_model = NonStationaryLocationStationModel(self.coordinates)
estimator = LinearMarginEstimator(self.dataset, margin_model)
estimator.fit()
# Checks that parameters returned are indeed different
@@ -71,7 +74,8 @@ class TestGevTemporal(unittest.TestCase):
self.assertNotEqual(mle_params_estimated_year5, mle_params_estimated_year3)
def fit_non_stationary_estimator(self, starting_point):
- margin_model = NonStationaryLocationTemporalModel(self.coordinates, starting_point=starting_point + self.start_year)
+ margin_model = NonStationaryLocationTemporalModel(self.coordinates,
+ starting_point=starting_point + self.start_year)
estimator = LinearMarginEstimator(self.dataset, margin_model)
estimator.fit()
return estimator