[contrasting project] add l_moments method for estimating stationary mle

extreme_fit/model/margin_model/linear_margin_model/abstract_temporal_linear_margin_model.py

@@ -21,6 +21,7 @@ class TemporalMarginFitMethod(Enum):
     extremes_fevd_bayesian = 1
     extremes_fevd_mle = 2
     extremes_fevd_gmle = 3
+    extremes_fevd_l_moments = 4
 
 
 class AbstractTemporalLinearMarginModel(LinearMarginModel):
@@ -53,7 +54,9 @@ class AbstractTemporalLinearMarginModel(LinearMarginModel):
             elif self.fit_method == TemporalMarginFitMethod.extremes_fevd_bayesian:
                 return self.extremes_fevd_bayesian_fit(x, df_coordinates_temp)
             elif self.fit_method in [TemporalMarginFitMethod.extremes_fevd_mle,
-                                     TemporalMarginFitMethod.extremes_fevd_gmle]:
+                                     TemporalMarginFitMethod.extremes_fevd_gmle,
+                                     TemporalMarginFitMethod.extremes_fevd_l_moments,
+                                     ]:
                 return self.extremes_fevd_mle_related_fit(x, df_coordinates_temp)
             else:
                 raise NotImplementedError
@@ -78,6 +81,9 @@ class AbstractTemporalLinearMarginModel(LinearMarginModel):
             method = "MLE"
         elif self.fit_method == TemporalMarginFitMethod.extremes_fevd_gmle:
             method = "GMLE"
+        elif self.fit_method == TemporalMarginFitMethod.extremes_fevd_l_moments:
+            method = "Lmoments"
+            assert self.margin_function_start_fit.is_a_stationary_model
         else:
             raise ValueError('wrong method')
         return self.run_fevd_fixed(df_coordinates_temp, method, x)

extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_mle_extremes.py

@@ -21,7 +21,10 @@ class ResultFromMleExtremes(AbstractResultFromExtremes):
     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']))}
+        if 'par' in d:
+            values = {i: param for i, param in enumerate(np.array(d['par']))}
+        else:
+            values = {i: np.array(v)[0] for i, v in enumerate(d.values())}
         return get_margin_coef_ordered_dict(self.gev_param_name_to_dim, values, self.type_for_mle)
 
     def _confidence_interval_method(self, common_kwargs, ci_method, return_period):

test/test_extreme_fit/test_distribution/test_gev/test_gev_temporal_extremes_l_moments.py

+import unittest
+
+import numpy as np
+import pandas as pd
+
+from extreme_trend.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 TestGevTemporalExtremesLMoments(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_l_moments
+
+    ref = {'loc': 0.0813843045950251, 'scale': 1.1791830110181365, 'shape': 0.6610403806908737}
+
+    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.0813843045950251, 'scale': 1.1791830110181365, 'shape': 0.6610403806908737}
+        for year in range(1, 3):
+            mle_params_estimated = estimator.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
+        with self.assertRaises(AssertionError):
+            fitted_linear_margin_estimator(NonStationaryLocationTemporalModel, self.coordinates,
+                                                       self.dataset,
+                                                       starting_year=0,
+                                                       fit_method=self.fit_method)
+
+
+if __name__ == '__main__':
+    unittest.main()