diff --git a/extreme_estimator/R_fit/max_stable_fit/abstract_max_stable_model.py b/extreme_estimator/R_fit/max_stable_fit/abstract_max_stable_model.py
index 75a4c3cccb823ee14ffb1472277d3226fd6dfd65..bd2c3723d7275024f4039062afb38ca15d73941d 100644
--- a/extreme_estimator/R_fit/max_stable_fit/abstract_max_stable_model.py
+++ b/extreme_estimator/R_fit/max_stable_fit/abstract_max_stable_model.py
@@ -1,3 +1,4 @@
+import rpy2
from enum import Enum
from rpy2.robjects import ListVector
@@ -15,19 +16,20 @@ class AbstractMaxStableModel(object):
self.user_params_sample = params_sample
self.r = get_loaded_r()
- def fitmaxstab(self, maxima: np.ndarray, coord: np.ndarray, ):
- # todo: find how to specify a startign point, understand how is it set by default
- # res = None
- # tries = 5
- # nb_tries = 0
- # while res is None and nb_tries < tries:
- # try:
- # res = self.r.fitmaxstab(maxima, coord, **self.cov_mod_param) # type: ListVector
- # except rpy2.rinterface.RRuntimeError:
- # pass
- # nb_tries += 1
- # print(nb_tries)
- res = self.r.fitmaxstab(np.transpose(maxima), coord, **self.cov_mod_param) # type: ListVector
+ def fitmaxstab(self, maxima: np.ndarray, coord: np.ndarray, fit_marge=False):
+ # Specify the fit params
+ fit_params = {
+ 'fit.marge': fit_marge,
+ 'start': self.r.list(**self.params_start_fit),
+ }
+ # Run the fitmaxstab in R
+ # todo: find how to specify the optim function to use
+ try:
+ res = self.r.fitmaxstab(np.transpose(maxima), coord, **self.cov_mod_param, **fit_params) # type: ListVector
+ except rpy2.rinterface.RRuntimeError as error:
+ raise Exception('Some R exception have been launched at RunTime: {}'.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')
fitted_values = {key: fitted_values.rx2(key)[0] for key in fitted_values.names}
return fitted_values
@@ -36,7 +38,7 @@ class AbstractMaxStableModel(object):
"""
Return an numpy of maxima. With rows being the stations and columns being the years of maxima
"""
- maxima = np.array(self.r.rmaxstab(nb_obs, coord, **self.cov_mod_param, **self.params_sample))
+ maxima = np.array(self.r.rmaxstab(nb_obs, coord, *list(self.cov_mod_param.values()), **self.params_sample))
return np.transpose(maxima)
@property
@@ -53,6 +55,7 @@ class AbstractMaxStableModel(object):
@staticmethod
def merge_params(default_params, input_params):
+ assert default_params is not None, 'some default_params need to be specified'
merged_params = default_params.copy()
if input_params is not None:
assert isinstance(default_params, dict) and isinstance(input_params, dict)
diff --git a/extreme_estimator/R_fit/max_stable_fit/max_stable_fit.R b/extreme_estimator/R_fit/max_stable_fit/max_stable_fit.R
index 287a186b83a731c2ff6b5aea44975952bc1b9dd9..5a3754e9a9cad8b83c1ba62256d693f83e04a789 100644
--- a/extreme_estimator/R_fit/max_stable_fit/max_stable_fit.R
+++ b/extreme_estimator/R_fit/max_stable_fit/max_stable_fit.R
@@ -10,22 +10,26 @@ if (call_main) {
n.site = 2
coord <- matrix(rnorm(2*n.site, sd = sqrt(.2)), ncol = 2)
- print(coord)
# Generate the data
- # data <- rmaxstab(n.obs, coord, "gauss", cov11 = 100, cov12 = 25, cov22 = 220)
+ 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)
- data <- rmaxstab(n.obs, coord, "cauchy", )
# Fit back the data
- print(data)
+ # print(data)
# res = fitmaxstab(data, coord, "gauss", fit.marge=FALSE, )
# res = fitmaxstab(data, coord, "brown")
- res = fitmaxstab(data, coord, "cauchy")
- # res = fitmaxstab(data, coord, "gauss", start=list(0,0,0))
- print(res)
- print(class(res))
- print(names(res))
+ # 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))
+ # print(names(res))
+ for (name in names(res)){
+ print(name)
+ print(res[name])
+ }
print(res['fitted.values'])
- print(res['convergence'])
+ # print(res['convergence'])
}
diff --git a/extreme_estimator/R_fit/max_stable_fit/max_stable_models.py b/extreme_estimator/R_fit/max_stable_fit/max_stable_models.py
index 2c77573f200a818e725e7c56147c0cea4d04e3d9..683899cea586068a5175a8d416820c15a002c42d 100644
--- a/extreme_estimator/R_fit/max_stable_fit/max_stable_models.py
+++ b/extreme_estimator/R_fit/max_stable_fit/max_stable_models.py
@@ -9,12 +9,12 @@ class Smith(AbstractMaxStableModel):
def __init__(self, params_start_fit=None, params_sample=None):
super().__init__(params_start_fit=params_start_fit, params_sample=params_sample)
self.cov_mod = 'gauss'
- self.default_params_start_fit = {}
- self.default_params_sample = {
- 'cov11': 100,
- 'cov12': 25,
- 'cov22': 220,
+ self.default_params_start_fit = {
+ 'cov11': 1,
+ 'cov12': 0,
+ 'cov22': 1
}
+ self.default_params_sample = self.default_params_start_fit.copy()
class BrownResnick(AbstractMaxStableModel):
@@ -22,7 +22,10 @@ class BrownResnick(AbstractMaxStableModel):
def __init__(self, params_start_fit=None, params_sample=None):
super().__init__(params_start_fit=params_start_fit, params_sample=params_sample)
self.cov_mod = 'brown'
- self.default_params_start_fit = {}
+ self.default_params_start_fit = {
+ 'range': 3,
+ 'smooth': 0.5,
+ }
self.default_params_sample = {
'range': 3,
'smooth': 0.5,
@@ -35,6 +38,7 @@ class Schlather(AbstractMaxStableModelWithCovarianceFunction):
super().__init__(params_start_fit, params_sample, covariance_function)
self.cov_mod = self.covariance_function.name
self.default_params_sample.update({})
+ self.default_params_start_fit = self.default_params_sample.copy()
class Geometric(AbstractMaxStableModelWithCovarianceFunction):
@@ -43,6 +47,7 @@ class Geometric(AbstractMaxStableModelWithCovarianceFunction):
super().__init__(params_start_fit, params_sample, covariance_function)
self.cov_mod = 'g' + self.covariance_function.name
self.default_params_sample .update({'sigma2': 0.5})
+ self.default_params_start_fit = self.default_params_sample.copy()
class ExtremalT(AbstractMaxStableModelWithCovarianceFunction):
@@ -51,6 +56,7 @@ class ExtremalT(AbstractMaxStableModelWithCovarianceFunction):
super().__init__(params_start_fit, params_sample, covariance_function)
self.cov_mod = 't' + self.covariance_function.name
self.default_params_sample .update({'DoF': 2})
+ self.default_params_start_fit = self.default_params_sample.copy()
class ISchlather(AbstractMaxStableModelWithCovarianceFunction):
@@ -59,3 +65,4 @@ class ISchlather(AbstractMaxStableModelWithCovarianceFunction):
super().__init__(params_start_fit, params_sample, covariance_function)
self.cov_mod = 'i' + self.covariance_function.name
self.default_params_sample .update({'alpha': 0.5})
+ self.default_params_start_fit = self.default_params_sample.copy()
diff --git a/test/extreme_estimator/test_max_stable_fit.py b/test/extreme_estimator/test_max_stable_fit.py
index cca832bfa3e1454d063df35c13550882867cd5ca..662c39ddee633451a698bd3931440ab6fb9169bb 100644
--- a/test/extreme_estimator/test_max_stable_fit.py
+++ b/test/extreme_estimator/test_max_stable_fit.py
@@ -16,7 +16,6 @@ class TestMaxStableFit(unittest.TestCase):
self.spatial_coord = CircleCoordinates.from_nb_points(nb_points=5, max_radius=1)
self.max_stable_models = []
for max_stable_class in self.MAX_STABLE_CLASSES:
- print(max_stable_class)
if issubclass(max_stable_class, AbstractMaxStableModelWithCovarianceFunction):
self.max_stable_models.extend([max_stable_class(covariance_function=covariance_function)
for covariance_function in CovarianceFunction])
@@ -27,9 +26,11 @@ class TestMaxStableFit(unittest.TestCase):
for max_stable_model in self.max_stable_models:
dataset = SimulatedDataset.from_max_stable_sampling(nb_obs=10, max_stable_model=max_stable_model,
spatial_coordinates=self.spatial_coord)
+ fitted_values = max_stable_model.fitmaxstab(maxima=dataset.maxima, coord=dataset.coord)
if display:
- print(dataset.head())
- max_stable_model.fitmaxstab(maxima=dataset.maxima, coord=dataset.coord)
+ print(type(max_stable_model))
+ print(dataset.df_dataset.head())
+ print(fitted_values)
self.assertTrue(True)