From b5f6f87811d3905be008f3a5c45e02e1ff0ad393 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Fri, 2 Oct 2020 13:27:52 +0200
Subject: [PATCH] [contrasting] add general test for
test_gev_spatio_temporal_polynomial_extremes_mle.py. refactor a bit
---
.../abstract_extract_eurocode_return_level.py | 4 +-
.../abstract_result_from_extremes.py | 4 +-
.../eurocode_return_level_uncertainties.py | 10 ++--
.../result_from_mle_extremes.py | 2 -
.../model/result_from_model_fit/utils.py | 1 -
.../one_fold_analysis/one_fold_fit.py | 17 ++++---
.../main_bayesian_mcmc.py | 2 +-
...spatio_temporal_polynomial_extremes_mle.py | 46 ++++++-------------
8 files changed, 34 insertions(+), 52 deletions(-)
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 6ddb75c8..aa9b8afb 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
@@ -69,8 +69,8 @@ class ExtractEurocodeReturnLevelFromCiMethod(AbstractExtractEurocodeReturnLevel)
class ExtractEurocodeReturnLevelFromMyBayesianExtremes(AbstractExtractEurocodeReturnLevel):
result_from_fit: ResultFromBayesianExtremes
- def __init__(self, estimator: LinearMarginEstimator, ci_method, temporal_covariate, quantile_level=EUROCODE_QUANTILE):
- super().__init__(estimator, ci_method, temporal_covariate, quantile_level)
+ def __init__(self, estimator: LinearMarginEstimator, ci_method, coordinate, quantile_level=EUROCODE_QUANTILE):
+ super().__init__(estimator, ci_method, coordinate, quantile_level)
assert isinstance(self.result_from_fit, ResultFromBayesianExtremes)
@property
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 8b62f357..4f5cf39b 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
@@ -61,15 +61,13 @@ class AbstractResultFromExtremes(AbstractResultFromModelFit):
'type': r.c("return.level")
}
if self.param_name_to_dim:
- print("here")
if isinstance(transformed_temporal_covariate, (int, float, np.int, np.float, np.int64)):
- print("here2")
d = {GevParams.greek_letter_from_param_name_confidence_interval(param_name) + '1': r.c(transformed_temporal_covariate) for
param_name in self.param_name_to_dim.keys()}
elif isinstance(transformed_temporal_covariate, np.ndarray):
d = OrderedDict()
linearity_in_shape = GevParams.SHAPE in self.param_name_to_dim
- nb_calls = 2 # or 4 (1 and 3 did not work for the test)
+ nb_calls = 4 # or 4 (1 and 3 did not work for the test)
for param_name in GevParams.PARAM_NAMES:
suffix = '0' if param_name in self.param_name_to_dim else ''
covariate = np.array([1] * nb_calls)
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 3d31dd08..545e9784 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,3 +1,7 @@
+from typing import Union
+
+import numpy as np
+
from extreme_data.eurocode_data.utils import EUROCODE_QUANTILE
from extreme_fit.estimator.margin_estimator.utils import fitted_linear_margin_estimator
from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_extract_eurocode_return_level import \
@@ -29,12 +33,12 @@ class EurocodeConfidenceIntervalFromExtremes(object):
@classmethod
def from_estimator_extremes(cls, estimator_extremes: LinearMarginEstimator,
ci_method: ConfidenceIntervalMethodFromExtremes,
- temporal_covariate: int,
+ coordinate: Union[int, np.ndarray],
):
if ci_method == ConfidenceIntervalMethodFromExtremes.my_bayes:
- extractor = ExtractEurocodeReturnLevelFromMyBayesianExtremes(estimator_extremes, ci_method, temporal_covariate, cls.quantile_level)
+ extractor = ExtractEurocodeReturnLevelFromMyBayesianExtremes(estimator_extremes, ci_method, coordinate, cls.quantile_level)
else:
- extractor = ExtractEurocodeReturnLevelFromCiMethod(estimator_extremes, ci_method, temporal_covariate, cls.quantile_level)
+ extractor = ExtractEurocodeReturnLevelFromCiMethod(estimator_extremes, ci_method, coordinate, cls.quantile_level)
return cls(extractor.mean_estimate, extractor.confidence_interval)
@classmethod
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 8cf89900..e24978c3 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
@@ -23,7 +23,6 @@ class ResultFromMleExtremes(AbstractResultFromExtremes):
d = self.get_python_dictionary(values)
if 'par' in d:
values = {i: param for i, param in enumerate(np.array(d['par']))}
- print(values)
else:
values = {i: np.array(v)[0] for i, v in enumerate(d.values())}
return get_margin_coef_ordered_dict(self.param_name_to_dim, values, self.type_for_mle,
@@ -38,7 +37,6 @@ class ResultFromMleExtremes(AbstractResultFromExtremes):
res = r('ci.fevd.mle_fixed')(self.result_from_fit, **mle_ci_parameters, **common_kwargs)
if self.is_non_stationary:
b = np.array(res)
- print(b)
a = b[0]
lower, mean_estimate, upper, _ = a
else:
diff --git a/extreme_fit/model/result_from_model_fit/utils.py b/extreme_fit/model/result_from_model_fit/utils.py
index dbdf9fe0..8a51ceee 100644
--- a/extreme_fit/model/result_from_model_fit/utils.py
+++ b/extreme_fit/model/result_from_model_fit/utils.py
@@ -41,5 +41,4 @@ def get_margin_coef_ordered_dict(param_name_to_dims, mle_values, type_for_mle="G
coef_name = coef_template.format(k + offset)
coef_dict[coef_name] = mle_values[i]
i += 1
- print(coef_dict)
return coef_dict
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 2fcf945c..c73cc089 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
@@ -273,11 +273,10 @@ class OneFoldFit(object):
has_not_opposite_sign = self.sign_of_change(estimator) * self.sign_of_change(new_estimator) >= 0
return has_not_opposite_sign
-
def sensitivity_of_fit_test_last_years(self, estimator: LinearMarginEstimator):
# Build the dataset without the maxima for each altitude
new_dataset = AbstractDataset.remove_last_maxima(self.dataset.observations,
- self.dataset.coordinates,
+ self.dataset.coordinates,
nb_years=10)
# Fit the new estimator
model_class = type(estimator.margin_model)
@@ -288,7 +287,7 @@ class OneFoldFit(object):
# Compare sign of change
has_not_opposite_sign = self.sign_of_change(estimator) * self.sign_of_change(new_estimator) >= 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))
+ # 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):
@@ -312,9 +311,9 @@ class OneFoldFit(object):
standard_gumbel_quantiles = [standard_gumbel_distribution.quantile(i / (n + 1)) for i in range(1, n + 1)]
return standard_gumbel_quantiles
-
- # def best_confidence_interval(self):
- # EurocodeConfidenceIntervalFromExtremes.from_estimator_extremes(self.best_estimator,
- # ci_method=ConfidenceIntervalMethodFromExtremes.ci_mle,
- # temporal_covariate=np.array([2019, self.altitude_plot]),)
-
+ @property
+ def best_confidence_interval(self) -> EurocodeConfidenceIntervalFromExtremes:
+ return EurocodeConfidenceIntervalFromExtremes.from_estimator_extremes(self.best_estimator,
+ ci_method=ConfidenceIntervalMethodFromExtremes.ci_mle,
+ coordinate=np.array(
+ [2019, self.altitude_plot]), )
diff --git a/projects/others/check_mcmc_convergence_for_return_levels/main_bayesian_mcmc.py b/projects/others/check_mcmc_convergence_for_return_levels/main_bayesian_mcmc.py
index 81818b43..429dcab6 100644
--- a/projects/others/check_mcmc_convergence_for_return_levels/main_bayesian_mcmc.py
+++ b/projects/others/check_mcmc_convergence_for_return_levels/main_bayesian_mcmc.py
@@ -107,7 +107,7 @@ def get_return_level_bayesian_example(nb_iterations_for_bayesian_fit):
nb_iterations_for_bayesian_fit=nb_iterations_for_bayesian_fit)
return_level_bayesian = ExtractEurocodeReturnLevelFromMyBayesianExtremes(estimator=fitted_estimator,
ci_method=ConfidenceIntervalMethodFromExtremes.my_bayes,
- temporal_covariate=2019)
+ coordinate=2019)
return return_level_bayesian
diff --git a/test/test_extreme_fit/test_estimator/test_spatio_temporal_estimator/test_gev_spatio_temporal_polynomial_extremes_mle.py b/test/test_extreme_fit/test_estimator/test_spatio_temporal_estimator/test_gev_spatio_temporal_polynomial_extremes_mle.py
index e29f7591..7959897d 100644
--- a/test/test_extreme_fit/test_estimator/test_spatio_temporal_estimator/test_gev_spatio_temporal_polynomial_extremes_mle.py
+++ b/test/test_extreme_fit/test_estimator/test_spatio_temporal_estimator/test_gev_spatio_temporal_polynomial_extremes_mle.py
@@ -14,6 +14,8 @@ from extreme_fit.model.margin_model.polynomial_margin_model.models_based_on_pari
from extreme_fit.model.margin_model.polynomial_margin_model.polynomial_margin_model import \
NonStationaryQuadraticLocationModel, \
NonStationaryQuadraticScaleModel, NonStationaryQuadraticLocationGumbelModel, NonStationaryQuadraticScaleGumbelModel
+from extreme_fit.model.margin_model.polynomial_margin_model.utils import \
+ ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS
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 \
@@ -76,45 +78,27 @@ class TestGevTemporalQuadraticExtremesMle(unittest.TestCase):
self.assertAlmostEqual(estimator.result_from_model_fit.aic, estimator.aic())
self.assertAlmostEqual(estimator.result_from_model_fit.bic, estimator.bic())
# Assert we can compute the return level
- covariate1_for_return_level = np.array([700, 0])
- covariate2_for_return_level = np.array([700, 1000])
+ covariate1_for_return_level = np.array([500, 0])
+ covariate2_for_return_level = np.array([500, 50])
covariate3_for_return_level = np.array([400, 0])
coordinates = [covariate1_for_return_level, covariate2_for_return_level, covariate3_for_return_level]
for coordinate in coordinates:
EurocodeConfidenceIntervalFromExtremes.quantile_level = 0.98
confidence_interval = EurocodeConfidenceIntervalFromExtremes.from_estimator_extremes(estimator,
ci_method=ConfidenceIntervalMethodFromExtremes.ci_mle,
- temporal_covariate=coordinate)
+ coordinate=coordinate)
gev_params = estimator.function_from_fit.get_params(coordinate)
- # print(gev_params, coordinate)
return_level = gev_params.return_level(return_period=50)
- # print("my return level", return_level)
- # print("their return level", confidence_interval.mean_estimate)
- # print("diff", confidence_interval.mean_estimate - return_level)
- # print('\n\n')
- self.assertAlmostEqual(return_level, confidence_interval.mean_estimate)
- self.assertFalse(np.isnan(confidence_interval.confidence_interval[0]))
- self.assertFalse(np.isnan(confidence_interval.confidence_interval[1]))
-
- def test_gev_spatio_temporal_margin_fit_1(self):
- self.function_test_gev_spatio_temporal_margin_fit_non_stationary(StationaryAltitudinal)
-
- #
- def test_gev_spatio_temporal_margin_fit_1_bis(self):
- self.function_test_gev_spatio_temporal_margin_fit_non_stationary(AltitudinalShapeLinearTimeStationary)
-
- def test_gev_spatio_temporal_margin_fit_2(self):
- # first model with both a time and altitude non stationarity
- self.function_test_gev_spatio_temporal_margin_fit_non_stationary(
- AltitudinalShapeConstantTimeLocationLinear)
-
- # def test_gev_spatio_temporal_margin_fit_3(self):
- # self.function_test_gev_spatio_temporal_margin_fit_non_stationary(
- # AltitudinalShapeLinearTimeLocationLinear)
- #
- # def test_gev_spatio_temporal_margin_fit_4(self):
- # self.function_test_gev_spatio_temporal_margin_fit_non_stationary(
- # AltitudinalShapeLinearTimeLocScaleLinear)
+ if np.isnan(return_level) or np.isnan(confidence_interval.mean_estimate):
+ self.assertTrue(np.isnan(return_level) and np.isnan(confidence_interval.mean_estimate))
+ else:
+ self.assertAlmostEqual(return_level, confidence_interval.mean_estimate)
+ self.assertFalse(np.isnan(confidence_interval.confidence_interval[0]))
+ self.assertFalse(np.isnan(confidence_interval.confidence_interval[1]))
+
+ def test_gev_spatio_temporal_all(self):
+ for model_class in ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS:
+ self.function_test_gev_spatio_temporal_margin_fit_non_stationary(model_class)
if __name__ == '__main__':
--
GitLab