From 65d8de60e3471194953fb3b9f0e24d54f6498cbb Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Mon, 30 Mar 2020 19:16:57 +0200
Subject: [PATCH] [contrasting project] add fitted_stationary_gev. remove gev
fit classes
---
.../visualization/study_visualizer.py | 6 +--
extreme_fit/distribution/gev/gev_fit.py | 15 ------
extreme_fit/distribution/gev/gevmle_fit.py | 29 -----------
extreme_fit/distribution/gev/ismev_gev_fit.py | 49 ------------------
.../estimator/margin_estimator/utils.py | 25 +++++++++
.../margin_function/linear_margin_function.py | 4 ++
.../presentation/statistical_model.py | 6 +--
.../test_gev/test_gevmle_fit.py | 44 ----------------
.../test_gev/test_stationary_gev_fit.py | 51 +++++++++++++++++++
9 files changed, 85 insertions(+), 144 deletions(-)
delete mode 100644 extreme_fit/distribution/gev/gev_fit.py
delete mode 100644 extreme_fit/distribution/gev/gevmle_fit.py
delete mode 100644 extreme_fit/distribution/gev/ismev_gev_fit.py
delete mode 100644 test/test_extreme_fit/test_distribution/test_gev/test_gevmle_fit.py
create mode 100644 test/test_extreme_fit/test_distribution/test_gev/test_stationary_gev_fit.py
diff --git a/extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py b/extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py
index a75e43b3..8a4fa435 100644
--- a/extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py
+++ b/extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py
@@ -12,6 +12,7 @@ import numpy as np
import pandas as pd
import seaborn as sns
+from extreme_fit.estimator.margin_estimator.utils import fitted_stationary_gev
from extreme_fit.model.result_from_model_fit.result_from_extremes.eurocode_return_level_uncertainties import \
EurocodeConfidenceIntervalFromExtremes, compute_eurocode_confidence_interval
from extreme_data.meteo_france_data.scm_models_data.abstract_extended_study import AbstractExtendedStudy
@@ -29,7 +30,6 @@ from extreme_fit.function.margin_function.abstract_margin_function import \
from extreme_fit.function.param_function.param_function import AbstractParamFunction
from extreme_fit.model.max_stable_model.abstract_max_stable_model import CovarianceFunction
from extreme_fit.distribution.gev.gev_params import GevParams
-from extreme_fit.distribution.gev.ismev_gev_fit import IsmevGevFit
from extreme_fit.distribution.gpd.gpd_params import GpdParams
from extreme_fit.distribution.gpd.gpdmle_fit import GpdMleFit
from spatio_temporal_dataset.coordinates.spatial_coordinates.abstract_spatial_coordinates import \
@@ -374,7 +374,7 @@ class StudyVisualizer(VisualizationParameters):
# Display the graph of the max on top
ax = plt.gca()
_, y = self.smooth_maxima_x_y(massif_names.index(massif_name))
- gev_param = IsmevGevFit(x_gev=y).gev_params
+ gev_param = fitted_stationary_gev(x_gev=y)
# Round up
# d = {k: self.round_sig(v, 2) for k, v in d.items()}
@@ -661,7 +661,7 @@ class StudyVisualizer(VisualizationParameters):
@cached_property
def massif_name_to_gev_mle_fitted(self) -> Dict[str, GevParams]:
- return {massif_name: IsmevGevFit(self.df_maxima_gev.loc[massif_name]).gev_params
+ return {massif_name: fitted_stationary_gev(self.df_maxima_gev.loc[massif_name])
for massif_name in self.study.study_massif_names}
@cached_property
diff --git a/extreme_fit/distribution/gev/gev_fit.py b/extreme_fit/distribution/gev/gev_fit.py
deleted file mode 100644
index ef81d978..00000000
--- a/extreme_fit/distribution/gev/gev_fit.py
+++ /dev/null
@@ -1,15 +0,0 @@
-import numpy as np
-
-from extreme_fit.distribution.gev.gev_params import GevParams
-
-
-class GevFit(object):
-
- def __init__(self, x_gev: np.ndarray, block_size=None):
- assert np.ndim(x_gev) == 1
- assert len(x_gev) > 1
- self.x_gev = x_gev
-
- @property
- def gev_params(self) -> GevParams:
- raise NotImplementedError
\ No newline at end of file
diff --git a/extreme_fit/distribution/gev/gevmle_fit.py b/extreme_fit/distribution/gev/gevmle_fit.py
deleted file mode 100644
index 41288004..00000000
--- a/extreme_fit/distribution/gev/gevmle_fit.py
+++ /dev/null
@@ -1,29 +0,0 @@
-import numpy as np
-
-from extreme_fit.distribution.gev.gev_fit import GevFit
-from extreme_fit.distribution.gev.gev_params import GevParams
-from extreme_fit.model.utils import r
-
-"""
-These two functions are “extremely light†functions to fit the GEV/GPD. These functions are mainlyuseful
-to compute starting values for the Schlather and Smith mode
-If more refined (univariate) analysis have to be performed, users should use more specialised pack-ages
- - e.g. POT, evd, ismev, . . . .
-"""
-
-
-def spatial_extreme_gevmle_fit(x_gev):
- res = r.gevmle(x_gev, method="Nelder")
- return dict(zip(res.names, res))
-
-
-class GevMleFit(GevFit):
-
- def __init__(self, x_gev: np.ndarray, block_size=None):
- super().__init__(x_gev, block_size)
- self._gev_params = spatial_extreme_gevmle_fit(x_gev)
- self.gev_params_object = GevParams.from_dict({**self._gev_params, 'block_size': block_size})
-
- @property
- def gev_params(self) -> GevParams:
- return self.gev_params_object
diff --git a/extreme_fit/distribution/gev/ismev_gev_fit.py b/extreme_fit/distribution/gev/ismev_gev_fit.py
deleted file mode 100644
index ca7997aa..00000000
--- a/extreme_fit/distribution/gev/ismev_gev_fit.py
+++ /dev/null
@@ -1,49 +0,0 @@
-import numpy as np
-
-from extreme_fit.distribution.gev.gev_fit import GevFit
-from extreme_fit.distribution.gev.gev_params import GevParams
-import rpy2.robjects as ro
-
-from extreme_fit.model.utils import r, get_null
-
-
-"""
-IsMev fit functions
-"""
-
-
-def ismev_gev_fit(x_gev, y, mul):
- """
- For non-stationary fitting it is recommended that the covariates within the generalized linear models are
- (at least approximately) centered and scaled (i.e.the columns of ydat should be approximately centered and scaled).
- """
- xdat = ro.FloatVector(x_gev)
- gev_fit = r('gev.fit')
- y = y if y is not None else get_null()
- mul = mul if mul is not None else get_null()
- res = gev_fit(xdat, y, mul)
-
- # r.assign('python_wrapping', True)
- # r.source(file="""/home/erwan/Documents/projects/spatiotemporalextremes/extreme_fit/distribution/gev/wrapper_ismev_gev_fit.R""")
- # y = np.arange(1, len(x_gev), 1).reshape((-11, 1))
- # res = r.gev_fit(data=xdat, y=y, trend=1)
- return dict(zip(res.names, res))
-
-class IsmevGevFit(GevFit):
- # todo: this could be modeled with a margin_function depending only on time
- # todo: I should remove the call to this object, in order to centralize all the calls
-
- def __init__(self, x_gev: np.ndarray, y=None, mul=None):
- super().__init__(x_gev)
- self.y = y
- self.mul = mul
- self.res = ismev_gev_fit(x_gev, y, mul)
-
- @property
- def gev_params(self) -> GevParams:
- assert self.y is None
- gev_params_dict = dict(zip(GevParams.PARAM_NAMES, self.res['mle']))
- return GevParams.from_dict(gev_params_dict)
-
-
-
diff --git a/extreme_fit/estimator/margin_estimator/utils.py b/extreme_fit/estimator/margin_estimator/utils.py
index 74e872a7..95dc74d0 100644
--- a/extreme_fit/estimator/margin_estimator/utils.py
+++ b/extreme_fit/estimator/margin_estimator/utils.py
@@ -1,4 +1,16 @@
+import pandas as pd
+
+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.temporal_linear_margin_models import StationaryTemporalModel
+from spatio_temporal_dataset.coordinates.temporal_coordinates.generated_temporal_coordinates import \
+ ConsecutiveTemporalCoordinates
+from spatio_temporal_dataset.coordinates.transformed_coordinates.transformation.abstract_transformation import \
+ CenteredScaledNormalization
+from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
+from spatio_temporal_dataset.spatio_temporal_observations.abstract_spatio_temporal_observations import \
+ AbstractSpatioTemporalObservations
+from spatio_temporal_dataset.spatio_temporal_observations.annual_maxima_observations import AnnualMaxima
def fitted_linear_margin_estimator(model_class, coordinates, dataset, starting_year, fit_method, **model_kwargs):
@@ -6,3 +18,16 @@ def fitted_linear_margin_estimator(model_class, coordinates, dataset, starting_y
estimator = LinearMarginEstimator(dataset, model)
estimator.fit()
return estimator
+
+
+def fitted_stationary_gev(x_gev, fit_method, model_class=StationaryTemporalModel, starting_year=None,
+ transformation_class=CenteredScaledNormalization) -> GevParams:
+ coordinates = ConsecutiveTemporalCoordinates.from_nb_temporal_steps(nb_temporal_steps=len(x_gev),
+ transformation_class=CenteredScaledNormalization)
+ df = pd.DataFrame([pd.Series(dict(zip(coordinates.index, x_gev)))]).transpose()
+ observations = AnnualMaxima(df_maxima_gev=df)
+ dataset = AbstractDataset(observations=observations, coordinates=coordinates)
+ estimator = fitted_linear_margin_estimator(model_class, coordinates, dataset, starting_year, fit_method)
+ first_coordinate = coordinates.coordinates_values()[0]
+ gev_param = estimator.function_from_fit.get_gev_params(first_coordinate)
+ return gev_param
diff --git a/extreme_fit/function/margin_function/linear_margin_function.py b/extreme_fit/function/margin_function/linear_margin_function.py
index 7deb1ab3..81a9ee8d 100644
--- a/extreme_fit/function/margin_function/linear_margin_function.py
+++ b/extreme_fit/function/margin_function/linear_margin_function.py
@@ -61,6 +61,10 @@ class LinearMarginFunction(ParametricMarginFunction):
coef_dict.update(coef.coef_dict(dims, self.idx_to_coefficient_name(self.coordinates)))
return coef_dict
+ @property
+ def is_a_stationary_model(self) -> bool:
+ return all([v == 'NULL' for v in self.form_dict.values()])
+
@property
def form_dict(self) -> Dict[str, str]:
form_dict = {}
diff --git a/projects/exceeding_snow_loads/presentation/statistical_model.py b/projects/exceeding_snow_loads/presentation/statistical_model.py
index 969eb23e..1a5e8dcd 100644
--- a/projects/exceeding_snow_loads/presentation/statistical_model.py
+++ b/projects/exceeding_snow_loads/presentation/statistical_model.py
@@ -1,8 +1,7 @@
import matplotlib.pyplot as plt
import numpy as np
-from extreme_fit.distribution.gev.gevmle_fit import GevMleFit
-from extreme_fit.distribution.gev.ismev_gev_fit import ismev_gev_fit
+from extreme_fit.estimator.margin_estimator.utils import fitted_stationary_gev
def binomial_observation():
@@ -34,8 +33,7 @@ def histogram_for_gev():
study = CrocusSnowLoadTotal(altitude=1800)
s = study.observations_annual_maxima.df_maxima_gev.loc['Vercors']
x_gev = s.values
- gev = GevMleFit(x_gev)
- gev_params = gev.gev_params
+ gev_params = fitted_stationary_gev(x_gev)
samples = gev_params.sample(10000)
nb = 10
epsilon = 0.0
diff --git a/test/test_extreme_fit/test_distribution/test_gev/test_gevmle_fit.py b/test/test_extreme_fit/test_distribution/test_gev/test_gevmle_fit.py
deleted file mode 100644
index 184be222..00000000
--- a/test/test_extreme_fit/test_distribution/test_gev/test_gevmle_fit.py
+++ /dev/null
@@ -1,44 +0,0 @@
-import unittest
-
-import numpy as np
-
-from extreme_fit.model.utils import r, set_seed_r
-from extreme_fit.distribution.gev.gev_params import GevParams
-from extreme_fit.distribution.gev.gevmle_fit import GevMleFit
-from extreme_fit.distribution.gev.ismev_gev_fit import IsmevGevFit
-
-
-class TestGevMleFit(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
- """)
-
- def test_gevmle_fit(self):
- estimator = GevMleFit(x_gev=np.array(r['x_gev']))
- mle_params_ref = {'loc': 0.0219, 'scale': 1.0347, 'shape': 0.8290}
- self.fit_estimator(estimator, ref=mle_params_ref)
-
- def test_ismev_gev_fit(self):
- estimator = IsmevGevFit(x_gev=np.array(r['x_gev']))
- ismev_ref = {'loc': 0.0219, 'scale': 1.0347, 'shape': 0.8295}
- self.fit_estimator(estimator, ismev_ref)
-
- # def test_s
-
- def fit_estimator(self, estimator, ref):
- # Compare the MLE estimated parameters to the reference
- mle_params_estimated = estimator.gev_params
- self.assertIsInstance(mle_params_estimated, GevParams)
- mle_params_estimated = mle_params_estimated.to_dict()
- for key in ref.keys():
- self.assertAlmostEqual(ref[key], mle_params_estimated[key], places=3)
-
-
-if __name__ == '__main__':
- unittest.main()
diff --git a/test/test_extreme_fit/test_distribution/test_gev/test_stationary_gev_fit.py b/test/test_extreme_fit/test_distribution/test_gev/test_stationary_gev_fit.py
new file mode 100644
index 00000000..737bab80
--- /dev/null
+++ b/test/test_extreme_fit/test_distribution/test_gev/test_stationary_gev_fit.py
@@ -0,0 +1,51 @@
+import unittest
+
+import numpy as np
+
+from extreme_fit.estimator.margin_estimator.utils import fitted_stationary_gev
+from extreme_fit.model.margin_model.linear_margin_model.abstract_temporal_linear_margin_model import \
+ TemporalMarginFitMethod
+from extreme_fit.model.utils import r, set_seed_r
+from extreme_fit.distribution.gev.gev_params import GevParams
+
+
+class TestStationaryGevFit(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
+ """)
+
+ def test_stationary_gev_fit_with_ismev(self):
+ params_estimated = fitted_stationary_gev(x_gev=np.array(r['x_gev']),
+ fit_method=TemporalMarginFitMethod.is_mev_gev_fit)
+ ref = {'loc': 0.0219, 'scale': 1.0347, 'shape': 0.8295}
+ self.fit_estimator(params_estimated, ref)
+
+ def test_stationary_gev_fit_with_mle(self):
+ params_estimated = fitted_stationary_gev(x_gev=np.array(r['x_gev']),
+ fit_method=TemporalMarginFitMethod.extremes_fevd_mle)
+ ref = {'loc': 0.02191974259369493, 'scale': 1.0347946062900268, 'shape': 0.829052520147379}
+ self.fit_estimator(params_estimated, ref)
+
+ def test_stationary_gev_fit_with_l_moments(self):
+ params_estimated = fitted_stationary_gev(x_gev=np.array(r['x_gev']),
+ fit_method=TemporalMarginFitMethod.extremes_fevd_l_moments)
+ ref = {'loc': 0.0813843045950251, 'scale': 1.1791830110181365, 'shape': 0.6610403806908737}
+ self.fit_estimator(params_estimated, ref)
+
+ def fit_estimator(self, params_estimated, ref):
+ # Compare the MLE estimated parameters to the reference
+ self.assertIsInstance(params_estimated, GevParams)
+ params_estimated = params_estimated.to_dict()
+ print(params_estimated)
+ for key in ref.keys():
+ self.assertAlmostEqual(ref[key], params_estimated[key], places=3)
+
+
+if __name__ == '__main__':
+ unittest.main()
--
GitLab