[paper 1] add nllh test to check computation

extreme_fit/estimator/margin_estimator/abstract_margin_estimator.py

@@ -3,11 +3,12 @@ from abc import ABC
 from cached_property import cached_property
 
 from extreme_fit.estimator.abstract_estimator import AbstractEstimator
-from extreme_fit.estimator.utils import load_margin_function
+from extreme_fit.estimator.utils import load_margin_function, compute_nllh
 from extreme_fit.model.margin_model.linear_margin_model.linear_margin_model import LinearMarginModel
 from extreme_fit.model.margin_model.margin_function.linear_margin_function import LinearMarginFunction
 from extreme_fit.model.result_from_model_fit.abstract_result_from_model_fit import AbstractResultFromModelFit
 from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
+from spatio_temporal_dataset.slicer.split import Split
 
 
 class AbstractMarginEstimator(AbstractEstimator, ABC):
@@ -26,13 +27,24 @@ class LinearMarginEstimator(AbstractMarginEstimator):
         self.margin_model = margin_model
 
     def _fit(self) -> AbstractResultFromModelFit:
-        maxima_gev_specialized = self.dataset.maxima_gev_for_spatial_extremes_package(self.train_split)
         df_coordinates_spat = self.dataset.coordinates.df_spatial_coordinates(self.train_split)
-        df_coordinates_temp = self.dataset.coordinates.df_temporal_coordinates_for_fit(split=self.train_split,
-                                                                                       starting_point=self.margin_model.starting_point)
-        return self.margin_model.fitmargin_from_maxima_gev(data=maxima_gev_specialized,
+        return self.margin_model.fitmargin_from_maxima_gev(data=self.maxima_gev_train,
                                                            df_coordinates_spat=df_coordinates_spat,
-                                                           df_coordinates_temp=df_coordinates_temp)
+                                                           df_coordinates_temp=self.coordinate_temp)
+
+    @property
+    def coordinate_temp(self):
+        return self.dataset.coordinates.df_temporal_coordinates_for_fit(split=self.train_split,
+                                                                        starting_point=self.margin_model.starting_point)
+
+    @property
+    def maxima_gev_train(self):
+        return self.dataset.maxima_gev_for_spatial_extremes_package(self.train_split)
+
+    @property
+    def nllh(self, split=Split.all):
+        assert split is Split.all
+        return compute_nllh(self, self.maxima_gev_train, self.coordinate_temp, self.margin_model)
 
     @cached_property
     def margin_function_from_fit(self) -> LinearMarginFunction:

extreme_fit/estimator/utils.py

+import numpy as np
+
 from extreme_fit.estimator.abstract_estimator import AbstractEstimator
 from extreme_fit.model.margin_model.linear_margin_model.linear_margin_model import LinearMarginModel
 from extreme_fit.model.margin_model.margin_function.linear_margin_function import LinearMarginFunction
 
 
-def load_margin_function(estimator: AbstractEstimator, margin_model: LinearMarginModel, margin_function_class=LinearMarginFunction, coef_dict=None):
+def load_margin_function(estimator: AbstractEstimator, margin_model: LinearMarginModel,
+                         margin_function_class=LinearMarginFunction, coef_dict=None):
     if coef_dict is None:
         coef_dict = estimator.result_from_model_fit.margin_coef_ordered_dict
     return margin_function_class.from_coef_dict(coordinates=estimator.dataset.coordinates,
                                                 gev_param_name_to_dims=margin_model.margin_function_start_fit.gev_param_name_to_dims,
                                                 coef_dict=coef_dict,
                                                 starting_point=margin_model.starting_point)
+
+
+def compute_nllh(estimator: AbstractEstimator, maxima, coordinate_temp, margin_model: LinearMarginModel,
+                 margin_function_class=LinearMarginFunction, coef_dict=None):
+    margin_function = load_margin_function(estimator, margin_model, margin_function_class, coef_dict)
+    nllh = 0
+    for maximum, year in zip(maxima[0], coordinate_temp.values):
+        gev_params = margin_function.get_gev_params(year, is_transformed=False)
+        p = gev_params.density(maximum)
+        nllh -= np.log(p)
+    return nllh

test/test_extreme_fit/test_distribution/test_gev/test_gev_temporal_extremes_mle.py

@@ -46,6 +46,7 @@ class TestGevTemporalExtremesMle(unittest.TestCase):
             mle_params_estimated = estimator.margin_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)
+            self.assertAlmostEqual(estimator.result_from_model_fit.nllh, estimator.nllh)
 
     def test_gev_temporal_margin_fit_non_stationary_location(self):
         # Create estimator
@@ -56,6 +57,7 @@ class TestGevTemporalExtremesMle(unittest.TestCase):
         mle_params_estimated_year1 = estimator.margin_function_from_fit.get_gev_params(np.array([1])).to_dict()
         mle_params_estimated_year3 = estimator.margin_function_from_fit.get_gev_params(np.array([3])).to_dict()
         self.assertNotEqual(mle_params_estimated_year1, mle_params_estimated_year3)
+        self.assertAlmostEqual(estimator.result_from_model_fit.nllh, estimator.nllh)
 
     def test_gev_temporal_margin_fit_non_stationary_location_and_scale(self):
         # Create estimator
@@ -67,6 +69,7 @@ class TestGevTemporalExtremesMle(unittest.TestCase):
         mle_params_estimated_year1 = estimator.margin_function_from_fit.get_gev_params(np.array([1])).to_dict()
         mle_params_estimated_year3 = estimator.margin_function_from_fit.get_gev_params(np.array([3])).to_dict()
         self.assertNotEqual(mle_params_estimated_year1, mle_params_estimated_year3)
+        self.assertAlmostEqual(estimator.result_from_model_fit.nllh, estimator.nllh)
 
 
 if __name__ == '__main__':