[quantile regression project] first stationary run for the gev simulation

projects/quantile_regression_vs_evt/AbstractSimulation.py

@@ -5,11 +5,14 @@ from collections import OrderedDict
 import numpy as np
 from cached_property import cached_property
 
+from extreme_fit.estimator.quantile_estimator.abstract_quantile_estimator import AbstractQuantileEstimator
 from extreme_fit.estimator.quantile_estimator.quantile_estimator_from_margin import QuantileEstimatorFromMargin
 from extreme_fit.estimator.quantile_estimator.quantile_estimator_from_regression import QuantileRegressionEstimator
 from extreme_fit.model.margin_model.linear_margin_model.abstract_temporal_linear_margin_model import \
     AbstractTemporalLinearMarginModel
 from extreme_fit.model.quantile_model.quantile_regression_model import AbstractQuantileRegressionModel
+from root_utils import get_display_name_from_object_type
+from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
 from spatio_temporal_dataset.coordinates.temporal_coordinates.generated_temporal_coordinates import \
     ConsecutiveTemporalCoordinates
 from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
@@ -31,22 +34,18 @@ class AbstractSimulation(object):
         self.time_series_lengths = time_series_lengths
         self.nb_time_series = nb_time_series
 
-    def generate_all_observation(self, nb_time_series, length, coordinates) -> List[AbstractSpatioTemporalObservations]:
+    def generate_all_observation(self, nb_time_series, length) -> List[AbstractSpatioTemporalObservations]:
         raise NotImplementedError
 
     @cached_property
     def time_serie_length_to_observation_list(self) -> Dict[int, List[AbstractSpatioTemporalObservations]]:
         d = OrderedDict()
         for length in self.time_series_lengths:
-            if self.multiprocessing:
-                raise NotImplementedError
-            else:
-                coordinates = self.time_serie_length_to_coordinates[length]
-                d[length] = self.generate_all_observation(self.nb_time_series, length, coordinates)
+            d[length] = self.generate_all_observation(self.nb_time_series, length)
         return d
 
     @cached_property
-    def time_serie_length_to_coordinates(self) -> Dict[int, Coordinates]:
+    def time_serie_length_to_coordinates(self) -> Dict[int, AbstractCoordinates]:
         d = OrderedDict()
         for length in self.time_series_lengths:
             d[length] = ConsecutiveTemporalCoordinates.from_nb_temporal_steps(length)
@@ -83,22 +82,25 @@ class AbstractSimulation(object):
         for model_class, d_sub in self.model_class_to_time_serie_length_to_estimator_fitted.items():
             length_to_error_values = OrderedDict()
             for length, estimators_fitted in d_sub.items():
-                errors = []
-                # Add the mean, and the quantile of the 95% confidence interval
-                error_values = [0.2, 1, 1.3]
+                errors = self.compute_errors(length, estimators_fitted)
+                error_values = [np.quantile(errors, q=0.025), np.mean(errors), np.quantile(errors, q=0.975)]
                 length_to_error_values[length] = error_values
             d[model_class] = length_to_error_values
-        print(d)
         return d
 
+    def compute_errors(self, length: int, estimators_fitted: List[AbstractQuantileEstimator]):
+        raise NotImplementedError
+
     def plot_error_for_last_year_quantile(self, show=True):
         ax = plt.gca()
         for model_class, length_to_error_values in self.model_class_to_error_last_year_quantile.items():
             lengths = list(length_to_error_values.keys())
             errors_values = np.array(list(length_to_error_values.values()))
-            print(errors_values.shape)
             mean_error = errors_values[:, 1]
-            ax.plot(lengths, mean_error, label=str(model_class))
+            label = get_display_name_from_object_type(model_class)
+            ax.plot(lengths, mean_error, label=label)
+            ax.set_xlabel('# Data')
+            ax.set_ylabel('Absolute error for the {} quantile at the last coordinate'.format(self.quantile))
             ax.legend()
         if show:
             plt.show()

projects/quantile_regression_vs_evt/GevSimulation.py

-from typing import List
+from collections import OrderedDict
+from typing import List, Dict
 
+import numpy as np
+from cached_property import cached_property
+
+from extreme_fit.distribution.gev.gev_params import GevParams
+from extreme_fit.estimator.quantile_estimator.abstract_quantile_estimator import AbstractQuantileEstimator
+from extreme_fit.model.margin_model.abstract_margin_model import AbstractMarginModel
+from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import StationaryTemporalModel
 from projects.quantile_regression_vs_evt.AbstractSimulation import AbstractSimulation
 from spatio_temporal_dataset.spatio_temporal_observations.abstract_spatio_temporal_observations import \
     AbstractSpatioTemporalObservations
@@ -8,11 +16,40 @@ from spatio_temporal_dataset.spatio_temporal_observations.annual_maxima_observat
 
 class GevSimulation(AbstractSimulation):
 
-    def __init__(self, margin_model_class_for_simulation, nb_time_series, quantile, time_series_lengths=None, multiprocessing=False, model_classes=None):
-        super().__init__(nb_time_series, quantile, time_series_lengths, multiprocessing, model_classes)
-        self.margin_model_class_for_simulation = margin_model_class_for_simulation
+    @cached_property
+    def time_series_lengths_to_margin_model(self) -> Dict[int, AbstractMarginModel]:
+        d = OrderedDict()
+        for length in self.time_series_lengths:
+            coordinates = self.time_serie_length_to_coordinates[length]
+            d[length] = self.create_model(coordinates)
+        return d
+
+    def create_model(self, coordinates):
+        raise NotImplementedError
+
+    def generate_all_observation(self, nb_time_series, length) -> List[AbstractSpatioTemporalObservations]:
+        coordinates = self.time_serie_length_to_coordinates[length]
+        margin_model = self.time_series_lengths_to_margin_model[length]
+        return [MarginAnnualMaxima.from_sampling(nb_obs=length, coordinates=coordinates, margin_model=margin_model)
+                for _ in range(nb_time_series)]
+
+    def compute_errors(self, length: int, estimators: List[AbstractQuantileEstimator]):
+        coordinates = self.time_serie_length_to_coordinates[length]
+        last_coordinate = coordinates.coordinates_values()[-1]
+        # Compute true value
+        margin_model = self.time_series_lengths_to_margin_model[length]
+        true_quantile = margin_model.margin_function_sample.get_gev_params(last_coordinate).quantile(self.quantile)
+        # Compute estimated values
+        estimated_quantiles = [estimator.function_from_fit.get_quantile(last_coordinate) for estimator in estimators]
+        return np.abs(np.array(estimated_quantiles) - true_quantile)
+
 
-    def generate_all_observation(self, nb_time_series, length, coordinates) -> List[AbstractSpatioTemporalObservations]:
-        margin_model = self.margin_model_class_for_simulation(coordinates)
-        return [MarginAnnualMaxima.from_sampling(nb_obs=length, coordinates=coordinates, margin_model=margin_model) for _ in range(nb_time_series)]
+class StationarySimulation(GevSimulation):
 
+    def create_model(self, coordinates):
+        gev_param_name_to_coef_list = {
+            GevParams.LOC: [0],
+            GevParams.SHAPE: [0],
+            GevParams.SCALE: [1],
+        }
+        return StationaryTemporalModel.from_coef_list(coordinates, gev_param_name_to_coef_list)

test/test_projects/test_quantile_regression/test_gev_simulations.py

@@ -2,17 +2,15 @@ import unittest
 
 from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import StationaryTemporalModel
 from extreme_fit.model.quantile_model.quantile_regression_model import ConstantQuantileRegressionModel
-from projects.quantile_regression_vs_evt.GevSimulation import GevSimulation
+from projects.quantile_regression_vs_evt.GevSimulation import GevSimulation, StationarySimulation
 
 
 class TestGevSimulations(unittest.TestCase):
     DISPLAY = False
 
     def test_stationary_run(self):
-        margin_model_class_for_simulation = StationaryTemporalModel
-        simulation = GevSimulation(margin_model_class_for_simulation,
-                                   nb_time_series=1, quantile=0.5, time_series_lengths=[50, 60],
-                                   model_classes=[StationaryTemporalModel, ConstantQuantileRegressionModel])
+        simulation = StationarySimulation(nb_time_series=1, quantile=0.5, time_series_lengths=[50, 60],
+                                          model_classes=[StationaryTemporalModel, ConstantQuantileRegressionModel])
         simulation.plot_error_for_last_year_quantile(self.DISPLAY)