From 2bc5363c89dbf5df59e20f6c4ab296440094199a Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Fri, 20 Mar 2020 16:29:57 +0100
Subject: [PATCH] [quantile regression project] first stationary run for the
gev simulation
---
.../AbstractSimulation.py | 28 +++++-----
.../GevSimulation.py | 51 ++++++++++++++++---
.../test_gev_simulations.py | 8 ++-
3 files changed, 62 insertions(+), 25 deletions(-)
diff --git a/projects/quantile_regression_vs_evt/AbstractSimulation.py b/projects/quantile_regression_vs_evt/AbstractSimulation.py
index 3a013a04..14dbb926 100644
--- a/projects/quantile_regression_vs_evt/AbstractSimulation.py
+++ b/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()
diff --git a/projects/quantile_regression_vs_evt/GevSimulation.py b/projects/quantile_regression_vs_evt/GevSimulation.py
index 317c8e83..9ce61302 100644
--- a/projects/quantile_regression_vs_evt/GevSimulation.py
+++ b/projects/quantile_regression_vs_evt/GevSimulation.py
@@ -1,5 +1,13 @@
-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)
diff --git a/test/test_projects/test_quantile_regression/test_gev_simulations.py b/test/test_projects/test_quantile_regression/test_gev_simulations.py
index acddcdea..455bd518 100644
--- a/test/test_projects/test_quantile_regression/test_gev_simulations.py
+++ b/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)
--
GitLab