diff --git a/experiment/fit_diagnosis/split_curve.py b/experiment/fit_diagnosis/split_curve.py
index a7511791734f5ea4503025f602a2608b0342e377..8dfebf30931573b25bd6db470fc6d28ddc7ece8d 100644
--- a/experiment/fit_diagnosis/split_curve.py
+++ b/experiment/fit_diagnosis/split_curve.py
@@ -73,10 +73,12 @@ class SplitCurve(object):
         # Create bins of data, each with an associated color corresponding to its error
 
         data = self.mean_error_dict[gev_value_name].values
+        data_min, data_max = data.min(), data.max()
         nb_bins = 10
-        limits = np.linspace(data.min(), data.max(), num=nb_bins + 1)
+        limits = np.linspace(data_min, data_max, num=nb_bins + 1)
         limits[-1] += 0.01
-        colors = cm.binary(limits)
+        # Binary color should
+        colors = cm.binary((limits - data_min/ (data_max - data_min)))
 
         # Display train/test points
         for split, marker in [(self.dataset.train_split, 'o'), (self.dataset.test_split, 'x')]:
diff --git a/experiment/fit_diagnosis/split_curve_alps.py b/experiment/fit_diagnosis/split_curve_alps.py
new file mode 100644
index 0000000000000000000000000000000000000000..13fbd82839d9de213a71d4be5008c4d24b372610
--- /dev/null
+++ b/experiment/fit_diagnosis/split_curve_alps.py
@@ -0,0 +1,46 @@
+from experiment.fit_diagnosis.split_curve import SplitCurve
+from extreme_estimator.estimator.full_estimator import FullEstimatorInASingleStepWithSmoothMargin
+from extreme_estimator.estimator.margin_estimator import SmoothMarginEstimator
+from extreme_estimator.extreme_models.margin_model.smooth_margin_model import ConstantMarginModel, \
+    LinearAllParametersAllDimsMarginModel
+from extreme_estimator.extreme_models.max_stable_model.max_stable_models import Smith
+from extreme_estimator.gev_params import GevParams
+from spatio_temporal_dataset.coordinates.spatial_coordinates.alps_station_2D_coordinates import \
+    AlpsStation2DCoordinates, AlpsStation2DCoordinatesBetweenZeroAndOne
+from spatio_temporal_dataset.coordinates.spatial_coordinates.coordinates_1D import LinSpaceSpatialCoordinates
+from spatio_temporal_dataset.dataset.simulation_dataset import FullSimulatedDataset
+
+
+class SplitCurveAlps(SplitCurve):
+
+    def __init__(self, nb_fit: int = 1):
+        super().__init__(nb_fit)
+        self.nb_obs = 2
+        # nb_points = len(AlpsStation2DCoordinates.from_csv())
+        nb_points = 10
+        self.coordinates = AlpsStation2DCoordinatesBetweenZeroAndOne.from_nb_points(nb_points, train_split_ratio=0.8)
+        # MarginModel Linear with respect to the shape (from 0.01 to 0.02)
+        params_sample = {
+            (GevParams.GEV_LOC, 0): 10,
+            (GevParams.GEV_SHAPE, 0): 1.0,
+            (GevParams.GEV_SCALE, 0): 10,
+        }
+        self.margin_model = ConstantMarginModel(coordinates=self.coordinates, params_sample=params_sample)
+        self.max_stable_model = Smith()
+
+    def load_dataset(self):
+        return FullSimulatedDataset.from_double_sampling(nb_obs=self.nb_obs, margin_model=self.margin_model,
+                                                         coordinates=self.coordinates,
+                                                         max_stable_model=self.max_stable_model)
+
+    def load_estimator(self, dataset):
+        max_stable_model = Smith()
+        margin_model_for_estimator = LinearAllParametersAllDimsMarginModel(dataset.coordinates)
+        # estimator = FullEstimatorInASingleStepWithSmoothMargin(dataset, margin_model_for_estimator, max_stable_model)
+        estimator = SmoothMarginEstimator(dataset, margin_model_for_estimator)
+        return estimator
+
+
+if __name__ == '__main__':
+    curve = SplitCurveAlps(nb_fit=2)
+    curve.fit()
diff --git a/experiment/fit_diagnosis/split_curve_example.py b/experiment/fit_diagnosis/split_curve_example.py
index 96e7c9bd2a6adc44e2197892fd533835081e7140..59732de5c62eca2987e12b234b9eb3ab1e740907 100644
--- a/experiment/fit_diagnosis/split_curve_example.py
+++ b/experiment/fit_diagnosis/split_curve_example.py
@@ -1,5 +1,6 @@
 from experiment.fit_diagnosis.split_curve import SplitCurve
 from extreme_estimator.estimator.full_estimator import FullEstimatorInASingleStepWithSmoothMargin
+from extreme_estimator.estimator.margin_estimator import SmoothMarginEstimator
 from extreme_estimator.extreme_models.margin_model.smooth_margin_model import ConstantMarginModel, \
     LinearAllParametersAllDimsMarginModel
 from extreme_estimator.extreme_models.max_stable_model.max_stable_models import Smith
@@ -32,14 +33,11 @@ class SplitCurveExample(SplitCurve):
     def load_estimator(self, dataset):
         max_stable_model = Smith()
         margin_model_for_estimator = LinearAllParametersAllDimsMarginModel(dataset.coordinates)
-        estimator = FullEstimatorInASingleStepWithSmoothMargin(dataset, margin_model_for_estimator, max_stable_model)
-        # estimator = SmoothMarginEstimator(dataset, margin_model_for_estimator)
+        # estimator = FullEstimatorInASingleStepWithSmoothMargin(dataset, margin_model_for_estimator, max_stable_model)
+        estimator = SmoothMarginEstimator(dataset, margin_model_for_estimator)
         return estimator
 
 
-
-
-
 if __name__ == '__main__':
     curve = SplitCurveExample(nb_fit=2)
     curve.fit()
diff --git a/extreme_estimator/estimator/abstract_estimator.py b/extreme_estimator/estimator/abstract_estimator.py
index 7d5bef88046cb14538b19345682d502e6bf4acc5..246359cb3bcca220e420925a21bfc225ee053e52 100644
--- a/extreme_estimator/estimator/abstract_estimator.py
+++ b/extreme_estimator/estimator/abstract_estimator.py
@@ -1,5 +1,8 @@
 import time
 
+from extreme_estimator.extreme_models.margin_model.margin_function.abstract_margin_function import \
+    AbstractMarginFunction
+from extreme_estimator.extreme_models.margin_model.margin_function.linear_margin_function import LinearMarginFunction
 from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
 from spatio_temporal_dataset.slicer.split import Split
 
@@ -15,10 +18,9 @@ class AbstractEstimator(object):
     def __init__(self, dataset: AbstractDataset):
         self.dataset = dataset  # type: AbstractDataset
         self.additional_information = dict()
-
-    @property
-    def train_split(self):
-        return self.dataset.train_split
+        self._params_fitted = None
+        self._margin_function_fitted = None
+        self._max_stable_model_fitted = None
 
     def fit(self):
         ts = time.time()
@@ -26,6 +28,31 @@ class AbstractEstimator(object):
         te = time.time()
         self.additional_information[self.DURATION] = int((te - ts) * 1000)
 
+    @property
+    def params_fitted(self):
+        assert self.is_fitted
+        return self._params_fitted
+
+    @property
+    def margin_function_fitted(self) -> AbstractMarginFunction:
+        assert self.is_fitted
+        assert self._margin_function_fitted is not None, 'No margin function has been fitted'
+        return self._margin_function_fitted
+
+    def extract_fitted_models_from_fitted_params(self, margin_function_to_fit, full_params_fitted):
+        coef_dict = {k: v for k, v in full_params_fitted.items() if 'Coeff' in k}
+        self._margin_function_fitted = LinearMarginFunction.from_coef_dict(coordinates=self.dataset.coordinates,
+                                                                           gev_param_name_to_linear_dims=margin_function_to_fit.gev_param_name_to_linear_dims,
+                                                                           coef_dict=coef_dict)
+
+    @property
+    def is_fitted(self):
+        return self._params_fitted is not None
+
+    @property
+    def train_split(self):
+        return self.dataset.train_split
+
     def scalars(self, true_max_stable_params: dict):
         error = self._error(true_max_stable_params)
         return {**error, **self.additional_information}
diff --git a/extreme_estimator/estimator/full_estimator.py b/extreme_estimator/estimator/full_estimator.py
index ecfcb66c71224279de9cb838374e0103a3388184..2ec2b86db93a90ffb428147bfe3a6c1c1a891e61 100644
--- a/extreme_estimator/estimator/full_estimator.py
+++ b/extreme_estimator/estimator/full_estimator.py
@@ -14,18 +14,6 @@ class AbstractFullEstimator(AbstractEstimator):
 
     def __init__(self, dataset: AbstractDataset):
         super().__init__(dataset)
-        self._margin_function_fitted = None
-        self._max_stable_model_fitted = None
-
-    @property
-    def margin_function_fitted(self) -> AbstractMarginFunction:
-        assert self._margin_function_fitted is not None, 'Error: estimator has not been fitted'
-        return self._margin_function_fitted
-
-    # @property
-    # def max_stable_fitted(self) -> AbstractMarginFunction:
-    #     assert self._margin_function_fitted is not None, 'Error: estimator has not been fitted'
-    #     return self._margin_function_fitted
 
 
 class SmoothMarginalsThenUnitaryMsp(AbstractFullEstimator):
@@ -43,7 +31,8 @@ class SmoothMarginalsThenUnitaryMsp(AbstractFullEstimator):
         # Compute the maxima_frech
         maxima_gev_train = self.dataset.maxima_gev(split=self.train_split)
         maxima_frech = AbstractMarginModel.gev2frech(maxima_gev=maxima_gev_train,
-                                                     coordinates_values=self.dataset.coordinates_values(self.train_split),
+                                                     coordinates_values=self.dataset.coordinates_values(
+                                                         self.train_split),
                                                      margin_function=self.margin_estimator.margin_function_fitted)
         # Update maxima frech field through the dataset object
         self.dataset.set_maxima_frech(maxima_frech, split=self.train_split)
@@ -68,7 +57,7 @@ class FullEstimatorInASingleStepWithSmoothMargin(AbstractFullEstimator):
 
     def _fit(self):
         # Estimate both the margin and the max-stable structure
-        self.full_params_fitted = self.max_stable_model.fitmaxstab(
+        self._params_fitted = self.max_stable_model.fitmaxstab(
             maxima_gev=self.dataset.maxima_gev(split=self.train_split),
             df_coordinates=self.dataset.df_coordinates(split=self.train_split),
             fit_marge=True,
@@ -76,10 +65,7 @@ class FullEstimatorInASingleStepWithSmoothMargin(AbstractFullEstimator):
             margin_start_dict=self.linear_margin_function_to_fit.coef_dict
         )
         # Create the fitted margin function
-        coef_dict = {k: v for k, v in self.full_params_fitted.items() if 'Coeff' in k}
-        self._margin_function_fitted = LinearMarginFunction.from_coef_dict(coordinates=self.dataset.coordinates,
-                                                                           gev_param_name_to_linear_dims=self.linear_margin_function_to_fit.gev_param_name_to_linear_dims,
-                                                                           coef_dict=coef_dict)
+        self.extract_fitted_models_from_fitted_params(self.linear_margin_function_to_fit, self._params_fitted)
 
 
 class PointwiseAndThenUnitaryMsp(AbstractFullEstimator):
diff --git a/extreme_estimator/estimator/margin_estimator.py b/extreme_estimator/estimator/margin_estimator.py
index 829cd65c7305c51cfa8421ba8392a98eb3b89250..a68539d769512324bed5bec8c66a46b7346a146e 100644
--- a/extreme_estimator/estimator/margin_estimator.py
+++ b/extreme_estimator/estimator/margin_estimator.py
@@ -43,6 +43,8 @@ class SmoothMarginEstimator(AbstractMarginEstimator):
 
     def _fit(self):
         maxima_gev = self.dataset.maxima_gev(split=self.train_split)
-        coordinate_values = self.dataset.coordinates_values(self.train_split)
-        self._margin_function_fitted = self.margin_model.fitmargin_from_maxima_gev(maxima_gev=maxima_gev,
-                                                                                   coordinates_values=coordinate_values)
+        df_coordinates = self.dataset.df_coordinates(self.train_split)
+        self._params_fitted = self.margin_model.fitmargin_from_maxima_gev(maxima_gev=maxima_gev,
+                                                                          df_coordinates=df_coordinates)
+        self.extract_fitted_models_from_fitted_params(self.margin_model.margin_function_start_fit, self._params_fitted)
+        assert isinstance(self.margin_function_fitted, AbstractMarginFunction)
diff --git a/extreme_estimator/extreme_models/abstract_model.py b/extreme_estimator/extreme_models/abstract_model.py
index 0b80cad354ab2db776095f810f9c626be7bf5c45..016362d0cd73c3bd81957346cca850108e9a0101 100644
--- a/extreme_estimator/extreme_models/abstract_model.py
+++ b/extreme_estimator/extreme_models/abstract_model.py
@@ -1,8 +1,9 @@
 class AbstractModel(object):
 
-    def __init__(self, params_start_fit=None, params_sample=None):
+    def __init__(self, use_start_value=True, params_start_fit=None, params_sample=None):
         self.default_params_start_fit = None
         self.default_params_sample = None
+        self.use_start_value = use_start_value
         self.user_params_start_fit = params_start_fit
         self.user_params_sample = params_sample
 
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 1134907035e7339c990750fd4a140ba5b13a57a3..e604abc6e5019b17235f8d01f3eb1047f7976935 100644
--- a/extreme_estimator/extreme_models/margin_model/abstract_margin_model.py
+++ b/extreme_estimator/extreme_models/margin_model/abstract_margin_model.py
@@ -1,4 +1,5 @@
 import numpy as np
+import pandas as pd
 
 from extreme_estimator.extreme_models.abstract_model import AbstractModel
 from extreme_estimator.extreme_models.margin_model.margin_function.abstract_margin_function \
@@ -10,8 +11,8 @@ from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoo
 
 class AbstractMarginModel(AbstractModel):
 
-    def __init__(self, coordinates: AbstractCoordinates, params_start_fit=None, params_sample=None):
-        super().__init__(params_start_fit, params_sample)
+    def __init__(self, coordinates: AbstractCoordinates, use_start_value=True, params_start_fit=None, params_sample=None):
+        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
@@ -68,7 +69,7 @@ class AbstractMarginModel(AbstractModel):
 
     # Fitting methods needs to be defined in child classes
 
-    def fitmargin_from_maxima_gev(self, maxima_gev: np.ndarray, coordinates_values: np.ndarray) \
+    def fitmargin_from_maxima_gev(self, maxima_gev: np.ndarray, df_coordinates: pd.DataFrame) \
             -> AbstractMarginFunction:
         pass
 
diff --git a/extreme_estimator/extreme_models/margin_model/margin_function/abstract_margin_function.py b/extreme_estimator/extreme_models/margin_model/margin_function/abstract_margin_function.py
index 1072d4a5d4749f32efffa3b6639fce674221ad13..bc593b7ef976f81cf1548693620fdbaec29ff0a5 100644
--- a/extreme_estimator/extreme_models/margin_model/margin_function/abstract_margin_function.py
+++ b/extreme_estimator/extreme_models/margin_model/margin_function/abstract_margin_function.py
@@ -24,6 +24,9 @@ class AbstractMarginFunction(object):
         self.color = 'skyblue'
         self.filter = None
 
+        self._grid_2D = None
+        self._grid_1D = None
+
     def get_gev_params(self, coordinate: np.ndarray) -> GevParams:
         """Main method that maps each coordinate to its GEV parameters"""
         pass
@@ -76,7 +79,7 @@ class AbstractMarginFunction(object):
 
     def visualize_1D(self, gev_value_name=GevParams.GEV_LOC, ax=None, show=True):
         x = self.coordinates.x_coordinates
-        grid, linspace = self.get_grid_values_1D(x)
+        grid, linspace = self.grid_1D(x)
         if ax is None:
             ax = plt.gca()
         if self.datapoint_display:
@@ -91,6 +94,11 @@ class AbstractMarginFunction(object):
         if show:
             plt.show()
 
+    def grid_1D(self, x):
+        if self._grid_1D is None:
+            self._grid_1D = self.get_grid_values_1D(x)
+        return self._grid_1D
+
     def get_grid_values_1D(self, x):
         # TODO: to avoid getting the value several times, I could cache the results
         if self.datapoint_display:
@@ -115,7 +123,7 @@ class AbstractMarginFunction(object):
     def visualize_2D(self, gev_value_name=GevParams.GEV_LOC, ax=None, show=True):
         x = self.coordinates.x_coordinates
         y = self.coordinates.y_coordinates
-        grid = self.get_grid_2D(x, y)
+        grid = self.grid_2D(x, y)
         if ax is None:
             ax = plt.gca()
         imshow_method = ax.imshow
@@ -133,6 +141,11 @@ class AbstractMarginFunction(object):
         if show:
             plt.show()
 
+    def grid_2D(self, x, y):
+        if self._grid_2D is None:
+            self._grid_2D = self.get_grid_2D(x, y)
+        return self._grid_2D
+
     def get_grid_2D(self, x, y):
         resolution = 100
         grid = []
diff --git a/extreme_estimator/extreme_models/margin_model/margin_function/utils.py b/extreme_estimator/extreme_models/margin_model/margin_function/utils.py
index 4f044b987f4cb0778dfaa87d86e05f95949332c4..b367ce6ab688bbf0f30d64957c0285b6fd4630ae 100644
--- a/extreme_estimator/extreme_models/margin_model/margin_function/utils.py
+++ b/extreme_estimator/extreme_models/margin_model/margin_function/utils.py
@@ -22,6 +22,6 @@ def error_dict_between_margin_functions(reference: AbstractMarginFunction, fitte
     fitted_values = fitted.gev_value_name_to_serie
     gev_param_name_to_error_serie = {}
     for gev_value_name in GevParams.GEV_VALUE_NAMES:
-        error = relative_abs_error(reference_values[gev_value_name], fitted_values[gev_value_name])
+        error = 100 * relative_abs_error(reference_values[gev_value_name], fitted_values[gev_value_name])
         gev_param_name_to_error_serie[gev_value_name] = error
     return gev_param_name_to_error_serie
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 ee076c367ca72db60f88bf275fa29a18db9fbedd..1655c1f0550a57cbd87c48e70e54a771ad4d7f33 100644
--- a/extreme_estimator/extreme_models/margin_model/smooth_margin_model.py
+++ b/extreme_estimator/extreme_models/margin_model/smooth_margin_model.py
@@ -1,10 +1,15 @@
+from typing import Dict
+
 import numpy as np
+import pandas as pd
 
 from extreme_estimator.extreme_models.margin_model.margin_function.abstract_margin_function import \
     AbstractMarginFunction
 from extreme_estimator.extreme_models.margin_model.abstract_margin_model import AbstractMarginModel
 from extreme_estimator.extreme_models.margin_model.margin_function.linear_margin_function import LinearMarginFunction
 from extreme_estimator.extreme_models.margin_model.param_function.linear_coef import LinearCoef
+from extreme_estimator.extreme_models.utils import safe_run_r_estimator, r, retrieve_fitted_values, get_coord, \
+    get_margin_formula
 from extreme_estimator.gev_params import GevParams
 
 
@@ -47,10 +52,6 @@ class LinearMarginModel(AbstractMarginModel):
             gev_param_name_to_linear_coef[gev_param_name] = linear_coef
         return gev_param_name_to_linear_coef
 
-    def fitmargin_from_maxima_gev(self, maxima_gev: np.ndarray,
-                                  coordinates_values: np.ndarray) -> AbstractMarginFunction:
-        return self.margin_function_start_fit
-
     @classmethod
     def from_coef_list(cls, coordinates, gev_param_name_to_coef_list):
         params = {}
@@ -59,6 +60,16 @@ class LinearMarginModel(AbstractMarginModel):
                 params[(gev_param_name, dim)] = coef
         return cls(coordinates, params_sample=params, params_start_fit=params)
 
+    def fitmargin_from_maxima_gev(self, maxima_gev: np.ndarray,
+                                  df_coordinates: pd.DataFrame) -> Dict[str, float]:
+        data = np.transpose(maxima_gev)
+        covariables = get_coord(df_coordinates)
+        fit_params = get_margin_formula(self.margin_function_start_fit.form_dict)
+        if self.use_start_value:
+            fit_params['start'] = r.list(**self.margin_function_start_fit.coef_dict)
+        res = safe_run_r_estimator(function=r.fitspatgev, data=data, covariables=covariables, **fit_params)
+        return retrieve_fitted_values(res)
+
 
 class ConstantMarginModel(LinearMarginModel):
 
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 ead7adce0026e9db417b4e30c740d15a0704837c..2c3b4314b0a9743f760c97f7c1a123dd6fd3c33d 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,22 +1,20 @@
-import pandas as pd
 from enum import Enum
 
 import numpy as np
-import rpy2
-from rpy2.rinterface import RRuntimeError, RRuntimeWarning
+import pandas as pd
 import rpy2.robjects as robjects
 
 from extreme_estimator.extreme_models.abstract_model import AbstractModel
-from extreme_estimator.extreme_models.utils import r, set_seed_r
+from extreme_estimator.extreme_models.utils import r, safe_run_r_estimator, retrieve_fitted_values, get_coord, \
+    get_margin_formula
 from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
 
 
 class AbstractMaxStableModel(AbstractModel):
 
-    def __init__(self, params_start_fit=None, params_sample=None):
-        super().__init__(params_start_fit, params_sample)
+    def __init__(self, use_start_value=True, params_start_fit=None, params_sample=None):
+        super().__init__(use_start_value, params_start_fit, params_sample)
         self.cov_mod = None
-        self.start_value_for_fitmaxstab = True
 
     @property
     def cov_mod_param(self):
@@ -48,8 +46,7 @@ class AbstractMaxStableModel(AbstractModel):
             assert AbstractCoordinates.COORDINATE_X in df_coordinates.columns
             df_coordinates[AbstractCoordinates.COORDINATE_Y] = df_coordinates[AbstractCoordinates.COORDINATE_X]
         # Give names to columns to enable a specification of the shape of each marginal parameter
-        coord = robjects.vectors.Matrix(df_coordinates.values)
-        coord.colnames = robjects.StrVector(list(df_coordinates.columns))
+        coord = get_coord(df_coordinates)
 
         #  Prepare the fit_params (a dictionary containing all additional parameters)
         fit_params = self.cov_mod_param.copy()
@@ -58,27 +55,17 @@ class AbstractMaxStableModel(AbstractModel):
         start_dict = self.remove_unused_parameters(start_dict, fitmaxstab_with_one_dimensional_data)
         if fit_marge:
             start_dict.update(margin_start_dict)
-            fit_params.update({k: robjects.Formula(v) for k, v in fit_marge_form_dict.items()})
+            margin_formulas = get_margin_formula(fit_marge_form_dict)
+            fit_params.update(margin_formulas)
         if fitmaxstab_with_one_dimensional_data:
             fit_params['iso'] = True
-        if self.start_value_for_fitmaxstab:
+        if self.use_start_value:
             fit_params['start'] = r.list(**start_dict)
         fit_params['fit.marge'] = fit_marge
 
         # Run the fitmaxstab in R
-        try:
-            res = r.fitmaxstab(data=data, coord=coord, **fit_params)  # type: robjects.ListVector
-        except (RRuntimeError, RRuntimeWarning) as e:
-            if isinstance(e, RRuntimeError):
-                raise Exception('Some R exception have been launched at RunTime: \n {}'.format(e.__repr__()))
-            if isinstance(e, RRuntimeWarning):
-                print(e.__repr__())
-                print('WARNING')
-        # 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
+        res = safe_run_r_estimator(function=r.fitmaxstab, data=data, coord=coord, **fit_params)
+        return retrieve_fitted_values(res)
 
     def rmaxstab(self, nb_obs: int, coordinates_values: np.ndarray) -> np.ndarray:
         """
@@ -102,8 +89,8 @@ class CovarianceFunction(Enum):
 
 class AbstractMaxStableModelWithCovarianceFunction(AbstractMaxStableModel):
 
-    def __init__(self, params_start_fit=None, params_sample=None, covariance_function: CovarianceFunction = None):
-        super().__init__(params_start_fit, params_sample)
+    def __init__(self, use_start_value=True, params_start_fit=None, params_sample=None, covariance_function: CovarianceFunction = None):
+        super().__init__(use_start_value, params_start_fit, params_sample)
         assert covariance_function is not None
         self.covariance_function = covariance_function
         self.default_params_sample = {
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
deleted file mode 100644
index 2217244abc75e58c41cf7f7516d41a8e2fe7d8b5..0000000000000000000000000000000000000000
--- a/extreme_estimator/extreme_models/max_stable_model/max_stable_fit.py
+++ /dev/null
@@ -1,55 +0,0 @@
-# test how to call the max stable method
-import pandas as pd
-
-import numpy as np
-
-from extreme_estimator.extreme_models.utils import 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 = R().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/max_stable_model/max_stable_models.py b/extreme_estimator/extreme_models/max_stable_model/max_stable_models.py
index 347f5ba31495a5d0de0ae1898281a24cabfa5563..0e4c7191ffbc2acd29c0a54e2e1530bf7a2375f4 100644
--- a/extreme_estimator/extreme_models/max_stable_model/max_stable_models.py
+++ b/extreme_estimator/extreme_models/max_stable_model/max_stable_models.py
@@ -6,8 +6,8 @@ from extreme_estimator.extreme_models.max_stable_model.abstract_max_stable_model
 
 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)
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
         self.cov_mod = 'gauss'
         self.default_params_start_fit = {
             'var': 1,
@@ -27,8 +27,8 @@ class Smith(AbstractMaxStableModel):
 
 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)
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
         self.cov_mod = 'brown'
         self.default_params_start_fit = {
             'range': 3,
@@ -42,8 +42,8 @@ class BrownResnick(AbstractMaxStableModel):
 
 class Schlather(AbstractMaxStableModelWithCovarianceFunction):
 
-    def __init__(self, params_start_fit=None, params_sample=None, covariance_function: CovarianceFunction = None):
-        super().__init__(params_start_fit, params_sample, covariance_function)
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
         self.cov_mod = self.covariance_function.name
         self.default_params_sample.update({})
         self.default_params_start_fit = self.default_params_sample.copy()
@@ -51,8 +51,8 @@ class Schlather(AbstractMaxStableModelWithCovarianceFunction):
 
 class Geometric(AbstractMaxStableModelWithCovarianceFunction):
 
-    def __init__(self, params_start_fit=None, params_sample=None, covariance_function: CovarianceFunction = None):
-        super().__init__(params_start_fit, params_sample, covariance_function)
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
         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()
@@ -60,8 +60,8 @@ class Geometric(AbstractMaxStableModelWithCovarianceFunction):
 
 class ExtremalT(AbstractMaxStableModelWithCovarianceFunction):
 
-    def __init__(self, params_start_fit=None, params_sample=None, covariance_function: CovarianceFunction = None):
-        super().__init__(params_start_fit, params_sample, covariance_function)
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
         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()
@@ -69,8 +69,8 @@ class ExtremalT(AbstractMaxStableModelWithCovarianceFunction):
 
 class ISchlather(AbstractMaxStableModelWithCovarianceFunction):
 
-    def __init__(self, params_start_fit=None, params_sample=None, covariance_function: CovarianceFunction = None):
-        super().__init__(params_start_fit, params_sample, covariance_function)
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
         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/extreme_estimator/extreme_models/utils.py b/extreme_estimator/extreme_models/utils.py
index 231544eb9c0bc4d627aa0541a1ca954a692c80e2..0ea8f3283edc5938917a0560727a97d0bd3040a5 100644
--- a/extreme_estimator/extreme_models/utils.py
+++ b/extreme_estimator/extreme_models/utils.py
@@ -1,8 +1,13 @@
 import os.path as op
 import random
 import sys
+from typing import Dict
 
+import pandas as pd
 import rpy2.robjects as ro
+from rpy2 import robjects
+from rpy2.rinterface import RRuntimeWarning
+from rpy2.rinterface._rinterface import RRuntimeError
 
 from rpy2.robjects import numpy2ri
 from rpy2.robjects import pandas2ri
@@ -23,6 +28,36 @@ def get_associated_r_file(python_filepath: str) -> str:
     assert op.exists(r_filepath), r_filepath
     return r_filepath
 
+
+def safe_run_r_estimator(function, **parameters):
+    try:
+        res = function(**parameters)  # type:
+    except (RRuntimeError, RRuntimeWarning) as e:
+        if isinstance(e, RRuntimeError):
+            raise Exception('Some R exception have been launched at RunTime: \n {}'.format(e.__repr__()))
+        if isinstance(e, RRuntimeWarning):
+            print(e.__repr__())
+            print('WARNING')
+    return res
+
+
+def retrieve_fitted_values(res: robjects.ListVector):
+    # 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
+
+
+def get_coord(df_coordinates: pd.DataFrame):
+    coord = robjects.vectors.Matrix(df_coordinates.values)
+    coord.colnames = robjects.StrVector(list(df_coordinates.columns))
+    return coord
+
+
+def get_margin_formula(fit_marge_form_dict) -> Dict:
+    return {k: robjects.Formula(v) for k, v in fit_marge_form_dict.items()}
+
 # def conversion_to_FloatVector(data):
 #     """Convert DataFrame or numpy array into FloatVector for r"""
 #     if isinstance(data, pd.DataFrame):
diff --git a/extreme_estimator/gev/gev_mle_fit.R b/extreme_estimator/gev/fit_gev.R
similarity index 100%
rename from extreme_estimator/gev/gev_mle_fit.R
rename to extreme_estimator/gev/fit_gev.R
diff --git a/extreme_estimator/gev/gev_mle_fit.py b/extreme_estimator/gev/fit_gev.py
similarity index 100%
rename from extreme_estimator/gev/gev_mle_fit.py
rename to extreme_estimator/gev/fit_gev.py
diff --git a/spatio_temporal_dataset/dataset/abstract_dataset.py b/spatio_temporal_dataset/dataset/abstract_dataset.py
index 4564bfa21f0620eb0e1eb8e4416fe07d2be20ae7..be1f4b9711852f11b54dac50304b5e752227d775 100644
--- a/spatio_temporal_dataset/dataset/abstract_dataset.py
+++ b/spatio_temporal_dataset/dataset/abstract_dataset.py
@@ -18,6 +18,11 @@ class AbstractDataset(object):
         self.observations = observations
         self.coordinates = coordinates
 
+    # @property
+    # def max_stable_fitted(self) -> AbstractMarginFunction:
+    #     assert self._margin_function_fitted is not None, 'Error: estimator has not been fitted'
+    #     return self._margin_function_fitted
+
     @property
     def slicer(self):
         return self.coordinates.slicer
@@ -42,9 +47,6 @@ class AbstractDataset(object):
         # todo: maybe I should add the split from the temporal observations
         return self.observations.df_maxima_gev.join(self.coordinates.df_merged)
 
-    def df_coordinates(self, split: Split = Split.all) -> pd.DataFrame:
-        return self.coordinates.df_coordinates(split=split)
-
     # Observation wrapper
 
     def maxima_gev(self, split: Split = Split.all) -> np.ndarray:
@@ -58,6 +60,9 @@ class AbstractDataset(object):
 
     # Coordinates wrapper
 
+    def df_coordinates(self, split: Split = Split.all) -> pd.DataFrame:
+        return self.coordinates.df_coordinates(split=split)
+
     def coordinates_values(self, split: Split = Split.all) -> np.ndarray:
         return self.coordinates.coordinates_values(split=split)
 
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 f50d875da103b6da6f8356c559c139a5a835a68d..e9f1ef207791ce6fc739d16a1ce83b3a3429ba75 100644
--- a/test/test_extreme_estimator/test_estimator/test_margin_estimators.py
+++ b/test/test_extreme_estimator/test_estimator/test_margin_estimators.py
@@ -28,6 +28,7 @@ class TestSmoothMarginEstimator(unittest.TestCase):
                 # Plot
                 if self.DISPLAY:
                     margin_model.margin_function_sample.visualize(show=True)
+                    estimator.margin_function_fitted.visualize(show=True)
         self.assertTrue(True)
 
 
diff --git a/test/test_extreme_estimator/test_extreme_models/__init__.py b/test/test_extreme_estimator/test_extreme_models/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/experiment/regression_margin/visualization_margin_model.py b/test/test_extreme_estimator/test_extreme_models/test_margin_model.py
similarity index 57%
rename from experiment/regression_margin/visualization_margin_model.py
rename to test/test_extreme_estimator/test_extreme_models/test_margin_model.py
index fe5ea8ba9a304c65d166bc0442ecb1ead0378598..c49a3bd355bf23bb624e6089d84c56cdcd324750 100644
--- a/experiment/regression_margin/visualization_margin_model.py
+++ b/test/test_extreme_estimator/test_extreme_models/test_margin_model.py
@@ -8,25 +8,22 @@ from spatio_temporal_dataset.coordinates.spatial_coordinates.coordinates_1D impo
 
 
 class VisualizationMarginModel(unittest.TestCase):
-    DISPLAY = True
-    nb_points = 50
+    DISPLAY = False
+    nb_points = 2
     margin_model = [LinearShapeDim1MarginModel, LinearAllParametersAllDimsMarginModel][-1]
 
-    @classmethod
-    def example_visualization_2D(cls):
-        spatial_coordinates = CircleSpatialCoordinates.from_nb_points(nb_points=cls.nb_points)
-        margin_model = cls.margin_model(coordinates=spatial_coordinates)
-        if cls.DISPLAY:
+    def test_example_visualization_2D(self):
+        spatial_coordinates = CircleSpatialCoordinates.from_nb_points(nb_points=self.nb_points)
+        margin_model = self.margin_model(coordinates=spatial_coordinates)
+        if self.DISPLAY:
             margin_model.margin_function_sample.visualize()
 
-    @classmethod
-    def example_visualization_1D(cls):
-        coordinates = LinSpaceSpatialCoordinates.from_nb_points(nb_points=cls.nb_points)
-        margin_model = cls.margin_model(coordinates=coordinates, params_sample={(GevParams.GEV_SHAPE, 1): 0.02})
-        if cls.DISPLAY:
-            margin_model.margin_function_sample.visualize()
+    def test_example_visualization_1D(self):
+        coordinates = LinSpaceSpatialCoordinates.from_nb_points(nb_points=self.nb_points)
+        margin_model = self.margin_model(coordinates=coordinates, params_sample={(GevParams.GEV_SHAPE, 1): 0.02})
+        margin_model.margin_function_sample.visualize(show=self.DISPLAY)
+        self.assertTrue(True)
 
 
 if __name__ == '__main__':
-    VisualizationMarginModel.example_visualization_1D()
-    # VisualizationMarginModel.example_visualization_2D()
+    unittest.main()
diff --git a/test/test_extreme_estimator/test_gev/test_gev_mle_fit.py b/test/test_extreme_estimator/test_gev/test_gev_mle_fit.py
index 02493e79dd7bf00260b2f640a37486bc9157dfce..9b948f469902f474fbe941d432ec522455a3c7c8 100644
--- a/test/test_extreme_estimator/test_gev/test_gev_mle_fit.py
+++ b/test/test_extreme_estimator/test_gev/test_gev_mle_fit.py
@@ -3,7 +3,7 @@ import unittest
 import numpy as np
 
 from extreme_estimator.extreme_models.utils import r, set_seed_r
-from extreme_estimator.gev.gev_mle_fit import GevMleFit
+from extreme_estimator.gev.fit_gev import GevMleFit
 
 
 class TestGevMleFit(unittest.TestCase):
diff --git a/test/test_unitary/test_fitmaxstab/test_fitmaxstab_with_margin.py b/test/test_unitary/test_fitmaxstab/test_fitmaxstab_with_margin.py
index c0ca47c8eb813101930c9a32a6415b11c6ac3a56..a0921cf48e0da112adfa3aed05925408aa7ec668 100644
--- a/test/test_unitary/test_fitmaxstab/test_fitmaxstab_with_margin.py
+++ b/test/test_unitary/test_fitmaxstab/test_fitmaxstab_with_margin.py
@@ -26,13 +26,12 @@ class TestMaxStableFitWithConstantMargin(TestUnitaryAbstract):
     @property
     def python_output(self):
         dataset = TestRMaxStabWithMarginConstant.python_code()
-        max_stable_model = Schlather(covariance_function=CovarianceFunction.whitmat)
+        max_stable_model = Schlather(covariance_function=CovarianceFunction.whitmat, use_start_value=False)
         margin_model = ConstantMarginModel(dataset.coordinates)
-        max_stable_model.start_value_for_fitmaxstab = False
         full_estimator = FullEstimatorInASingleStepWithSmoothMargin(dataset, margin_model,
                                                                     max_stable_model)
         full_estimator.fit()
-        return full_estimator.full_params_fitted
+        return full_estimator.params_fitted
 
     def test_max_stable_fit_with_constant_margin(self):
         self.compare()
@@ -55,13 +54,12 @@ class TestMaxStableFitWithLinearMargin(TestUnitaryAbstract):
     @property
     def python_output(self):
         dataset = TestRMaxStabWithMarginConstant.python_code()
-        max_stable_model = Schlather(covariance_function=CovarianceFunction.whitmat)
+        max_stable_model = Schlather(covariance_function=CovarianceFunction.whitmat, use_start_value=False)
         margin_model = LinearMarginModelExample(dataset.coordinates)
-        max_stable_model.start_value_for_fitmaxstab = False
         full_estimator = FullEstimatorInASingleStepWithSmoothMargin(dataset, margin_model,
                                                                     max_stable_model)
         full_estimator.fit()
-        return full_estimator.full_params_fitted
+        return full_estimator.params_fitted
 
     def test_max_stable_fit_with_linear_margin(self):
         self.compare()
diff --git a/test/test_unitary/test_fitmaxstab/test_fitmaxstab_without_margin.py b/test/test_unitary/test_fitmaxstab/test_fitmaxstab_without_margin.py
index 9256c0ae32e8cd8ac874fbdc5e2b025535671975..baa4b03c1e1c79a3db2f6067810e37c1116f1f17 100644
--- a/test/test_unitary/test_fitmaxstab/test_fitmaxstab_without_margin.py
+++ b/test/test_unitary/test_fitmaxstab/test_fitmaxstab_without_margin.py
@@ -21,8 +21,7 @@ class TestMaxStableFitWithoutMargin(TestUnitaryAbstract):
     def python_output(self):
         coordinates, max_stable_model = TestRMaxStab.python_code()
         dataset = MaxStableDataset.from_sampling(nb_obs=40, max_stable_model=max_stable_model, coordinates=coordinates)
-        max_stable_model = Schlather(covariance_function=CovarianceFunction.whitmat)
-        max_stable_model.start_value_for_fitmaxstab = False
+        max_stable_model = Schlather(covariance_function=CovarianceFunction.whitmat, use_start_value=False)
         max_stable_estimator = MaxStableEstimator(dataset, max_stable_model)
         max_stable_estimator.fit()
         return max_stable_estimator.max_stable_params_fitted