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 6ddb75c8f1bc10d693e6c697f510d46a974ac6ba..aa9b8afb76d6f61c76623a4041a7a546b7b750a6 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 8b62f3571e421b64a210477babe25f8cf5885d55..4f5cf39bed754d32e861d6b34963b448a6c9a5ef 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 3d31dd08d44fad9204bb11cfefa3f957721bd3cb..545e9784238862baf7548443d1bdfaa9c5508c31 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 8cf8990083759d9c1bae433cdc0308f0da44b64d..e24978c3413a43424e80a5f7c7e6cdd0b58af55b 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 dbdf9fe0085e2858b69a2477988fa9c29101a487..8a51ceee46a53bbe774db444b38bbf7adb680864 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 2fcf945c86f49b9b5fd92c49ccef652cfefbe01b..c73cc0890444be47a6b1168ade69d7cf09579f70 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 81818b43c67466fabe9c67da25344a3bc4ea7768..429dcab67e7344ba33548ef32473dbd101f7ee77 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 e29f75912eb9eeccaced9862357367b737055b59..7959897dfa44377665f3bd91723a81e93c4b8cee 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__':