diff --git a/extreme_estimator/estimator/full_estimator.py b/extreme_estimator/estimator/full_estimator.py
index d6df3dd7aba67c2e1efbc4bfd24a576970e2ac38..bd3e1618185eadd8a3f0983a7911ac42101e412d 100644
--- a/extreme_estimator/estimator/full_estimator.py
+++ b/extreme_estimator/estimator/full_estimator.py
@@ -1,4 +1,7 @@
from extreme_estimator.extreme_models.margin_model.abstract_margin_model import AbstractMarginModel
+from extreme_estimator.extreme_models.margin_model.margin_function.independent_margin_function import \
+ LinearMarginFunction
+from extreme_estimator.extreme_models.margin_model.smooth_margin_model import LinearMarginModel
from extreme_estimator.extreme_models.max_stable_model.abstract_max_stable_model import AbstractMaxStableModel
from extreme_estimator.estimator.abstract_estimator import AbstractEstimator
from extreme_estimator.estimator.margin_estimator import SmoothMarginEstimator
@@ -12,7 +15,7 @@ class AbstractFullEstimator(AbstractEstimator):
class SmoothMarginalsThenUnitaryMsp(AbstractFullEstimator):
- def __init__(self, dataset: AbstractDataset, margin_model: AbstractMarginModel,
+ def __init__(self, dataset: AbstractDataset, margin_model: LinearMarginModel,
max_stable_model: AbstractMaxStableModel):
super().__init__(dataset)
# Instantiate the two associated estimators
@@ -36,9 +39,25 @@ class FullEstimatorInASingleStep(AbstractFullEstimator):
pass
-class FullEstimatorInASingleStepWithSmoothMarginals(AbstractFullEstimator):
+class FullEstimatorInASingleStepWithSmoothMargin(AbstractFullEstimator):
"""The method of Gaume, check if its method is in a single step or not"""
- pass
+
+ def __init__(self, dataset: AbstractDataset, margin_model: LinearMarginModel,
+ max_stable_model: AbstractMaxStableModel):
+ super().__init__(dataset)
+ self.max_stable_model = max_stable_model
+ self.smooth_margin_function_to_fit = margin_model.margin_function_start_fit
+
+ def _fit(self):
+ # todo: specify the shape of the margin
+ # Estimate the margin
+ self.max_stable_params_fitted = self.max_stable_model.fitmaxstab(
+ maxima_frech=self.dataset.maxima_frech,
+ df_coordinates=self.dataset.df_coordinates,
+ fit_marge=True,
+ fit_marge_form_dict=self.smooth_margin_function_to_fit.fit_marge_form_dict,
+ margin_start_dict=self.smooth_margin_function_to_fit.margin_start_dict
+ )
class PointwiseAndThenUnitaryMsp(AbstractFullEstimator):
diff --git a/extreme_estimator/estimator/max_stable_estimator.py b/extreme_estimator/estimator/max_stable_estimator.py
index b663ac670d7650d4740639eeaf72722f61c8233c..4ca92d18ee48b44a5bee8d581326b78d19a231ea 100644
--- a/extreme_estimator/estimator/max_stable_estimator.py
+++ b/extreme_estimator/estimator/max_stable_estimator.py
@@ -19,7 +19,8 @@ class MaxStableEstimator(AbstractMaxStableEstimator):
assert self.dataset.maxima_frech is not None
self.max_stable_params_fitted = self.max_stable_model.fitmaxstab(
maxima_frech=self.dataset.maxima_frech,
- coord=self.dataset.coordinates)
+ df_coordinates=self.dataset.spatial_coordinates.df_coordinates)
+ # coord=self.dataset.coordinates)
def _error(self, true_max_stable_params: dict):
absolute_errors = {param_name: np.abs(param_true_value - self.max_stable_params_fitted[param_name])
diff --git a/extreme_estimator/extreme_models/margin_model/abstract_margin_model.py b/extreme_estimator/extreme_models/margin_model/abstract_margin_model.py
index fea9bfccbc5e6fa747f92c57e24a5ae3ae9d9621..4d4f62625a0f70948521894435857ae920ce4bcd 100644
--- a/extreme_estimator/extreme_models/margin_model/abstract_margin_model.py
+++ b/extreme_estimator/extreme_models/margin_model/abstract_margin_model.py
@@ -28,6 +28,7 @@ class AbstractMarginModel(AbstractModel):
@classmethod
def convert_maxima(cls, convertion_r_function, maxima: np.ndarray, coordinates: np.ndarray,
margin_function: AbstractMarginFunction):
+ assert isinstance(coordinates, np.ndarray)
assert len(maxima) == len(coordinates)
converted_maxima = []
for x, coordinate in zip(maxima, coordinates):
@@ -58,9 +59,7 @@ class AbstractMarginModel(AbstractModel):
maxima_gev.append(x_gev)
return np.array(maxima_gev)
- # Fitting methods
+ # Fitting methods needs to be defined in child classes
def fitmargin_from_maxima_gev(self, maxima_gev: np.ndarray, coordinates: np.ndarray) -> AbstractMarginFunction:
pass
-
- # Define the method to sample/fit a single gev
diff --git a/extreme_estimator/extreme_models/margin_model/margin_function/independent_margin_function.py b/extreme_estimator/extreme_models/margin_model/margin_function/independent_margin_function.py
index 18fbd21ecd9fe8b6bc2ae2b55b7950cd017f1ee0..a9f416fa3228082c02985e0529dd00fdbed30d18 100644
--- a/extreme_estimator/extreme_models/margin_model/margin_function/independent_margin_function.py
+++ b/extreme_estimator/extreme_models/margin_model/margin_function/independent_margin_function.py
@@ -41,7 +41,7 @@ class ConstantParamFunction(ParamFunction):
class LinearOneAxisParamFunction(ParamFunction):
- def __init__(self, linear_axis, coordinates_axis, start, end=0.0):
+ def __init__(self, linear_axis: int, coordinates_axis: np.ndarray, start: float, end: float = 0.0):
self.linear_axis = linear_axis
self.t_min = coordinates_axis.min()
self.t_max = coordinates_axis.max()
@@ -63,6 +63,10 @@ class LinearMarginFunction(IndependentMarginFunction):
def __init__(self, spatial_coordinates: AbstractSpatialCoordinates, default_params: GevParams,
gev_param_name_to_linear_axis: Dict[str, int]):
super().__init__(spatial_coordinates, default_params)
+
+ # By default, set a constant function for all the params
+ self.gev_param_to_str_form = {gev_param_name: '1' for gev_param_name in GevParams.GEV_PARAM_NAMES}
+
self.param_to_linear_dims = gev_param_name_to_linear_axis
assert all([axis < np.ndim(spatial_coordinates.coordinates) for axis in gev_param_name_to_linear_axis.values()])
# Initialize gev_parameter_to_param_function
@@ -76,3 +80,22 @@ class LinearMarginFunction(IndependentMarginFunction):
param_function = LinearOneAxisParamFunction(linear_axis=linear_axis, coordinates_axis=coordinates_axis,
start=self.default_params[gev_param_name])
self.gev_param_name_to_param_function[gev_param_name] = param_function
+
+
+ @property
+ def fit_marge_form_dict(self):
+ """
+ Example:
+ loc.form = loc ~ E
+ scale.form = scale ~ N
+ shape.form = shape ~ E+N
+ :return:
+ """
+ # todo specify all functions
+ return {gev_param_name + '.form': gev_param_name + ' ~ ' + self.gev_param_to_str_form[gev_param_name]
+ for gev_param_name in GevParams.GEV_PARAM_NAMES}
+
+ @property
+ def margin_start_dict(self):
+ """Example for the names: locCoeff1 locCoeff2 """
+ return {gev_param_name + 'Coeff1': 0 for gev_param_name in GevParams.GEV_PARAM_NAMES}
diff --git a/extreme_estimator/extreme_models/margin_model/smooth_margin_model.py b/extreme_estimator/extreme_models/margin_model/smooth_margin_model.py
index 6a33d939a567af114610bfc45b2fba0653d5c9ff..b15b8c6dac1c75fc6b4d68e839d13523bb0a30ea 100644
--- a/extreme_estimator/extreme_models/margin_model/smooth_margin_model.py
+++ b/extreme_estimator/extreme_models/margin_model/smooth_margin_model.py
@@ -18,6 +18,9 @@ class LinearMarginModel(AbstractMarginModel):
default_params=GevParams.from_dict(self.params_start_fit),
gev_param_name_to_linear_axis=gev_param_name_to_linear_axis)
+ def fitmargin_from_maxima_gev(self, maxima_gev: np.ndarray, coordinates: np.ndarray) -> AbstractMarginFunction:
+ return self.margin_function_start_fit
+
class ConstantMarginModel(LinearMarginModel):
@@ -33,8 +36,8 @@ class LinearShapeAxis0MarginModel(LinearMarginModel):
def load_margin_functions(self, margin_function_class: type = None, gev_param_name_to_linear_axis=None):
super().load_margin_functions({GevParams.GEV_SHAPE: 0})
- def fitmargin_from_maxima_gev(self, maxima_gev: np.ndarray, coordinates: np.ndarray) -> AbstractMarginFunction:
- return self.margin_function_start_fit
+ # def fitmargin_from_maxima_gev(self, maxima_gev: np.ndarray, coordinates: np.ndarray) -> AbstractMarginFunction:
+ # pass
if __name__ == '__main__':
diff --git a/extreme_estimator/extreme_models/max_stable_model/abstract_max_stable_model.py b/extreme_estimator/extreme_models/max_stable_model/abstract_max_stable_model.py
index 8c5a5f31e59ca89d7c3442cc7ea74bfc127ea401..1043f5bc51a3c9879a57ce9f26341b3e9e0323fc 100644
--- a/extreme_estimator/extreme_models/max_stable_model/abstract_max_stable_model.py
+++ b/extreme_estimator/extreme_models/max_stable_model/abstract_max_stable_model.py
@@ -1,7 +1,9 @@
+import pandas as pd
from enum import Enum
import numpy as np
import rpy2
+import rpy2.robjects as ro
from rpy2.robjects import ListVector
from extreme_estimator.extreme_models.abstract_model import AbstractModel
@@ -17,20 +19,34 @@ class AbstractMaxStableModel(AbstractModel):
def cov_mod_param(self):
return {'cov.mod': self.cov_mod}
- def fitmaxstab(self, maxima_frech: np.ndarray, coord: np.ndarray, fit_marge=False):
- assert all([isinstance(arr, np.ndarray) for arr in [maxima_frech, coord]])
+ def fitmaxstab(self, maxima_frech: np.ndarray, df_coordinates: pd.DataFrame, fit_marge=False,
+ fit_marge_form_dict=None, margin_start_dict=None):
+ assert isinstance(maxima_frech, np.ndarray)
+ assert isinstance(df_coordinates, pd.DataFrame)
+ if fit_marge:
+ assert fit_marge_form_dict is not None
+ assert margin_start_dict is not None
+
+ # Add the colnames to df_coordinates DataFrame to enable specification of the margin functions
+ df_coordinates = df_coordinates.copy()
+ df_coordinates.colnames = ro.StrVector(list(df_coordinates.columns))
+ # Transform the formula string representation into robjects.Formula("y ~ x")
# Specify the fit params
- fit_params = {
- 'fit.marge': fit_marge,
- 'start': self.r.list(**self.params_start_fit),
- }
+ fit_params = self.cov_mod_param.copy()
+ start_dict = self.params_start_fit
+ if fit_marge:
+ start_dict.update(margin_start_dict)
+ fit_params.update({k: ro.Formula(v) for k, v in fit_marge_form_dict.items()})
+ fit_params['start'] = self.r.list(**start_dict)
+ fit_params['fit.marge'] = fit_marge
# Run the fitmaxstab in R
# todo: find how to specify the optim function to use
+ coord = df_coordinates.values
+
try:
- res = self.r.fitmaxstab(np.transpose(maxima_frech), coord, **self.cov_mod_param,
- **fit_params) # type: ListVector
+ res = self.r.fitmaxstab(data=np.transpose(maxima_frech), coord=coord, **fit_params) # type: ListVector
except rpy2.rinterface.RRuntimeError as error:
- raise Exception('Some R exception have been launched at RunTime: {}'.format(error.__repr__()))
+ raise Exception('Some R exception have been launched at RunTime: \n {}'.format(error.__repr__()))
# todo: maybe if the convergence was not successful I could try other starting point several times
# Retrieve the resulting fitted values
fitted_values = res.rx2('fitted.values')
diff --git a/extreme_estimator/extreme_models/max_stable_model/max_stable_fit.R b/extreme_estimator/extreme_models/max_stable_model/max_stable_fit.R
index 5c0491e8b4b948f6df4f378097f7c7ab4fc91c2f..4daba7f105b904c578c7df2d1f429feced3f5235 100644
--- a/extreme_estimator/extreme_models/max_stable_model/max_stable_fit.R
+++ b/extreme_estimator/extreme_models/max_stable_model/max_stable_fit.R
@@ -9,24 +9,29 @@ if (call_main) {
n.obs = 50
n.site = 2
coord <- matrix(rnorm(2*n.site, sd = sqrt(.2)), ncol = 2)
+ colnames(coord) = c("E", "N")
# Generate the data
data <- rmaxstab(n.obs, coord, "gauss", cov11 = 100, cov12 = 25, cov22 = 220)
# data <- rmaxstab(n.obs, coord, "brown", range = 3, smooth = 0.5)
# data <- rmaxstab(n.obs, coord, "whitmat", nugget = 0.0, range = 3, smooth = 0.5)
# Fit back the data
- # print(data)
+ # print(data)n
# res = fitmaxstab(data, coord, "gauss", fit.marge=FALSE, )
# res = fitmaxstab(data, coord, "brown")
# res = fitmaxstab(data, coord, "whitmat", start=)
- namedlist = list(cov11 = 1.0, cov12 = 1.2, cov22 = 2.2)
- res = fitmaxstab(data=data, coord=coord, cov.mod="gauss", start=namedlist)
- # res = fitmaxstab(data, coord, "whitmat", par())
- # print(res)
- # print(class(res))
- gev2frech(x, loc, scale, shape, emp=FALSE)
- frech2gev(x, loc, scale, shape)
- # print(names(res))
+
+ loc.form = loc ~ 1
+ scale.form = scale ~ 1
+ shape.form = shape ~ 1
+
+
+ namedlist = list(cov11 = 1.0, cov12 = 1.2, cov22 = 2.2, locCoeff1=1.0, scaleCoeff1=1.0, shapeCoeff1=1.0)
+ res = fitmaxstab(data=data, coord=coord, cov.mod="gauss", start=namedlist, fit.marge=TRUE, loc.form=loc.form, scale.form=scale.form,shape.form=shape.form)
+
+ # namedlist = list(cov11 = 1.0, cov12 = 1.2, cov22 = 2.2)
+ # res = fitmaxstab(data=data, coord=coord, cov.mod="gauss", start=namedlist)
+
for (name in names(res)){
print(name)
print(res[name])
diff --git a/extreme_estimator/extreme_models/max_stable_model/max_stable_fit.py b/extreme_estimator/extreme_models/max_stable_model/max_stable_fit.py
new file mode 100644
index 0000000000000000000000000000000000000000..cfb879503863c680f066a1d627c6d40d84f96390
--- /dev/null
+++ b/extreme_estimator/extreme_models/max_stable_model/max_stable_fit.py
@@ -0,0 +1,57 @@
+# test how to call the max stable method
+import pandas as pd
+
+import numpy as np
+
+from extreme_estimator.extreme_models.utils import get_loaded_r
+from extreme_estimator.gev.gev_mle_fit import GevMleFit
+import rpy2.robjects.numpy2ri as rpyn
+
+
+import rpy2.robjects as robjects
+
+def max_stable_fit():
+ r = get_loaded_r()
+ r("""
+ set.seed(42)
+ n.obs = 50
+ n.site = 2
+ coord <- matrix(rnorm(2*n.site, sd = sqrt(.2)), ncol = 2)
+ colnames(coord) = c("E", "N")
+
+ # Generate the data
+ data <- rmaxstab(n.obs, coord, "gauss", cov11 = 100, cov12 = 25, cov22 = 220)
+
+
+ loc.form = loc ~ 1
+ scale.form = scale ~ 1
+ shape.form = shape ~ 1
+
+ namedlist = list(cov11 = 1.0, cov12 = 1.2, cov22 = 2.2, locCoeff1=1.0, scaleCoeff1=1.0, shapeCoeff1=1.0)
+ # res = fitmaxstab(data=data, coord=coord, cov.mod="gauss", start=namedlist, fit.marge=TRUE, loc.form=loc.form, scale.form=scale.form,shape.form=shape.form)
+ """)
+
+ # coord = np.array(r['coord'])
+ data = r['data']
+ coord = pd.DataFrame(r['coord'])
+ coord.colnames = robjects.StrVector(['E', 'N'])
+
+ print(r['loc.form'])
+ print(type(r['loc.form']))
+ # x_gev = rpyn.numpy2ri(x_gev)
+ print(coord)
+
+ print(coord.colnames)
+
+ # res = r.fitmaxstab(data=data, coord=coord, covmod="gauss", start=namedlist, fit.marge=TRUE, loc.form=loc.form, scale.form=scale.form,shape.form=shape.form)
+
+
+
+ # m2.colnames = R.StrVector(['x', 'y'])
+ # res = r.fitmaxstab()
+ # mle_params = dict(r.attr(res, 'coef').items())
+ # print(mle_params)
+
+
+if __name__ == '__main__':
+ max_stable_fit()
diff --git a/extreme_estimator/extreme_models/utils.py b/extreme_estimator/extreme_models/utils.py
index 87a89fb947fb1504a02892e901c285e5c25441a6..513d50c7506a84f768a2a33082d9065a9b66f26f 100644
--- a/extreme_estimator/extreme_models/utils.py
+++ b/extreme_estimator/extreme_models/utils.py
@@ -1,13 +1,14 @@
import os.path as op
+
import rpy2.robjects as ro
-import pandas as pd
-import numpy as np
-import rpy2.robjects.numpy2ri as npr
+from rpy2.robjects import numpy2ri
+from rpy2.robjects import pandas2ri
def get_loaded_r() -> ro.R:
r = ro.r
- ro.numpy2ri.activate()
+ numpy2ri.activate()
+ pandas2ri.activate()
r.library('SpatialExtremes')
return r
diff --git a/spatio_temporal_dataset/dataset/abstract_dataset.py b/spatio_temporal_dataset/dataset/abstract_dataset.py
index 33faaa6d5771fe11fa0618bd25b4b1be1b580440..a63ed182d556fd27c8202fe9774d6fa92878a19a 100644
--- a/spatio_temporal_dataset/dataset/abstract_dataset.py
+++ b/spatio_temporal_dataset/dataset/abstract_dataset.py
@@ -33,6 +33,10 @@ class AbstractDataset(object):
# Merge dataframes with the maxima and with the coordinates
return self.temporal_observations.df_maxima_gev.join(self.spatial_coordinates.df_coordinates)
+ @property
+ def df_coordinates(self):
+ return self.spatial_coordinates.df_coordinates
+
@property
def coordinates(self):
return self.spatial_coordinates.coordinates
diff --git a/test/test_extreme_estimator/test_estimator/test_full_estimators.py b/test/test_extreme_estimator/test_estimator/test_full_estimators.py
index 49df744f0ded4bbefef5f84e2825a2e02bf2e4d4..938462e19b31ab450044b5e462dfc17702480593 100644
--- a/test/test_extreme_estimator/test_estimator/test_full_estimators.py
+++ b/test/test_extreme_estimator/test_estimator/test_full_estimators.py
@@ -1,7 +1,8 @@
import unittest
from itertools import product
-from extreme_estimator.estimator.full_estimator import SmoothMarginalsThenUnitaryMsp
+from extreme_estimator.estimator.full_estimator import SmoothMarginalsThenUnitaryMsp, \
+ FullEstimatorInASingleStepWithSmoothMargin
from spatio_temporal_dataset.dataset.simulation_dataset import FullSimulatedDataset
from spatio_temporal_dataset.spatial_coordinates.generated_coordinates import CircleCoordinatesRadius1
from test.test_extreme_estimator.test_estimator.test_margin_estimators import TestSmoothMarginEstimator
@@ -10,16 +11,16 @@ from test.test_extreme_estimator.test_estimator.test_max_stable_estimators impor
class TestFullEstimators(unittest.TestCase):
DISPLAY = False
- FULL_ESTIMATORS = [SmoothMarginalsThenUnitaryMsp]
+ FULL_ESTIMATORS = [SmoothMarginalsThenUnitaryMsp, FullEstimatorInASingleStepWithSmoothMargin][:]
def setUp(self):
super().setUp()
self.spatial_coordinates = CircleCoordinatesRadius1.from_nb_points(nb_points=5, max_radius=1)
self.max_stable_models = TestMaxStableEstimators.load_max_stable_models()
- self.margin_models = TestSmoothMarginEstimator.load_margin_models(spatial_coordinates=self.spatial_coordinates)
+ self.smooth_margin_models = TestSmoothMarginEstimator.load_smooth_margin_models(spatial_coordinates=self.spatial_coordinates)
def test_full_estimators(self):
- for margin_model, max_stable_model in product(self.margin_models, self.max_stable_models):
+ for margin_model, max_stable_model in product(self.smooth_margin_models, self.max_stable_models):
dataset = FullSimulatedDataset.from_double_sampling(nb_obs=10, margin_model=margin_model,
spatial_coordinates=self.spatial_coordinates,
max_stable_model=max_stable_model)
diff --git a/test/test_extreme_estimator/test_estimator/test_margin_estimators.py b/test/test_extreme_estimator/test_estimator/test_margin_estimators.py
index 9dcea3dcaef2861970b3cf0e4244f8bc0f70b94c..fe77b35ac48b9bcaa19c5f8faced344a5b8b345c 100644
--- a/test/test_extreme_estimator/test_estimator/test_margin_estimators.py
+++ b/test/test_extreme_estimator/test_estimator/test_margin_estimators.py
@@ -12,19 +12,19 @@ from spatio_temporal_dataset.spatial_coordinates.generated_coordinates import Ci
class TestSmoothMarginEstimator(unittest.TestCase):
DISPLAY = False
MARGIN_TYPES = [ConstantMarginModel, LinearShapeAxis0MarginModel][1:]
- MARGIN_ESTIMATORS = [SmoothMarginEstimator]
+ SMOOTH_MARGIN_ESTIMATORS = [SmoothMarginEstimator]
def setUp(self):
super().setUp()
self.spatial_coordinates = CircleCoordinatesRadius1.from_nb_points(nb_points=5, max_radius=1)
- self.margin_models = self.load_margin_models(spatial_coordinates=self.spatial_coordinates)
+ self.smooth_margin_models = self.load_smooth_margin_models(spatial_coordinates=self.spatial_coordinates)
@classmethod
- def load_margin_models(cls, spatial_coordinates):
+ def load_smooth_margin_models(cls, spatial_coordinates):
return [margin_class(spatial_coordinates=spatial_coordinates) for margin_class in cls.MARGIN_TYPES]
def test_dependency_estimators(self):
- for margin_model in self.margin_models:
+ for margin_model in self.smooth_margin_models:
dataset = MarginDataset.from_sampling(nb_obs=10, margin_model=margin_model,
spatial_coordinates=self.spatial_coordinates)
# Fit estimator