From 9c88536f3488bc2c7f07bffb66659f7892447cc3 Mon Sep 17 00:00:00 2001 From: Le Roux Erwan <erwan.le-roux@irstea.fr> Date: Mon, 23 Mar 2020 21:47:48 +0100 Subject: [PATCH] [quantile regression project] add annual maxima osbervation from daily exponential. add sammple_r_function as argument to rmargin_from_nb_obs. add abstract extreme params --- .../distribution/abstract_extreme_params.py | 20 +++++++++++++ extreme_fit/distribution/abstract_params.py | 13 ++------- extreme_fit/distribution/exp_params.py | 17 ++++++++--- extreme_fit/distribution/gev/gev_params.py | 3 +- extreme_fit/distribution/gpd/gpd_params.py | 3 +- .../abstract_margin_function.py | 3 +- .../combined_margin_function.py | 2 +- .../independent_margin_function.py | 14 ++++----- .../margin_function/linear_margin_function.py | 6 ++-- .../parametric_margin_function.py | 7 +++-- .../margin_model/abstract_margin_model.py | 8 +++-- .../abstract_temporal_linear_margin_model.py | 6 ++-- .../linear_margin_model.py | 29 ++++++++++--------- .../margin_model/parametric_margin_model.py | 5 ++-- .../daily_observations.py | 3 +- .../test_spatio_temporal_observations.py | 15 +++++++--- 16 files changed, 98 insertions(+), 56 deletions(-) create mode 100644 extreme_fit/distribution/abstract_extreme_params.py diff --git a/extreme_fit/distribution/abstract_extreme_params.py b/extreme_fit/distribution/abstract_extreme_params.py new file mode 100644 index 00000000..53cd1147 --- /dev/null +++ b/extreme_fit/distribution/abstract_extreme_params.py @@ -0,0 +1,20 @@ +from extreme_fit.distribution.abstract_params import AbstractParams + + +class AbstractExtremeParams(AbstractParams): + pass + + # Extreme parameters + SCALE = 'scale' + LOC = 'loc' + SHAPE = 'shape' + + def __init__(self, loc: float, scale: float, shape: float): + self.location = loc + self.scale = scale + self.shape = shape + # By default, scale cannot be negative + # (sometimes it happens, when we want to find a quantile for every point of a 2D map + # then it can happen that a corner point that was not used for fitting correspond to a negative scale, + # in the case we set all the parameters as equal to np.nan, and we will not display those points) + self.has_undefined_parameters = self.scale <= 0 \ No newline at end of file diff --git a/extreme_fit/distribution/abstract_params.py b/extreme_fit/distribution/abstract_params.py index f367a820..631f8656 100644 --- a/extreme_fit/distribution/abstract_params.py +++ b/extreme_fit/distribution/abstract_params.py @@ -17,21 +17,14 @@ class AbstractParams(object): QUANTILE_COLORS = ['orange', 'red', 'darkviolet'] # Summary SUMMARY_NAMES = PARAM_NAMES + QUANTILE_NAMES + + # Simple parameters + RATE = 'rate' # Extreme parameters SCALE = 'scale' LOC = 'loc' SHAPE = 'shape' - def __init__(self, loc: float, scale: float, shape: float): - self.location = loc - self.scale = scale - self.shape = shape - # By default, scale cannot be negative - # (sometimes it happens, when we want to find a quantile for every point of a 2D map - # then it can happen that a corner point that was not used for fitting correspond to a negative scale, - # in the case we set all the parameters as equal to np.nan, and we will not display those points) - self.has_undefined_parameters = self.scale <= 0 - @classmethod def from_dict(cls, params: dict): return cls(**params) diff --git a/extreme_fit/distribution/exp_params.py b/extreme_fit/distribution/exp_params.py index 6a003ee2..5ffb82f9 100644 --- a/extreme_fit/distribution/exp_params.py +++ b/extreme_fit/distribution/exp_params.py @@ -1,9 +1,18 @@ +import numpy as np + from extreme_fit.distribution.abstract_params import AbstractParams -class ExpParams(object): - PARAM_NAMES = [AbstractParams.SCALE] +class ExpParams(AbstractParams): + PARAM_NAMES = [AbstractParams.RATE] - def __init__(self, scale: float): - self.scale = scale + def __init__(self, rate) -> None: + self.rate = rate + self.has_undefined_parameters = self.rate < 0 + @property + def param_values(self): + if self.has_undefined_parameters: + return [np.nan for _ in range(1)] + else: + return [self.rate] diff --git a/extreme_fit/distribution/gev/gev_params.py b/extreme_fit/distribution/gev/gev_params.py index db84b878..a116f245 100644 --- a/extreme_fit/distribution/gev/gev_params.py +++ b/extreme_fit/distribution/gev/gev_params.py @@ -5,13 +5,14 @@ from typing import List from cached_property import cached_property from mpmath import euler, pi +from extreme_fit.distribution.abstract_extreme_params import AbstractExtremeParams from extreme_fit.distribution.abstract_params import AbstractParams from extreme_fit.model.utils import r import numpy as np from scipy.special import gamma -class GevParams(AbstractParams): +class GevParams(AbstractExtremeParams): # Parameters PARAM_NAMES = [AbstractParams.LOC, AbstractParams.SCALE, AbstractParams.SHAPE] # Summary diff --git a/extreme_fit/distribution/gpd/gpd_params.py b/extreme_fit/distribution/gpd/gpd_params.py index c440415f..ed6b1ac9 100644 --- a/extreme_fit/distribution/gpd/gpd_params.py +++ b/extreme_fit/distribution/gpd/gpd_params.py @@ -1,8 +1,9 @@ +from extreme_fit.distribution.abstract_extreme_params import AbstractExtremeParams from extreme_fit.distribution.abstract_params import AbstractParams from extreme_fit.model.utils import r -class GpdParams(AbstractParams): +class GpdParams(AbstractExtremeParams): # TODO: understand better why the gpdfit return 2 parameters, alors que d'un autre cote d autres definitions de la distribution parlent d un parametre location # Parameters diff --git a/extreme_fit/function/margin_function/abstract_margin_function.py b/extreme_fit/function/margin_function/abstract_margin_function.py index df30892b..860d25c5 100644 --- a/extreme_fit/function/margin_function/abstract_margin_function.py +++ b/extreme_fit/function/margin_function/abstract_margin_function.py @@ -20,8 +20,9 @@ class AbstractMarginFunction(AbstractFunction): VISUALIZATION_RESOLUTION = 100 VISUALIZATION_TEMPORAL_STEPS = 2 - def __init__(self, coordinates: AbstractCoordinates): + def __init__(self, coordinates: AbstractCoordinates, params_class: type = GevParams): super().__init__(coordinates) + self.params_class = params_class self.mask_2D = None # Visualization parameters diff --git a/extreme_fit/function/margin_function/combined_margin_function.py b/extreme_fit/function/margin_function/combined_margin_function.py index 08737bc1..20b5e474 100644 --- a/extreme_fit/function/margin_function/combined_margin_function.py +++ b/extreme_fit/function/margin_function/combined_margin_function.py @@ -18,7 +18,7 @@ class CombinedMarginFunction(AbstractMarginFunction): def get_gev_params(self, coordinate: np.ndarray) -> GevParams: gev_params_list = [margin_function.get_gev_params(coordinate) for margin_function in self.margin_functions] mean_gev_params = np.mean(np.array([gev_param.to_array() for gev_param in gev_params_list]), axis=0) - gev_param = GevParams(*mean_gev_params) + gev_param = self.params_class(*mean_gev_params) return gev_param @classmethod diff --git a/extreme_fit/function/margin_function/independent_margin_function.py b/extreme_fit/function/margin_function/independent_margin_function.py index ccb6dba8..43169d71 100644 --- a/extreme_fit/function/margin_function/independent_margin_function.py +++ b/extreme_fit/function/margin_function/independent_margin_function.py @@ -14,9 +14,9 @@ class IndependentMarginFunction(AbstractMarginFunction): IndependentMarginFunction: each parameter of the GEV are modeled independently """ - def __init__(self, coordinates: AbstractCoordinates): + def __init__(self, coordinates: AbstractCoordinates, params_class: type = GevParams): """Attribute 'gev_param_name_to_param_function' maps each GEV parameter to its corresponding function""" - super().__init__(coordinates) + super().__init__(coordinates, params_class) self.gev_param_name_to_param_function = None # type: Union[None, Dict[str, AbstractParamFunction]] def get_gev_params(self, coordinate: np.ndarray, is_transformed: bool = True) -> GevParams: @@ -24,13 +24,11 @@ class IndependentMarginFunction(AbstractMarginFunction): # Since all the coordinates are usually transformed by default # then we assume that the input coordinate are transformed by default assert self.gev_param_name_to_param_function is not None - assert len(self.gev_param_name_to_param_function) == 3 + assert len(self.gev_param_name_to_param_function) == len(self.params_class.PARAM_NAMES) transformed_coordinate = coordinate if is_transformed else self.transform(coordinate) - gev_params = {} - for gev_param_name in GevParams.PARAM_NAMES: - param_function = self.gev_param_name_to_param_function[gev_param_name] - gev_params[gev_param_name] = param_function.get_param_value(transformed_coordinate) - return GevParams.from_dict(gev_params) + params = {param_name: param_function.get_param_value(transformed_coordinate) + for param_name, param_function in self.gev_param_name_to_param_function.items()} + return self.params_class.from_dict(params) diff --git a/extreme_fit/function/margin_function/linear_margin_function.py b/extreme_fit/function/margin_function/linear_margin_function.py index d89a7a59..3cdc1f84 100644 --- a/extreme_fit/function/margin_function/linear_margin_function.py +++ b/extreme_fit/function/margin_function/linear_margin_function.py @@ -28,9 +28,11 @@ class LinearMarginFunction(ParametricMarginFunction): COEF_CLASS = LinearCoef def __init__(self, coordinates: AbstractCoordinates, gev_param_name_to_dims: Dict[str, List[int]], - gev_param_name_to_coef: Dict[str, AbstractCoef], starting_point: Union[None, int] = None): + gev_param_name_to_coef: Dict[str, AbstractCoef], starting_point: Union[None, int] = None, + params_class: type = GevParams): self.gev_param_name_to_coef = None # type: Union[None, Dict[str, LinearCoef]] - super().__init__(coordinates, gev_param_name_to_dims, gev_param_name_to_coef, starting_point) + super().__init__(coordinates, gev_param_name_to_dims, gev_param_name_to_coef, starting_point, + params_class) def load_specific_param_function(self, gev_param_name) -> AbstractParamFunction: return LinearParamFunction(dims=self.gev_param_name_to_dims[gev_param_name], diff --git a/extreme_fit/function/margin_function/parametric_margin_function.py b/extreme_fit/function/margin_function/parametric_margin_function.py index 57037216..60a1bbb8 100644 --- a/extreme_fit/function/margin_function/parametric_margin_function.py +++ b/extreme_fit/function/margin_function/parametric_margin_function.py @@ -31,10 +31,11 @@ class ParametricMarginFunction(IndependentMarginFunction): COEF_CLASS = None def __init__(self, coordinates: AbstractCoordinates, gev_param_name_to_dims: Dict[str, List[int]], - gev_param_name_to_coef: Dict[str, AbstractCoef], starting_point: Union[None, int] = None): + gev_param_name_to_coef: Dict[str, AbstractCoef], starting_point: Union[None, int] = None, + params_class: type = GevParams): # Starting point for the trend is the same for all the parameters self.starting_point = starting_point - super().__init__(coordinates) + super().__init__(coordinates, params_class) self.gev_param_name_to_dims = gev_param_name_to_dims # type: Dict[str, List[int]] # Check the dimension are well-defined with respect to the coordinates @@ -48,7 +49,7 @@ class ParametricMarginFunction(IndependentMarginFunction): # Build gev_parameter_to_param_function dictionary self.gev_param_name_to_param_function = {} # type: Dict[str, AbstractParamFunction] # Map each gev_param_name to its corresponding param_function - for gev_param_name in GevParams.PARAM_NAMES: + for gev_param_name in self.params_class.PARAM_NAMES: # By default, if dims are not specified, a constantParamFunction is chosen if self.gev_param_name_to_dims.get(gev_param_name) is None: param_function = ConstantParamFunction(constant=self.gev_param_name_to_coef[gev_param_name].intercept) diff --git a/extreme_fit/model/margin_model/abstract_margin_model.py b/extreme_fit/model/margin_model/abstract_margin_model.py index b4ac6918..79638873 100644 --- a/extreme_fit/model/margin_model/abstract_margin_model.py +++ b/extreme_fit/model/margin_model/abstract_margin_model.py @@ -20,12 +20,14 @@ class AbstractMarginModel(AbstractModel, ABC): """ def __init__(self, coordinates: AbstractCoordinates, use_start_value=False, - params_start_fit=None, params_sample=None): + params_start_fit=None, params_sample=None, + params_class=GevParams): super().__init__(use_start_value, params_start_fit, params_sample) assert isinstance(coordinates, AbstractCoordinates), type(coordinates) self.coordinates = coordinates self.margin_function_sample = None # type: AbstractMarginFunction self.margin_function_start_fit = None # type: AbstractMarginFunction + self.params_class = params_class self.load_margin_functions() def load_margin_functions(self): @@ -34,9 +36,9 @@ class AbstractMarginModel(AbstractModel, ABC): def default_load_margin_functions(self, margin_function_class): # todo: check it i could remove these attributes self.margin_function_sample = margin_function_class(coordinates=self.coordinates, - default_params=GevParams.from_dict(self.params_sample)) + default_params=self.params_class.from_dict(self.params_sample)) self.margin_function_start_fit = margin_function_class(coordinates=self.coordinates, - default_params=GevParams.from_dict( + default_params=self.params_class.from_dict( self.params_start_fit)) # Conversion class methods diff --git a/extreme_fit/model/margin_model/linear_margin_model/abstract_temporal_linear_margin_model.py b/extreme_fit/model/margin_model/linear_margin_model/abstract_temporal_linear_margin_model.py index e5e2ee45..69f6531e 100644 --- a/extreme_fit/model/margin_model/linear_margin_model/abstract_temporal_linear_margin_model.py +++ b/extreme_fit/model/margin_model/linear_margin_model/abstract_temporal_linear_margin_model.py @@ -29,8 +29,10 @@ class AbstractTemporalLinearMarginModel(LinearMarginModel): fit_method=TemporalMarginFitMethod.is_mev_gev_fit, nb_iterations_for_bayesian_fit=5000, params_start_fit_bayesian=None, - type_for_MLE="GEV"): - super().__init__(coordinates, use_start_value, params_start_fit, params_sample, starting_point) + type_for_MLE="GEV", + params_class=GevParams): + super().__init__(coordinates, use_start_value, params_start_fit, params_sample, starting_point, + params_class) self.type_for_mle = type_for_MLE self.params_start_fit_bayesian = params_start_fit_bayesian self.nb_iterations_for_bayesian_fit = nb_iterations_for_bayesian_fit diff --git a/extreme_fit/model/margin_model/linear_margin_model/linear_margin_model.py b/extreme_fit/model/margin_model/linear_margin_model/linear_margin_model.py index 08ac9900..422bb806 100644 --- a/extreme_fit/model/margin_model/linear_margin_model/linear_margin_model.py +++ b/extreme_fit/model/margin_model/linear_margin_model/linear_margin_model.py @@ -7,12 +7,12 @@ from extreme_fit.distribution.gev.gev_params import GevParams class LinearMarginModel(ParametricMarginModel): @classmethod - def from_coef_list(cls, coordinates, gev_param_name_to_coef_list, **kwargs): + def from_coef_list(cls, coordinates, gev_param_name_to_coef_list, params_class=GevParams, **kwargs): params = {} for param_name, coef_list in gev_param_name_to_coef_list.items(): for idx, coef in enumerate(coef_list, -1): params[(param_name, idx)] = coef - return cls(coordinates, params_sample=params, params_start_fit=params, **kwargs) + return cls(coordinates, params_sample=params, params_start_fit=params, params_class=params_class, **kwargs) def load_margin_functions(self, gev_param_name_to_dims=None): assert gev_param_name_to_dims is not None, 'LinearMarginModel cannot be used for sampling/fitting \n' \ @@ -27,34 +27,37 @@ class LinearMarginModel(ParametricMarginModel): self.margin_function_sample = LinearMarginFunction(coordinates=self.coordinates, gev_param_name_to_coef=coef_sample, gev_param_name_to_dims=gev_param_name_to_dims, - starting_point=self.starting_point) + starting_point=self.starting_point, + params_class=self.params_class) # Load start fit coef coef_start_fit = self.gev_param_name_to_linear_coef(param_name_and_dim_to_coef=self.params_start_fit) self.margin_function_start_fit = LinearMarginFunction(coordinates=self.coordinates, gev_param_name_to_coef=coef_start_fit, gev_param_name_to_dims=gev_param_name_to_dims, - starting_point=self.starting_point) + starting_point=self.starting_point, + params_class=self.params_class) @property def default_param_name_and_dim_to_coef(self) -> dict: default_intercept = 1 default_slope = 0.01 gev_param_name_and_dim_to_coef = {} - for gev_param_name in GevParams.PARAM_NAMES: - gev_param_name_and_dim_to_coef[(gev_param_name, -1)] = default_intercept + for param_name in self.params_class.PARAM_NAMES: + gev_param_name_and_dim_to_coef[(param_name, -1)] = default_intercept for dim in self.coordinates.coordinates_dims: - gev_param_name_and_dim_to_coef[(gev_param_name, dim)] = default_slope + gev_param_name_and_dim_to_coef[(param_name, dim)] = default_slope return gev_param_name_and_dim_to_coef def gev_param_name_to_linear_coef(self, param_name_and_dim_to_coef): - gev_param_name_to_linear_coef = {} - for gev_param_name in GevParams.PARAM_NAMES: - idx_to_coef = {idx: param_name_and_dim_to_coef[(gev_param_name, idx)] for idx in + param_name_to_linear_coef = {} + param_names = list(set([e[0] for e in param_name_and_dim_to_coef.keys()])) + for param_name in param_names: + idx_to_coef = {idx: param_name_and_dim_to_coef[(param_name, idx)] for idx in [-1] + self.coordinates.coordinates_dims} - linear_coef = LinearCoef(gev_param_name=gev_param_name, idx_to_coef=idx_to_coef) - gev_param_name_to_linear_coef[gev_param_name] = linear_coef - return gev_param_name_to_linear_coef + linear_coef = LinearCoef(gev_param_name=param_name, idx_to_coef=idx_to_coef) + param_name_to_linear_coef[param_name] = linear_coef + return param_name_to_linear_coef class ConstantMarginModel(LinearMarginModel): diff --git a/extreme_fit/model/margin_model/parametric_margin_model.py b/extreme_fit/model/margin_model/parametric_margin_model.py index 61279f4a..724366c4 100644 --- a/extreme_fit/model/margin_model/parametric_margin_model.py +++ b/extreme_fit/model/margin_model/parametric_margin_model.py @@ -3,6 +3,7 @@ from abc import ABC import numpy as np import pandas as pd +from extreme_fit.distribution.gev.gev_params import GevParams from extreme_fit.function.margin_function.parametric_margin_function import \ ParametricMarginFunction from extreme_fit.model.result_from_model_fit.result_from_spatial_extreme import ResultFromSpatialExtreme @@ -15,14 +16,14 @@ from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoo class ParametricMarginModel(AbstractMarginModel, ABC): def __init__(self, coordinates: AbstractCoordinates, use_start_value=False, params_start_fit=None, - params_sample=None, starting_point=None): + params_sample=None, starting_point=None, params_class=GevParams): """ :param starting_point: starting coordinate for the temporal trend """ self.starting_point = starting_point self.margin_function_sample = None # type: ParametricMarginFunction self.margin_function_start_fit = None # type: ParametricMarginFunction - super().__init__(coordinates, use_start_value, params_start_fit, params_sample) + super().__init__(coordinates, use_start_value, params_start_fit, params_sample, params_class) def fitmargin_from_maxima_gev(self, data: np.ndarray, df_coordinates_spat: pd.DataFrame, df_coordinates_temp: pd.DataFrame) -> ResultFromSpatialExtreme: diff --git a/spatio_temporal_dataset/spatio_temporal_observations/daily_observations.py b/spatio_temporal_dataset/spatio_temporal_observations/daily_observations.py index 9483f3b3..b0044f33 100644 --- a/spatio_temporal_dataset/spatio_temporal_observations/daily_observations.py +++ b/spatio_temporal_dataset/spatio_temporal_observations/daily_observations.py @@ -21,6 +21,7 @@ class DailyExp(AbstractSpatioTemporalObservations): def from_sampling(cls, nb_obs: int, coordinates: AbstractCoordinates, margin_model: AbstractMarginModel): exponential_values = margin_model.rmargin_from_nb_obs(nb_obs=nb_obs, - coordinates_values=coordinates.coordinates_values()) + coordinates_values=coordinates.coordinates_values(), + sample_r_function='rexp') df_exponential_values = pd.DataFrame(data=exponential_values, index=coordinates.index) return cls(df_maxima_gev=df_exponential_values) diff --git a/test/test_spatio_temporal_dataset/test_spatio_temporal_observations.py b/test/test_spatio_temporal_dataset/test_spatio_temporal_observations.py index 716a5856..4786e54c 100644 --- a/test/test_spatio_temporal_dataset/test_spatio_temporal_observations.py +++ b/test/test_spatio_temporal_dataset/test_spatio_temporal_observations.py @@ -4,6 +4,7 @@ import numpy as np import pandas as pd from extreme_fit.distribution.abstract_params import AbstractParams +from extreme_fit.distribution.exp_params import ExpParams from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import StationaryTemporalModel from extreme_fit.model.utils import set_seed_for_test from spatio_temporal_dataset.coordinates.temporal_coordinates.generated_temporal_coordinates import \ @@ -33,18 +34,24 @@ class TestDailyObservations(unittest.TestCase): set_seed_for_test(seed=42) self.coordinates = ConsecutiveTemporalCoordinates.from_nb_temporal_steps(nb_temporal_steps=10) gev_param_name_to_coef_list = { - AbstractParams.SCALE: [1], + AbstractParams.RATE: [1], } - self.margin_model = StationaryTemporalModel.from_coef_list(self.coordinates, gev_param_name_to_coef_list) + self.margin_model = StationaryTemporalModel.from_coef_list(self.coordinates, gev_param_name_to_coef_list, + params_class=ExpParams) + + def test_instance_exp_params(self): + last_coordinate = self.coordinates.coordinates_values()[-1] + params = self.margin_model.margin_function_sample.get_gev_params(last_coordinate) + self.assertIsInstance(params, ExpParams) def test_exponential_observations(self): obs = DailyExp.from_sampling(nb_obs=1, coordinates=self.coordinates, margin_model=self.margin_model) - self.assertAlmostEqual(obs.df_maxima_gev.mean()[0], 4.692385276235156) + self.assertAlmostEqual(obs.df_maxima_gev.mean()[0], 0.574829320536985) def test_annual_maxima_observations_from_daily_observations(self): obs = DailyExpAnnualMaxima.from_sampling(1, self.coordinates, self.margin_model) - self.assertAlmostEqual(obs.df_maxima_gev.mean()[0], 1183.8468374768636) + self.assertAlmostEqual(obs.df_maxima_gev.mean()[0], 6.523848726794694) if __name__ == '__main__': -- GitLab