[contrasting project] add fitted_stationary_gev. remove gev fit classes

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

extreme_fit/distribution/gev/gev_fit.py

-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

extreme_fit/distribution/gev/gevmle_fit.py

-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

extreme_fit/distribution/gev/ismev_gev_fit.py

-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)
-
-
-

extreme_fit/estimator/margin_estimator/utils.py

+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

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 = {}

projects/exceeding_snow_loads/presentation/statistical_model.py

 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

test/test_extreme_fit/test_distribution/test_gev/test_gevmle_fit.py

-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()

test/test_extreme_fit/test_distribution/test_gev/test_stationary_gev_fit.py

+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()