[projection snowfall] add test_gev_temporal_polynomial_evgam.py and pass the tests.

extreme_fit/model/margin_model/linear_margin_model/abstract_temporal_linear_margin_model.py

@@ -76,6 +76,8 @@ class AbstractTemporalLinearMarginModel(LinearMarginModel):
         df = pd.DataFrame({maxima_column_name: np.array(x)})
         df = pd.concat([df, df_coordinates_spat, df_coordinates_temp], axis=1)
         data = get_r_dataframe_from_python_dataframe(df)
+        if self.type_for_mle is not "GEV":
+            raise NotImplementedError
         res = safe_run_r_estimator(function=r('evgam'),
                                    formula=formula,
                                    data=data,

extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_evgam.py

 import numpy as np
 from rpy2 import robjects
+from rpy2.robjects.pandas2ri import ri2py_dataframe
 
+from extreme_fit.distribution.gev.gev_params import GevParams
 from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_result_from_extremes import \
     AbstractResultFromExtremes
 from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import \
@@ -18,16 +20,39 @@ class ResultFromEvgam(AbstractResultFromExtremes):
         self.type_for_mle = type_for_mle
 
     @property
-    def aic(self):
-        """Compute the aic from the list of parameters in the results,
+    def nllh(self):
+        """Compute the nllh from the list of parameters in the results,
          find a way to comptue it directly from the result to compare"""
-        location = self.name_to_value['location']
-        a = np.array(location)
-        print(a)
-        print(len(a))
-        print('here2')
-        # 'location', 'logscale', 'shape'
-        raise NotImplementedError
+        param_names = ['location', 'logscale', 'shape']
+        parameters = [np.array(self.get_python_dictionary(self.name_to_value[param_name])['fitted'])
+                      for param_name in param_names]
+        # Add maxima
+        parameters.append(self.maxima)
+        for p in parameters:
+            assert len(p) == len(self.maxima)
+        # Compute nllh
+        nllh = 0
+        for loc, log_scale, shape, maximum in zip(*parameters):
+            gev_params = GevParams(loc, np.exp(log_scale), shape)
+            p = gev_params.density(maximum)
+            nllh += -np.log(p)
+        return nllh
+
+    @property
+    def maxima(self):
+        return np.array(self.get_python_dictionary(self.name_to_value['location'])['model'][0])
+
+    @property
+    def nb_parameters(self):
+        return len(np.array(self.name_to_value['coefficients']))
+
+    @property
+    def aic(self):
+        return 2 * self.nllh + 2 * self.nb_parameters
+
+    @property
+    def bic(self):
+        return 2 * self.nllh + np.log(len(self.maxima)) * self.nb_parameters
 
     @property
     def log_scale(self):

test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_polynomial_evgam.py

+import unittest
+
+import numpy as np
+import pandas as pd
+
+from extreme_fit.distribution.gev.gev_params import GevParams
+from extreme_fit.model.margin_model.polynomial_margin_model.polynomial_margin_model import NonStationaryQuadraticLocationModel, \
+    NonStationaryQuadraticScaleModel, NonStationaryQuadraticLocationGumbelModel, NonStationaryQuadraticScaleGumbelModel
+from extreme_trend.trend_test.abstract_gev_trend_test import fitted_linear_margin_estimator
+from extreme_fit.model.margin_model.utils import \
+    MarginFitMethod
+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 TestGevTemporalPolynomialEvgam(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 = MarginFitMethod.evgam
+
+    def function_test_gev_temporal_margin_fit_non_stationary_quadratic(self, model_class, quadratic_param):
+        # Create estimator
+        estimator = fitted_linear_margin_estimator(model_class,
+                                                   self.coordinates, self.dataset,
+                                                   starting_year=0,
+                                                   fit_method=self.fit_method)
+        # Checks that parameters returned are indeed different
+        mle_params_estimated_year1 = estimator.function_from_fit.get_params(np.array([1])).to_dict()
+        mle_params_estimated_year3 = estimator.function_from_fit.get_params(np.array([21])).to_dict()
+        mle_params_estimated_year5 = estimator.function_from_fit.get_params(np.array([41])).to_dict()
+        self.assertNotEqual(mle_params_estimated_year1, mle_params_estimated_year3)
+        self.assertNotEqual(mle_params_estimated_year3, mle_params_estimated_year5)
+        # Assert the relationship for the location is indeed quadratic
+        diff1 = mle_params_estimated_year1[quadratic_param] - mle_params_estimated_year3[quadratic_param]
+        diff2 = mle_params_estimated_year3[quadratic_param] - mle_params_estimated_year5[quadratic_param]
+        self.assertNotAlmostEqual(diff1, diff2)
+        # Assert that indicators are correctly computed
+        self.assertAlmostEqual(estimator.result_from_model_fit.nllh, estimator.nllh())
+        self.assertAlmostEqual(estimator.result_from_model_fit.aic, estimator.aic())
+        self.assertAlmostEqual(estimator.result_from_model_fit.bic, estimator.bic())
+
+    def test_gev_temporal_margin_fit_non_stationary_quadratic_location(self):
+        self.function_test_gev_temporal_margin_fit_non_stationary_quadratic(NonStationaryQuadraticLocationModel,
+                                                                            quadratic_param=GevParams.LOC)
+
+    def test_gev_temporal_margin_fit_non_stationary_quadratic_scale(self):
+        self.function_test_gev_temporal_margin_fit_non_stationary_quadratic(NonStationaryQuadraticScaleModel,
+                                                                            quadratic_param=GevParams.SCALE)
+
+
+if __name__ == '__main__':
+    unittest.main()