diff --git a/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/study_visualizer.py b/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/study_visualizer.py
index be39035a8e6bdd82c54f5a65a6f1cf93a9b97839..80016f960b8087cfeee8f91c8172ab30adfe13cc 100644
--- a/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/study_visualizer.py
+++ b/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/study_visualizer.py
@@ -722,7 +722,7 @@ class StudyVisualizer(VisualizationParameters):
         estimator.fit()
 
         # Set visualization attributes for margin_fct
-        margin_fct = estimator.margin_function_fitted
+        margin_fct = estimator.margin_function_from_fit
         margin_fct._visualization_x_limits = self.study.visualization_x_limits
         margin_fct._visualization_y_limits = self.study.visualization_y_limits
         margin_fct.mask_2D = self.study.mask_french_alps
diff --git a/experiment/meteo_france_data/stations_data/comparison_analysis.py b/experiment/meteo_france_data/stations_data/comparison_analysis.py
index dab4fe2bd4a73ad5659cbf6ee130eecb564c01fd..78bdff75f4edb7455cde1f0208fb35366fb500a9 100644
--- a/experiment/meteo_france_data/stations_data/comparison_analysis.py
+++ b/experiment/meteo_france_data/stations_data/comparison_analysis.py
@@ -280,11 +280,11 @@ class ComparisonAnalysis(object):
                                                                        margin_model=margin_model,
                                                                        max_stable_model=max_stable_model)
                 estimator.fit()
-                print(estimator.result_from_fit.margin_coef_dict)
-                print(estimator.result_from_fit.other_coef_dict)
+                print(estimator.result_from_model_fit.margin_coef_dict)
+                print(estimator.result_from_model_fit.other_coef_dict)
                 estimators.append(estimator)
             # Compare the sign of them margin coefficient for the estimators
-            coefs = [{k: v for k, v in e.result_from_fit.margin_coef_dict.items() if 'Coeff1' not in k} for e in
+            coefs = [{k: v for k, v in e.result_from_model_fit.margin_coef_dict.items() if 'Coeff1' not in k} for e in
                      estimators]
             different_sign = [k for k, v in coefs[0].items() if np.sign(coefs[1][k]) != np.sign(v)]
             print('All linear coefficient have the same sign: {}, different_signs for: {}'.format(
diff --git a/experiment/meteo_france_data/stations_data/visualization/comparisons_visualization/comparisons_visualization.py b/experiment/meteo_france_data/stations_data/visualization/comparisons_visualization/comparisons_visualization.py
index 5ee81c4f3b970006e3b8a686dd380e3e4659775e..af0fd9b579eedf432b5b5826a7e3e16618f7ddeb 100644
--- a/experiment/meteo_france_data/stations_data/visualization/comparisons_visualization/comparisons_visualization.py
+++ b/experiment/meteo_france_data/stations_data/visualization/comparisons_visualization/comparisons_visualization.py
@@ -16,13 +16,7 @@ from experiment.meteo_france_data.stations_data.comparison_analysis import Compa
     REANALYSE_STR, ALTITUDE_COLUMN_NAME, STATION_COLUMN_NAME
 from experiment.trend_analysis.univariate_test.abstract_univariate_test import AbstractUnivariateTest
 from experiment.trend_analysis.univariate_test.utils import compute_gev_change_point_test_results
-from extreme_estimator.extreme_models.result_from_fit import ResultFromIsmev
-from extreme_estimator.extreme_models.utils import r, safe_run_r_estimator, ro
 from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
-from utils import classproperty
-
-
-
 
 path = r'/home/erwan/Documents/projects/spatiotemporalextremes/experiment/meteo_france_data/stations_data/csv'
 path_last_csv_file = op.join(path, 'example.csv')
diff --git a/experiment/regression_margin/regression_margin.py b/experiment/regression_margin/regression_margin.py
index 18cc9a577e9dd387c280f2a39500dcd495ba18f4..b16bd28a1980a686b17219158fb30ccda856afb7 100644
--- a/experiment/regression_margin/regression_margin.py
+++ b/experiment/regression_margin/regression_margin.py
@@ -57,7 +57,7 @@ for i in range(nb_estimator):
     margin_model_for_estimator = margin_model_for_estimator_class(coordinates)
     full_estimator = FullEstimatorInASingleStepWithSmoothMargin(dataset, margin_model_for_estimator, max_stable_model)
     full_estimator.fit()
-    full_estimator.margin_function_fitted.visualize_function(axes=axes, show=False)
+    full_estimator.margin_function_from_fit.visualize_function(axes=axes, show=False)
 plt.show()
 
 # Display all the margin on the same graph for comparison
diff --git a/experiment/robustness_plot/estimation_robustness/max_stable_process_plot.py b/experiment/robustness_plot/estimation_robustness/max_stable_process_plot.py
index b98a7b44b733c9c3bd66c6ae25218063ed7e458f..3be6ef87e05a860e01ce8afc54b719dcf77d4689 100644
--- a/experiment/robustness_plot/estimation_robustness/max_stable_process_plot.py
+++ b/experiment/robustness_plot/estimation_robustness/max_stable_process_plot.py
@@ -40,7 +40,11 @@ class MaxStableProcessPlot(object):
                                                  coordinates=spatial_coordinates)
         estimator = MaxStableEstimator(dataset, max_stable_model)
         estimator.fit()
-        return estimator.scalars(max_stable_model.params_sample)
+        # Estimator was computing some error as follows:
+        # absolute_errors = {param_name: np.abs(param_true_value - self.max_stable_params_fitted[param_name])
+        #                    for param_name, param_true_value in true_max_stable_params.items()}
+        # mean_absolute_error = np.mean(np.array(list(absolute_errors.values())))
+        # return estimator.scalars(max_stable_model.params_sample)
 
 
 class SingleMaxStableProcessPlot(SinglePlot, MaxStableProcessPlot):
diff --git a/experiment/robustness_plot/estimation_robustness/spatial_robustness/alps_msp_robustness.py b/experiment/robustness_plot/estimation_robustness/spatial_robustness/alps_msp_robustness.py
index 8e148738f9ab68ccc4a09a00948e41b074b77eb7..e8a7e84b44a0ad7afa420b24f1b168511615c23d 100644
--- a/experiment/robustness_plot/estimation_robustness/spatial_robustness/alps_msp_robustness.py
+++ b/experiment/robustness_plot/estimation_robustness/spatial_robustness/alps_msp_robustness.py
@@ -1,4 +1,3 @@
-from extreme_estimator.estimator.abstract_estimator import AbstractEstimator
 from extreme_estimator.extreme_models.max_stable_model.max_stable_models import Smith, BrownResnick
 from experiment.robustness_plot.estimation_robustness.max_stable_process_plot import MultipleMaxStableProcessPlot, MaxStableProcessPlot
 from experiment.robustness_plot.single_plot import SinglePlot
@@ -45,7 +44,7 @@ def multiple_spatial_robustness_alps():
 
     # Put only the parameter that will vary
     spatial_robustness.robustness_grid_plot(**{
-        SinglePlot.OrdinateItem.name: [AbstractEstimator.MAE_ERROR, AbstractEstimator.DURATION],
+        # SinglePlot.OrdinateItem.name: [AbstractEstimator.MAE_ERROR, AbstractEstimator.DURATION],
         MaxStableProcessPlot.NbStationItem.name: nb_stations,
         MaxStableProcessPlot.NbObservationItem.name: nb_observation,
         MaxStableProcessPlot.MaxStableModelItem.name: msp_models,
diff --git a/experiment/robustness_plot/estimation_robustness/unidimensional_robustness/unidimensional_robustness.py b/experiment/robustness_plot/estimation_robustness/unidimensional_robustness/unidimensional_robustness.py
index 937a65c9c2c7ff26a66f12a11bf378aa25032786..ad55216acee31364f4cea4fe212da48a00e97862 100644
--- a/experiment/robustness_plot/estimation_robustness/unidimensional_robustness/unidimensional_robustness.py
+++ b/experiment/robustness_plot/estimation_robustness/unidimensional_robustness/unidimensional_robustness.py
@@ -1,4 +1,3 @@
-from extreme_estimator.estimator.abstract_estimator import AbstractEstimator
 from extreme_estimator.extreme_models.max_stable_model.max_stable_models import Smith, BrownResnick
 from experiment.robustness_plot.estimation_robustness.max_stable_process_plot import MultipleMaxStableProcessPlot, MaxStableProcessPlot
 from experiment.robustness_plot.single_plot import SinglePlot
diff --git a/experiment/simulation/abstract_simulation.py b/experiment/simulation/abstract_simulation.py
index 22f6ee32c0bd8fbec2320441d06c5cd64529b856..8cd1c03c40ba9d89a999fa35e2c5e59f9066ccea 100644
--- a/experiment/simulation/abstract_simulation.py
+++ b/experiment/simulation/abstract_simulation.py
@@ -59,7 +59,7 @@ class AbstractSimulation(object):
             estimator = estimator_class.from_dataset(dataset)  # type: AbstractEstimator
             # Fit the estimator and get the margin_function
             estimator.fit()
-            margin_function_fitted_list.append(estimator.margin_function_fitted)
+            margin_function_fitted_list.append(estimator.margin_function_from_fit)
 
         # Individual error dict
         self.dump_fitted_margins_pickle(estimator_class, margin_function_fitted_list)
diff --git a/experiment/trend_analysis/non_stationary_trends.py b/experiment/trend_analysis/non_stationary_trends.py
index 0334c3c23b9fb63b9e20861583af48a540bd73c2..6f532ed29ae1509afc6081ccc44792452db350fd 100644
--- a/experiment/trend_analysis/non_stationary_trends.py
+++ b/experiment/trend_analysis/non_stationary_trends.py
@@ -59,7 +59,7 @@ class AbstractNonStationaryTrendTest(object):
             text = 'Fittig {} with margin: {} for starting_point={}\n'.format(estimator_name,
                                                                               margin_model_name,
                                                                               starting_point)
-            text += 'Fit took {}s and was {}'.format(round(duration, 1), estimator.result_from_fit.convergence)
+            text += 'Fit took {}s and was {}'.format(round(duration, 1), estimator.result_from_model_fit.convergence)
             print(text)
         return estimator
 
@@ -68,14 +68,14 @@ class AbstractNonStationaryTrendTest(object):
 
     def get_metric(self, starting_point):
         estimator = self.get_estimator(starting_point)
-        metric = estimator.result_from_fit.__getattribute__(self.RESULT_ATTRIBUTE_METRIC)
+        metric = estimator.result_from_model_fit.__getattribute__(self.RESULT_ATTRIBUTE_METRIC)
         assert isinstance(metric, float)
         return metric
 
     def get_mu_coefs(self, starting_point):
         # for the non stationary model gives the mu1 parameters that was fitted
         estimator = self.get_estimator(starting_point)
-        margin_function = estimator.margin_function_fitted  # type: LinearMarginFunction
+        margin_function = estimator.margin_function_from_fit  # type: LinearMarginFunction
         assert isinstance(margin_function, LinearMarginFunction)
         mu_coefs = [margin_function.mu_intercept, margin_function.mu1_temporal_trend]
         if self.has_spatial_coordinates:
diff --git a/experiment/trend_analysis/univariate_test/abstract_gev_trend_test.py b/experiment/trend_analysis/univariate_test/abstract_gev_trend_test.py
index f9a59b27a53ccfa304bf3adde7edb8250f6c29a6..edea9fdd17260502cb1060c7aa1dcc16801121f6 100644
--- a/experiment/trend_analysis/univariate_test/abstract_gev_trend_test.py
+++ b/experiment/trend_analysis/univariate_test/abstract_gev_trend_test.py
@@ -89,21 +89,21 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
         if self.crashed:
             return np.nan
         else:
-            return self.constrained_estimator.result_from_fit.deviance
+            return self.constrained_estimator.result_from_model_fit.deviance
 
     @property
     def unconstrained_model_deviance(self):
         if self.crashed:
             return np.nan
         else:
-            return self.unconstrained_estimator.result_from_fit.deviance
+            return self.unconstrained_estimator.result_from_model_fit.deviance
 
     @property
     def unconstained_nllh(self):
         if self.crashed:
             return np.nan
         else:
-            return self.unconstrained_estimator.result_from_fit.nllh
+            return self.unconstrained_estimator.result_from_model_fit.nllh
 
     # Evolution of the GEV parameters and corresponding quantiles
 
@@ -112,20 +112,20 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
         return np.sign(self.test_trend_slope_strength)
 
     def get_non_stationary_linear_coef(self, gev_param_name: str):
-        return self.unconstrained_estimator.margin_function_fitted.get_coef(gev_param_name,
-                                                                            AbstractCoordinates.COORDINATE_T)
+        return self.unconstrained_estimator.margin_function_from_fit.get_coef(gev_param_name,
+                                                                              AbstractCoordinates.COORDINATE_T)
 
     @cached_property
     def non_stationary_constant_gev_params(self) -> GevParams:
         # Constant parameters correspond to the gev params in 1958
-        return self.unconstrained_estimator.margin_function_fitted.get_gev_params(coordinate=np.array([1958]),
-                                                                                  is_transformed=False)
+        return self.unconstrained_estimator.margin_function_from_fit.get_gev_params(coordinate=np.array([1958]),
+                                                                                    is_transformed=False)
 
     @cached_property
     def stationary_constant_gev_params(self) -> GevParams:
         # Constant parameters correspond to any gev params
-        return self.constrained_estimator.margin_function_fitted.get_gev_params(coordinate=np.array([1958]),
-                                                                                is_transformed=False)
+        return self.constrained_estimator.margin_function_from_fit.get_gev_params(coordinate=np.array([1958]),
+                                                                                  is_transformed=False)
 
 
     @property
diff --git a/extreme_estimator/estimator/abstract_estimator.py b/extreme_estimator/estimator/abstract_estimator.py
index 2becbbfbcc6780d799f3d0df73bb0f2c16e1d1ec..8d13038b973191e26b53afa1d2a72147faff8909 100644
--- a/extreme_estimator/estimator/abstract_estimator.py
+++ b/extreme_estimator/estimator/abstract_estimator.py
@@ -1,12 +1,10 @@
-import time
+from typing import Union
 
 from cached_property import cached_property
 
-from extreme_estimator.extreme_models.margin_model.linear_margin_model.linear_margin_model import LinearMarginModel
-from extreme_estimator.extreme_models.result_from_fit import ResultFromFit
+from extreme_estimator.extreme_models.result_from_model_fit.abstract_result_from_model_fit import AbstractResultFromModelFit
 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
 
 
@@ -14,34 +12,38 @@ class AbstractEstimator(object):
 
     def __init__(self, dataset: AbstractDataset):
         self.dataset = dataset  # type: AbstractDataset
-        self._result_from_fit = None  # type: ResultFromFit
+        self._result_from_fit = None  # type: Union[None, AbstractResultFromModelFit]
 
+    # Class constructor
     @classmethod
     def from_dataset(cls, dataset: AbstractDataset):
         raise NotImplementedError
 
+    # Fit estimator
+
     def fit(self):
+        self._result_from_fit = self._fit()
+
+    def _fit(self) -> AbstractResultFromModelFit:
         raise NotImplementedError
 
+    # Results from model fit
+
     @property
-    def result_from_fit(self) -> ResultFromFit:
-        assert self._result_from_fit is not None, 'Fit has not be done'
+    def result_from_model_fit(self) -> AbstractResultFromModelFit:
+        assert self._result_from_fit is not None, 'Estimator has not been fitted'
         return self._result_from_fit
 
     @cached_property
-    def margin_function_fitted(self) -> AbstractMarginFunction:
-        return self.extract_function_fitted()
-
-    def extract_function_fitted(self) -> AbstractMarginFunction:
+    def margin_function_from_fit(self) -> AbstractMarginFunction:
         raise NotImplementedError
 
-    def extract_function_fitted_from_the_model_shape(self, margin_model: LinearMarginModel):
-        return LinearMarginFunction.from_coef_dict(coordinates=self.dataset.coordinates,
-                                                   gev_param_name_to_dims=margin_model.margin_function_start_fit.gev_param_name_to_dims,
-                                                   coef_dict=self.result_from_fit.margin_coef_dict,
-                                                   starting_point=margin_model.starting_point)
+    # Short cut properties
 
-    #
     @property
     def train_split(self):
         return self.dataset.train_split
+
+
+
+
diff --git a/extreme_estimator/estimator/full_estimator/abstract_full_estimator.py b/extreme_estimator/estimator/full_estimator/abstract_full_estimator.py
index 85fcb13594b69054d05a6bcef6e54aa8aa9e9fb5..8185ba1072f7003c6586f535644c600535262c4b 100644
--- a/extreme_estimator/estimator/full_estimator/abstract_full_estimator.py
+++ b/extreme_estimator/estimator/full_estimator/abstract_full_estimator.py
@@ -3,6 +3,7 @@ from cached_property import cached_property
 from extreme_estimator.estimator.abstract_estimator import AbstractEstimator
 from extreme_estimator.estimator.margin_estimator.abstract_margin_estimator import LinearMarginEstimator
 from extreme_estimator.estimator.max_stable_estimator.abstract_max_stable_estimator import MaxStableEstimator
+from extreme_estimator.estimator.utils import load_margin_function
 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.linear_margin_model.linear_margin_model import LinearMarginModel
@@ -33,12 +34,18 @@ class SmoothMarginalsThenUnitaryMsp(AbstractFullEstimator):
         maxima_frech = AbstractMarginModel.gev2frech(maxima_gev=maxima_gev_train,
                                                      coordinates_values=self.dataset.coordinates_values(
                                                          self.train_split),
-                                                     margin_function=self.margin_estimator.margin_function_fitted)
+                                                     margin_function=self.margin_estimator.margin_function_from_fit)
         # Update maxima frech field through the dataset object
         self.dataset.set_maxima_frech(maxima_frech, split=self.train_split)
         # Estimate the max stable parameters
         self.max_stable_estimator.fit()
 
+    """
+    To clean things, I could let an abstract Estimator whose only function is fit
+    then create a sub class abstract Estimator One Model, where the function _fit exists and always return something
+    then here full estimator class will be a sub class of abstract estimator, with two  abstract Estimator One Model as attributes
+    """
+
 
 class FullEstimatorInASingleStep(AbstractFullEstimator):
     pass
@@ -67,9 +74,9 @@ class FullEstimatorInASingleStepWithSmoothMargin(AbstractFullEstimator):
         return self.dataset.coordinates.df_temporal_coordinates_for_fit(split=self.train_split,
                                                                         starting_point=self.linear_margin_model.starting_point)
 
-    def fit(self):
+    def _fit(self):
         # Estimate both the margin and the max-stable structure
-        self._result_from_fit = self.max_stable_model.fitmaxstab(
+        return self.max_stable_model.fitmaxstab(
             data_gev=self.dataset.maxima_gev_for_spatial_extremes_package(self.train_split),
             df_coordinates_spat=self.df_coordinates_spat,
             df_coordinates_temp=self.df_coordinates_temp,
@@ -78,12 +85,9 @@ class FullEstimatorInASingleStepWithSmoothMargin(AbstractFullEstimator):
             margin_start_dict=self.linear_margin_model.margin_function_start_fit.coef_dict
         )
 
-    def extract_function_fitted(self):
-        return self.extract_function_fitted_from_the_model_shape(self.linear_margin_model)
-
     @cached_property
-    def margin_function_fitted(self) -> LinearMarginFunction:
-        return super().margin_function_fitted
+    def margin_function_from_fit(self) -> LinearMarginFunction:
+        return load_margin_function(self, self.linear_margin_model)
 
 
 class PointwiseAndThenUnitaryMsp(AbstractFullEstimator):
diff --git a/extreme_estimator/estimator/margin_estimator/abstract_margin_estimator.py b/extreme_estimator/estimator/margin_estimator/abstract_margin_estimator.py
index afdf32c8cae2e684d57975f1503fb6638a640581..be13d011c1939975ac4035b88ca2d21d7fa8794f 100644
--- a/extreme_estimator/estimator/margin_estimator/abstract_margin_estimator.py
+++ b/extreme_estimator/estimator/margin_estimator/abstract_margin_estimator.py
@@ -3,8 +3,10 @@ from abc import ABC
 from cached_property import cached_property
 
 from extreme_estimator.estimator.abstract_estimator import AbstractEstimator
+from extreme_estimator.estimator.utils import load_margin_function
 from extreme_estimator.extreme_models.margin_model.linear_margin_model.linear_margin_model import LinearMarginModel
 from extreme_estimator.extreme_models.margin_model.margin_function.linear_margin_function import LinearMarginFunction
+from extreme_estimator.extreme_models.result_from_model_fit.abstract_result_from_model_fit import AbstractResultFromModelFit
 from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
 
 
@@ -14,29 +16,24 @@ class AbstractMarginEstimator(AbstractEstimator, ABC):
         super().__init__(dataset)
         assert self.dataset.maxima_gev() is not None
 
+
 class LinearMarginEstimator(AbstractMarginEstimator):
     """# with different type of marginals: cosntant, linear...."""
 
-    def _error(self, true_max_stable_params: dict):
-        pass
-
     def __init__(self, dataset: AbstractDataset, margin_model: LinearMarginModel):
         super().__init__(dataset)
         assert isinstance(margin_model, LinearMarginModel)
         self.margin_model = margin_model
 
-    def fit(self):
+    def _fit(self) -> AbstractResultFromModelFit:
         maxima_gev_specialized = self.dataset.maxima_gev_for_spatial_extremes_package(self.train_split)
         df_coordinates_spat = self.dataset.coordinates.df_spatial_coordinates(self.train_split)
         df_coordinates_temp = self.dataset.coordinates.df_temporal_coordinates_for_fit(split=self.train_split,
                                                                                        starting_point=self.margin_model.starting_point)
-        self._result_from_fit = self.margin_model.fitmargin_from_maxima_gev(data=maxima_gev_specialized,
+        return self.margin_model.fitmargin_from_maxima_gev(data=maxima_gev_specialized,
                                                                             df_coordinates_spat=df_coordinates_spat,
                                                                             df_coordinates_temp=df_coordinates_temp)
 
     @cached_property
-    def margin_function_fitted(self) -> LinearMarginFunction:
-        return super().margin_function_fitted
-
-    def extract_function_fitted(self) -> LinearMarginFunction:
-        return self.extract_function_fitted_from_the_model_shape(self.margin_model)
+    def margin_function_from_fit(self) -> LinearMarginFunction:
+        return load_margin_function(self, self.margin_model)
diff --git a/extreme_estimator/estimator/max_stable_estimator/abstract_max_stable_estimator.py b/extreme_estimator/estimator/max_stable_estimator/abstract_max_stable_estimator.py
index 8a9bed2490d1c386769a3a4f6d496247a12b9594..a824b422a35cbe4cfb9b192caf90a672f1171abb 100644
--- a/extreme_estimator/estimator/max_stable_estimator/abstract_max_stable_estimator.py
+++ b/extreme_estimator/estimator/max_stable_estimator/abstract_max_stable_estimator.py
@@ -10,25 +10,19 @@ class AbstractMaxStableEstimator(AbstractEstimator):
     def __init__(self, dataset: AbstractDataset, max_stable_model: AbstractMaxStableModel):
         super().__init__(dataset=dataset)
         self.max_stable_model = max_stable_model
-        # Fit parameters
-        self.max_stable_params_fitted = None
 
+    @property
+    def max_stable_params_fitted(self):
+        raise NotImplementedError
 
 class MaxStableEstimator(AbstractMaxStableEstimator):
 
-    def fit(self):
+    def _fit(self):
         assert self.dataset.maxima_frech(split=self.train_split) is not None
-        self._result_from_fit = self.max_stable_model.fitmaxstab(
+        return self.max_stable_model.fitmaxstab(
             data_frech=self.dataset.maxima_frech_for_spatial_extremes_package(split=self.train_split),
             df_coordinates_spat=self.dataset.df_coordinates(split=self.train_split))
-        self.max_stable_params_fitted = self.result_from_fit.all_parameters
 
-    def scalars(self, true_max_stable_params: dict):
-        error = self._error(true_max_stable_params)
-        return {**error}
-
-    def _error(self, true_max_stable_params: dict):
-        absolute_errors = {param_name: np.abs(param_true_value - self.max_stable_params_fitted[param_name])
-                           for param_name, param_true_value in true_max_stable_params.items()}
-        mean_absolute_error = np.mean(np.array(list(absolute_errors.values())))
-        return {**absolute_errors, **{self.MAE_ERROR: mean_absolute_error}}
+    @property
+    def max_stable_params_fitted(self):
+        return self.result_from_model_fit.all_parameters
diff --git a/extreme_estimator/estimator/utils.py b/extreme_estimator/estimator/utils.py
new file mode 100644
index 0000000000000000000000000000000000000000..f14fd7f4b49410b519a31f2d19f7b6937c134dc9
--- /dev/null
+++ b/extreme_estimator/estimator/utils.py
@@ -0,0 +1,10 @@
+from extreme_estimator.estimator.abstract_estimator import AbstractEstimator
+from extreme_estimator.extreme_models.margin_model.linear_margin_model.linear_margin_model import LinearMarginModel
+from extreme_estimator.extreme_models.margin_model.margin_function.linear_margin_function import LinearMarginFunction
+
+
+def load_margin_function(estimator: AbstractEstimator, margin_model: LinearMarginModel, margin_function_class=LinearMarginFunction):
+    return margin_function_class.from_coef_dict(coordinates=estimator.dataset.coordinates,
+                                                gev_param_name_to_dims=margin_model.margin_function_start_fit.gev_param_name_to_dims,
+                                                coef_dict=estimator.result_from_model_fit.margin_coef_dict,
+                                                starting_point=margin_model.starting_point)
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 7e44d0dc3dce6b50589658b369a3120cd6eb22d4..7db9d3d8b223cb14d3a21908364c1ae22452217b 100644
--- a/extreme_estimator/extreme_models/margin_model/abstract_margin_model.py
+++ b/extreme_estimator/extreme_models/margin_model/abstract_margin_model.py
@@ -6,7 +6,7 @@ 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 \
     import AbstractMarginFunction
-from extreme_estimator.extreme_models.result_from_fit import ResultFromFit
+from extreme_estimator.extreme_models.result_from_model_fit.abstract_result_from_model_fit import AbstractResultFromModelFit
 from extreme_estimator.extreme_models.utils import r
 from extreme_estimator.margin_fits.gev.gev_params import GevParams
 from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
@@ -79,7 +79,7 @@ class AbstractMarginModel(AbstractModel, ABC):
     # Fitting methods needs to be defined in child classes
 
     def fitmargin_from_maxima_gev(self, maxima_gev: np.ndarray, df_coordinates_spatial: pd.DataFrame,
-                                  df_coordinates_temporal: pd.DataFrame) -> ResultFromFit:
+                                  df_coordinates_temporal: pd.DataFrame) -> AbstractResultFromModelFit:
         raise NotImplementedError
 
 
diff --git a/extreme_estimator/extreme_models/margin_model/linear_margin_model/abstract_temporal_linear_margin_model.py b/extreme_estimator/extreme_models/margin_model/linear_margin_model/abstract_temporal_linear_margin_model.py
index 5f3eefffc31aa3f2643ab95787ac90defff69ce1..67c2ffea56df0964b039a125a5fb62214e3676d8 100644
--- a/extreme_estimator/extreme_models/margin_model/linear_margin_model/abstract_temporal_linear_margin_model.py
+++ b/extreme_estimator/extreme_models/margin_model/linear_margin_model/abstract_temporal_linear_margin_model.py
@@ -2,10 +2,11 @@ import numpy as np
 import pandas as pd
 
 from extreme_estimator.extreme_models.margin_model.linear_margin_model.linear_margin_model import LinearMarginModel
-from extreme_estimator.extreme_models.result_from_fit import ResultFromFit, ResultFromIsmev, ResultFromExtremes
+from extreme_estimator.extreme_models.result_from_model_fit.abstract_result_from_model_fit import AbstractResultFromModelFit
+from extreme_estimator.extreme_models.result_from_model_fit.result_from_extremes import ResultFromExtremes
+from extreme_estimator.extreme_models.result_from_model_fit.result_from_ismev import ResultFromIsmev
 from extreme_estimator.extreme_models.utils import r, ro, get_null
 from extreme_estimator.extreme_models.utils import safe_run_r_estimator
-from extreme_estimator.margin_fits.gev.gev_params import GevParams
 from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
 
 
@@ -15,33 +16,47 @@ class AbstractTemporalLinearMarginModel(LinearMarginModel):
     EXTREMES_FEVD_BAYESIAN_FIT_METHOD_STR = 'extRemes.fevd.Bayesian'
 
     def __init__(self, coordinates: AbstractCoordinates, use_start_value=False, params_start_fit=None,
-                 params_sample=None, starting_point=None, fit_method='isMev.gev.fit'):
+                 params_sample=None, starting_point=None, fit_method=ISMEV_GEV_FIT_METHOD_STR):
         super().__init__(coordinates, use_start_value, params_start_fit, params_sample, starting_point)
         assert fit_method in [self.ISMEV_GEV_FIT_METHOD_STR, self.EXTREMES_FEVD_BAYESIAN_FIT_METHOD_STR]
         self.fit_method = fit_method
 
     def fitmargin_from_maxima_gev(self, data: np.ndarray, df_coordinates_spat: pd.DataFrame,
-                                  df_coordinates_temp: pd.DataFrame) -> ResultFromFit:
+                                  df_coordinates_temp: pd.DataFrame) -> AbstractResultFromModelFit:
         assert data.shape[1] == len(df_coordinates_temp.values)
+        x = ro.FloatVector(data[0])
+        y = df_coordinates_temp.values
         if self.fit_method == self.ISMEV_GEV_FIT_METHOD_STR:
-            return self.ismev_gev_fit(data, df_coordinates_temp)
+            return self.ismev_gev_fit(x, y)
         if self.fit_method == self.EXTREMES_FEVD_BAYESIAN_FIT_METHOD_STR:
-            return self.extremes_fevd_bayesian_fit(data, df_coordinates_temp)
+            return self.extremes_fevd_bayesian_fit(x, y)
 
     # Gev Fit with isMev package
 
-    def ismev_gev_fit(self, data: np.ndarray, df_coordinates_temp: pd.DataFrame) -> ResultFromIsmev:
+    def ismev_gev_fit(self, x, y) -> ResultFromIsmev:
         res = safe_run_r_estimator(function=r('gev.fit'), use_start=self.use_start_value,
-                                   xdat=ro.FloatVector(data[0]), y=df_coordinates_temp.values, mul=self.mul,
+                                   xdat=x, y=y, mul=self.mul,
                                    sigl=self.sigl, shl=self.shl)
         return ResultFromIsmev(res, self.margin_function_start_fit.gev_param_name_to_dims)
 
     # Gev fit with extRemes package
 
-    def extremes_fevd_bayesian_fit(self, data, df_coordinates_temp) -> ResultFromExtremes:
-        res = safe_run_r_estimator(function=r('fevd_fixed'), use_start=self.use_start_value,
-                                   xdat=ro.FloatVector(data[0]), y=df_coordinates_temp.values, mul=self.mul,
-                                   sigl=self.sigl, shl=self.shl)
+    def extremes_fevd_bayesian_fit(self, x, y) -> ResultFromExtremes:
+        # (x_gev, method='Bayesian', priorFun="fevdPriorMyMy", priorParams=list(q=c(6), p = c(
+        #     9)), iter = 5000, verbose = TRUE, use.phi = FALSE)
+
+        r_type_argument_kwargs = {'use.phi': False}
+        # location.fun = ~1,
+        # scale.fun = ~1, shape.fun = ~1
+        # , priorParams = list(q=c(6), p=c(9))
+        res = safe_run_r_estimator(function=r('fevd_fixed'),
+                                   x=x,
+                                   data=y,
+                                   method='Bayesian',
+                                   # priorFun="fevdPriorCustom",
+                                   iter=5000,
+                                   **r_type_argument_kwargs
+                                   )
         return ResultFromExtremes(res, self.margin_function_start_fit.gev_param_name_to_dims)
 
     # Default arguments for all methods
@@ -61,8 +76,3 @@ class AbstractTemporalLinearMarginModel(LinearMarginModel):
     @property
     def siglink(self):
         return r('identity')
-
-
-
-
-
diff --git a/extreme_estimator/extreme_models/margin_model/parametric_margin_model.py b/extreme_estimator/extreme_models/margin_model/parametric_margin_model.py
index e9241c4586cf96df592a4e30a41657ad5cd9ddab..be503f032ec11e4f467014660ea5d30f61dc6d25 100644
--- a/extreme_estimator/extreme_models/margin_model/parametric_margin_model.py
+++ b/extreme_estimator/extreme_models/margin_model/parametric_margin_model.py
@@ -5,15 +5,11 @@ import pandas as pd
 
 from extreme_estimator.extreme_models.margin_model.margin_function.parametric_margin_function import \
     ParametricMarginFunction
-from extreme_estimator.extreme_models.result_from_fit import ResultFromFit, ResultFromSpatialExtreme
+from extreme_estimator.extreme_models.result_from_model_fit.result_from_spatial_extreme import ResultFromSpatialExtreme
 from extreme_estimator.extreme_models.margin_model.abstract_margin_model import AbstractMarginModel
 from extreme_estimator.extreme_models.utils import safe_run_r_estimator, r, get_coord, \
     get_margin_formula
 from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
-from spatio_temporal_dataset.coordinates.spatio_temporal_coordinates.abstract_spatio_temporal_coordinates import \
-    AbstractSpatioTemporalCoordinates
-from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_coordinates import \
-    AbstractTemporalCoordinates
 
 
 class ParametricMarginModel(AbstractMarginModel, ABC):
@@ -29,7 +25,7 @@ class ParametricMarginModel(AbstractMarginModel, ABC):
         super().__init__(coordinates, use_start_value, params_start_fit, params_sample)
 
     def fitmargin_from_maxima_gev(self, data: np.ndarray, df_coordinates_spat: pd.DataFrame,
-                                  df_coordinates_temp: pd.DataFrame) -> ResultFromFit:
+                                  df_coordinates_temp: pd.DataFrame) -> ResultFromSpatialExtreme:
         assert data.shape[1] == len(df_coordinates_spat)
 
         # Margin formula for fitspatgev
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 57cd069e2974791b8bc26b9eb5ed2dabb29fbb22..b1fc374d860c02b881108ef3ed5a302881d3b8ba 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
@@ -6,7 +6,7 @@ from rpy2.rinterface import RRuntimeWarning
 from rpy2.rinterface._rinterface import RRuntimeError
 
 from extreme_estimator.extreme_models.abstract_model import AbstractModel
-from extreme_estimator.extreme_models.result_from_fit import ResultFromFit, ResultFromSpatialExtreme
+from extreme_estimator.extreme_models.result_from_model_fit.result_from_spatial_extreme import ResultFromSpatialExtreme
 from extreme_estimator.extreme_models.utils import r, safe_run_r_estimator, get_coord, \
     get_margin_formula, SafeRunException
 from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
@@ -24,7 +24,7 @@ class AbstractMaxStableModel(AbstractModel):
 
     def fitmaxstab(self, df_coordinates_spat: pd.DataFrame, df_coordinates_temp: pd.DataFrame = None,
                    data_frech: np.ndarray = None, data_gev: np.ndarray = None,
-                   fit_marge=False, fit_marge_form_dict=None, margin_start_dict=None) -> ResultFromFit:
+                   fit_marge=False, fit_marge_form_dict=None, margin_start_dict=None) -> ResultFromSpatialExtreme:
         assert isinstance(df_coordinates_spat, pd.DataFrame)
         if fit_marge:
             assert fit_marge_form_dict is not None
diff --git a/extreme_estimator/extreme_models/result_from_fit.py b/extreme_estimator/extreme_models/result_from_fit.py
deleted file mode 100644
index 212ce05f883090f199461ec99b18bbaff08f856a..0000000000000000000000000000000000000000
--- a/extreme_estimator/extreme_models/result_from_fit.py
+++ /dev/null
@@ -1,168 +0,0 @@
-from typing import Dict
-
-import numpy as np
-from rpy2 import robjects
-
-from extreme_estimator.extreme_models.margin_model.param_function.linear_coef import LinearCoef
-from extreme_estimator.margin_fits.gev.gev_params import GevParams
-from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
-
-
-def convertFloatVector_to_float(f):
-    return np.array(f)[0]
-
-
-class ResultFromFit(object):
-
-    def __init__(self, result_from_fit: robjects.ListVector) -> None:
-        if hasattr(result_from_fit, 'names'):
-            self.name_to_value = {name: result_from_fit.rx2(name) for name in result_from_fit.names}
-        else:
-            self.name_to_value = {}
-
-    @property
-    def names(self):
-        return self.name_to_value.keys()
-
-    @property
-    def all_parameters(self):
-        raise NotImplementedError
-
-    @property
-    def margin_coef_dict(self):
-        raise NotImplementedError
-
-    @property
-    def other_coef_dict(self):
-        raise NotImplementedError
-
-    @property
-    def nllh(self):
-        raise NotImplementedError
-
-    @property
-    def deviance(self):
-        raise NotImplementedError
-
-    @property
-    def convergence(self) -> str:
-        raise NotImplementedError
-
-
-class ResultFromIsmev(ResultFromFit):
-
-    def __init__(self, result_from_fit: robjects.ListVector, gev_param_name_to_dim=None) -> None:
-        super().__init__(result_from_fit)
-        self.gev_param_name_to_dim = gev_param_name_to_dim
-
-    @property
-    def margin_coef_dict(self):
-        assert self.gev_param_name_to_dim is not None
-        # Build the Coeff dict from gev_param_name_to_dim
-        coef_dict = {}
-        i = 0
-        mle_values = self.name_to_value['mle']
-        for gev_param_name in GevParams.PARAM_NAMES:
-            # Add intercept
-            intercept_coef_name = LinearCoef.coef_template_str(gev_param_name, LinearCoef.INTERCEPT_NAME).format(1)
-            coef_dict[intercept_coef_name] = mle_values[i]
-            i += 1
-            # Add a potential linear temporal trend
-            if gev_param_name in self.gev_param_name_to_dim:
-                temporal_coef_name = LinearCoef.coef_template_str(gev_param_name,
-                                                                  AbstractCoordinates.COORDINATE_T).format(1)
-                coef_dict[temporal_coef_name] = mle_values[i]
-                i += 1
-        return coef_dict
-
-    @property
-    def all_parameters(self):
-        return self.margin_coef_dict
-
-    @property
-    def nllh(self):
-        return convertFloatVector_to_float(self.name_to_value['nllh'])
-
-    @property
-    def deviance(self):
-        return - 2 * self.nllh
-
-    @property
-    def convergence(self) -> str:
-        return convertFloatVector_to_float(self.name_to_value['conv']) == 0
-
-
-class ResultFromExtremes(ResultFromFit):
-
-    def __init__(self, result_from_fit: robjects.ListVector, gev_param_name_to_dim=None) -> None:
-        super().__init__(result_from_fit)
-        self.gev_param_name_to_dim = gev_param_name_to_dim
-
-    @property
-    def margin_coef_dict(self):
-        assert self.gev_param_name_to_dim is not None
-        # Build the Coeff dict from gev_param_name_to_dim
-        coef_dict = {}
-        i = 0
-        mle_values = self.name_to_value['mle']
-        for gev_param_name in GevParams.PARAM_NAMES:
-            # Add intercept
-            intercept_coef_name = LinearCoef.coef_template_str(gev_param_name, LinearCoef.INTERCEPT_NAME).format(1)
-            coef_dict[intercept_coef_name] = mle_values[i]
-            i += 1
-            # Add a potential linear temporal trend
-            if gev_param_name in self.gev_param_name_to_dim:
-                temporal_coef_name = LinearCoef.coef_template_str(gev_param_name,
-                                                                  AbstractCoordinates.COORDINATE_T).format(1)
-                coef_dict[temporal_coef_name] = mle_values[i]
-                i += 1
-        return coef_dict
-
-    @property
-    def all_parameters(self):
-        return self.margin_coef_dict
-
-    @property
-    def nllh(self):
-        return convertFloatVector_to_float(self.name_to_value['nllh'])
-
-    @property
-    def deviance(self):
-        return - 2 * self.nllh
-
-    @property
-    def convergence(self) -> str:
-        return convertFloatVector_to_float(self.name_to_value['conv']) == 0
-
-class ResultFromSpatialExtreme(ResultFromFit):
-    """
-    Handler from any result with the result of a fit functions from the package Spatial Extreme
-    """
-    FITTED_VALUES_NAME = 'fitted.values'
-    CONVERGENCE_NAME = 'convergence'
-
-    @property
-    def deviance(self):
-        return np.array(self.name_to_value['deviance'])[0]
-
-    @property
-    def convergence(self) -> str:
-        convergence_value = self.name_to_value[self.CONVERGENCE_NAME]
-        return convergence_value[0]
-
-    @property
-    def is_convergence_successful(self) -> bool:
-        return self.convergence == "successful"
-
-    @property
-    def all_parameters(self) -> Dict[str, float]:
-        fitted_values = self.name_to_value[self.FITTED_VALUES_NAME]
-        return {key: fitted_values.rx2(key)[0] for key in fitted_values.names}
-
-    @property
-    def margin_coef_dict(self):
-        return {k: v for k, v in self.all_parameters.items() if LinearCoef.COEFF_STR in k}
-
-    @property
-    def other_coef_dict(self):
-        return {k: v for k, v in self.all_parameters.items() if LinearCoef.COEFF_STR not in k}
diff --git a/extreme_estimator/extreme_models/result_from_model_fit/__init__.py b/extreme_estimator/extreme_models/result_from_model_fit/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/extreme_estimator/extreme_models/result_from_model_fit/abstract_result_from_model_fit.py b/extreme_estimator/extreme_models/result_from_model_fit/abstract_result_from_model_fit.py
new file mode 100644
index 0000000000000000000000000000000000000000..5560e1c862877450498c1f9807fb03cc71465456
--- /dev/null
+++ b/extreme_estimator/extreme_models/result_from_model_fit/abstract_result_from_model_fit.py
@@ -0,0 +1,45 @@
+import numpy as np
+from rpy2 import robjects
+
+
+
+class AbstractResultFromModelFit(object):
+
+    def __init__(self, result_from_fit: robjects.ListVector) -> None:
+        if hasattr(result_from_fit, 'names'):
+            self.name_to_value = {name: result_from_fit.rx2(name) for name in result_from_fit.names}
+        else:
+            self.name_to_value = {}
+
+    @property
+    def names(self):
+        return self.name_to_value.keys()
+
+    @property
+    def all_parameters(self):
+        raise NotImplementedError
+
+    @property
+    def margin_coef_dict(self):
+        raise NotImplementedError
+
+    @property
+    def other_coef_dict(self):
+        raise NotImplementedError
+
+    @property
+    def nllh(self):
+        raise NotImplementedError
+
+    @property
+    def deviance(self):
+        raise NotImplementedError
+
+    @property
+    def convergence(self) -> str:
+        raise NotImplementedError
+
+
+
+
+
diff --git a/extreme_estimator/extreme_models/result_from_model_fit/result_from_extremes.py b/extreme_estimator/extreme_models/result_from_model_fit/result_from_extremes.py
new file mode 100644
index 0000000000000000000000000000000000000000..0103b8bec3e822dc8b8c447fddcd8c289cd61fa6
--- /dev/null
+++ b/extreme_estimator/extreme_models/result_from_model_fit/result_from_extremes.py
@@ -0,0 +1,52 @@
+from rpy2 import robjects
+
+from extreme_estimator.extreme_models.result_from_model_fit.abstract_result_from_model_fit import \
+    AbstractResultFromModelFit
+
+
+class ResultFromExtremes(AbstractResultFromModelFit):
+
+    def __init__(self, result_from_fit: robjects.ListVector, gev_param_name_to_dim=None) -> None:
+        super().__init__(result_from_fit)
+        self.gev_param_name_to_dim = gev_param_name_to_dim
+        print(list(self.name_to_value.keys()))
+
+    # @property
+    # def
+
+    # @property
+    # def margin_coef_dict(self):
+    #     assert self.gev_param_name_to_dim is not None
+    #     # Build the Coeff dict from gev_param_name_to_dim
+    #     coef_dict = {}
+    #     i = 0
+    #     mle_values = self.name_to_value['mle']
+    #     for gev_param_name in GevParams.PARAM_NAMES:
+    #         # Add intercept
+    #         intercept_coef_name = LinearCoef.coef_template_str(gev_param_name, LinearCoef.INTERCEPT_NAME).format(1)
+    #         coef_dict[intercept_coef_name] = mle_values[i]
+    #         i += 1
+    #         # Add a potential linear temporal trend
+    #         if gev_param_name in self.gev_param_name_to_dim:
+    #             temporal_coef_name = LinearCoef.coef_template_str(gev_param_name,
+    #                                                               AbstractCoordinates.COORDINATE_T).format(1)
+    #             coef_dict[temporal_coef_name] = mle_values[i]
+    #             i += 1
+    #     return coef_dict
+
+    # @property
+    # def all_parameters(self):
+    #     return self.margin_coef_dict
+    #
+    # @property
+    # def nllh(self):
+    #     return convertFloatVector_to_float(self.name_to_value['nllh'])
+    #
+    # @property
+    # def deviance(self):
+    #     return - 2 * self.nllh
+    #
+    # @property
+    # def convergence(self) -> str:
+    #     return convertFloatVector_to_float(self.name_to_value['conv']) == 0
+
diff --git a/extreme_estimator/extreme_models/result_from_model_fit/result_from_ismev.py b/extreme_estimator/extreme_models/result_from_model_fit/result_from_ismev.py
new file mode 100644
index 0000000000000000000000000000000000000000..0300d86ecaa5eff51e90d59b28982ced053402e3
--- /dev/null
+++ b/extreme_estimator/extreme_models/result_from_model_fit/result_from_ismev.py
@@ -0,0 +1,54 @@
+from typing import Dict
+
+import numpy as np
+from rpy2 import robjects
+
+from extreme_estimator.extreme_models.margin_model.param_function.linear_coef import LinearCoef
+from extreme_estimator.extreme_models.result_from_model_fit.abstract_result_from_model_fit import \
+    AbstractResultFromModelFit
+from extreme_estimator.extreme_models.result_from_model_fit.utils import convertFloatVector_to_float
+from extreme_estimator.margin_fits.gev.gev_params import GevParams
+from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
+
+
+class ResultFromIsmev(AbstractResultFromModelFit):
+
+    def __init__(self, result_from_fit: robjects.ListVector, gev_param_name_to_dim=None) -> None:
+        super().__init__(result_from_fit)
+        self.gev_param_name_to_dim = gev_param_name_to_dim
+
+    @property
+    def margin_coef_dict(self):
+        assert self.gev_param_name_to_dim is not None
+        # Build the Coeff dict from gev_param_name_to_dim
+        coef_dict = {}
+        i = 0
+        mle_values = self.name_to_value['mle']
+        for gev_param_name in GevParams.PARAM_NAMES:
+            # Add intercept
+            intercept_coef_name = LinearCoef.coef_template_str(gev_param_name, LinearCoef.INTERCEPT_NAME).format(1)
+            coef_dict[intercept_coef_name] = mle_values[i]
+            i += 1
+            # Add a potential linear temporal trend
+            if gev_param_name in self.gev_param_name_to_dim:
+                temporal_coef_name = LinearCoef.coef_template_str(gev_param_name,
+                                                                  AbstractCoordinates.COORDINATE_T).format(1)
+                coef_dict[temporal_coef_name] = mle_values[i]
+                i += 1
+        return coef_dict
+
+    @property
+    def all_parameters(self):
+        return self.margin_coef_dict
+
+    @property
+    def nllh(self):
+        return convertFloatVector_to_float(self.name_to_value['nllh'])
+
+    @property
+    def deviance(self):
+        return - 2 * self.nllh
+
+    @property
+    def convergence(self) -> str:
+        return convertFloatVector_to_float(self.name_to_value['conv']) == 0
diff --git a/extreme_estimator/extreme_models/result_from_model_fit/result_from_spatial_extreme.py b/extreme_estimator/extreme_models/result_from_model_fit/result_from_spatial_extreme.py
new file mode 100644
index 0000000000000000000000000000000000000000..bf98c66558dca817ebb3586df7a597bbc2fb83b2
--- /dev/null
+++ b/extreme_estimator/extreme_models/result_from_model_fit/result_from_spatial_extreme.py
@@ -0,0 +1,41 @@
+from typing import Dict
+
+import numpy as np
+
+from extreme_estimator.extreme_models.margin_model.param_function.linear_coef import LinearCoef
+from extreme_estimator.extreme_models.result_from_model_fit.abstract_result_from_model_fit import \
+    AbstractResultFromModelFit
+
+
+class ResultFromSpatialExtreme(AbstractResultFromModelFit):
+    """
+    Handler from any result with the result of a fit functions from the package Spatial Extreme
+    """
+    FITTED_VALUES_NAME = 'fitted.values'
+    CONVERGENCE_NAME = 'convergence'
+
+    @property
+    def deviance(self):
+        return np.array(self.name_to_value['deviance'])[0]
+
+    @property
+    def convergence(self) -> str:
+        convergence_value = self.name_to_value[self.CONVERGENCE_NAME]
+        return convergence_value[0]
+
+    @property
+    def is_convergence_successful(self) -> bool:
+        return self.convergence == "successful"
+
+    @property
+    def all_parameters(self) -> Dict[str, float]:
+        fitted_values = self.name_to_value[self.FITTED_VALUES_NAME]
+        return {key: fitted_values.rx2(key)[0] for key in fitted_values.names}
+
+    @property
+    def margin_coef_dict(self):
+        return {k: v for k, v in self.all_parameters.items() if LinearCoef.COEFF_STR in k}
+
+    @property
+    def other_coef_dict(self):
+        return {k: v for k, v in self.all_parameters.items() if LinearCoef.COEFF_STR not in k}
diff --git a/extreme_estimator/extreme_models/result_from_model_fit/utils.py b/extreme_estimator/extreme_models/result_from_model_fit/utils.py
new file mode 100644
index 0000000000000000000000000000000000000000..ebd20e933d72386e9f91564710cffbe02f379a78
--- /dev/null
+++ b/extreme_estimator/extreme_models/result_from_model_fit/utils.py
@@ -0,0 +1,5 @@
+import numpy as np
+
+
+def convertFloatVector_to_float(f):
+    return np.array(f)[0]
diff --git a/extreme_estimator/extreme_models/utils.py b/extreme_estimator/extreme_models/utils.py
index c2fe6ea0c62534f09001e462712ce4fc2560d955..d085837c6c8ff3590cd2ed8c5fb841f628224088 100644
--- a/extreme_estimator/extreme_models/utils.py
+++ b/extreme_estimator/extreme_models/utils.py
@@ -18,13 +18,23 @@ from rpy2.rinterface._rinterface import RRuntimeError
 from rpy2.robjects import numpy2ri
 from rpy2.robjects import pandas2ri
 
+from utils import get_root_path
+
+# Load R variables
 r = ro.R()
 numpy2ri.activate()
 pandas2ri.activate()
 r.library('SpatialExtremes')
+# Desactivate temporarily warnings
 default_filters = warnings.filters.copy()
 warnings.filterwarnings("ignore")
+# Load ismev
 r.library('ismev')
+# Load fevd fixed
+fevd_fixed_filepath = op.join(get_root_path(), 'extreme_estimator', 'margin_fits', 'gev', 'fevd_fixed.R')
+assert op.exists(fevd_fixed_filepath)
+r.source(fevd_fixed_filepath)
+# Reactivate warning
 warnings.filters = default_filters
 
 
diff --git a/extreme_estimator/margin_fits/margin_fits_utils.py b/extreme_estimator/margin_fits/margin_fits_utils.py
index bc41ff956447f852acf85242a95958a7b49f0a05..182ff5570fc329fe6941515a585b2ecb9851be5a 100644
--- a/extreme_estimator/margin_fits/margin_fits_utils.py
+++ b/extreme_estimator/margin_fits/margin_fits_utils.py
@@ -1,8 +1,6 @@
-import numpy as np
 import rpy2.robjects as ro
 
-from extreme_estimator.extreme_models.result_from_fit import ResultFromIsmev
-from extreme_estimator.extreme_models.utils import r, get_null, safe_run_r_estimator
+from extreme_estimator.extreme_models.utils import r, get_null
 
 """
 These two functions are “extremely light” functions to fit the GEV/GPD. These functions are mainlyuseful 
diff --git a/test/test_experiment/test_coordinate_sensitivity.py b/test/test_experiment/test_coordinate_sensitivity.py
index f6178e13beeb219dc31c41cc1f3a784b32439968..081d5ce00fc46ec9a49152c2ce79237f5b9b1452 100644
--- a/test/test_experiment/test_coordinate_sensitivity.py
+++ b/test/test_experiment/test_coordinate_sensitivity.py
@@ -33,15 +33,15 @@ class TestCoordinateSensitivity(unittest.TestCase):
                 if self.DISPLAY:
                     print('Stationary')
                     stationary_est = trend_test.get_estimator(starting_point=None)
-                    print(stationary_est.margin_function_fitted.coordinates.df_all_coordinates)
-                    print(stationary_est.result_from_fit.convergence)
-                    print(stationary_est.margin_function_fitted.coef_dict)
+                    print(stationary_est.margin_function_from_fit.coordinates.df_all_coordinates)
+                    print(stationary_est.result_from_model_fit.convergence)
+                    print(stationary_est.margin_function_from_fit.coef_dict)
                     print('Non Stationary')
                     non_stationary_est = trend_test.get_estimator(starting_point=1960)
-                    print(non_stationary_est.result_from_fit.convergence)
+                    print(non_stationary_est.result_from_model_fit.convergence)
                     non_stationary_est = trend_test.get_estimator(starting_point=1990)
-                    print(non_stationary_est.result_from_fit.convergence)
-                    print(non_stationary_est.margin_function_fitted.coef_dict)
+                    print(non_stationary_est.result_from_model_fit.convergence)
+                    print(non_stationary_est.margin_function_from_fit.coef_dict)
                     print(get_display_name_from_object_type(transformation_class), 'mu1s: ', mu1s)
                     print('\n')
                 self.assertTrue(0.0 not in mu1s)
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 11f2813ff3676d53e4f4f1f7e9a9cb3974210057..7d6ac6a66ca3b4cbee533f9904fb230b5cea20c2 100644
--- a/test/test_extreme_estimator/test_estimator/test_margin_estimators.py
+++ b/test/test_extreme_estimator/test_estimator/test_margin_estimators.py
@@ -30,11 +30,11 @@ class TestSmoothMarginEstimator(unittest.TestCase):
                 # Fit estimator
                 estimator = LinearMarginEstimator(dataset=dataset, margin_model=margin_model)
                 estimator.fit()
-                print(estimator.result_from_fit.name_to_value.keys())
+                print(estimator.result_from_model_fit.name_to_value.keys())
                 # Plot
                 if self.DISPLAY:
                     margin_model.margin_function_sample.visualize_function(show=True)
-                    estimator.margin_function_fitted.visualize_function(show=True)
+                    estimator.margin_function_from_fit.visualize_function(show=True)
         self.assertTrue(True)
 
 
diff --git a/test/test_extreme_estimator/test_extreme_models/test_margin_temporal.py b/test/test_extreme_estimator/test_extreme_models/test_margin_temporal.py
index affe7ad8e7241f9432a972df7ce092ecddce4fcc..7e2c5f6e08f53beb30c1fb6393f25f99b921aaef 100644
--- a/test/test_extreme_estimator/test_extreme_models/test_margin_temporal.py
+++ b/test/test_extreme_estimator/test_extreme_models/test_margin_temporal.py
@@ -33,7 +33,7 @@ class TestMarginTemporal(unittest.TestCase):
         ref = {'loc': 1.3456595684773085, 'scale': 1.090369430386199, 'shape': 0.6845422250749476}
         for year in range(1, 3):
             coordinate = np.array([0.0, 0.0, year])
-            mle_params_estimated = estimator.margin_function_fitted.get_gev_params(coordinate).to_dict()
+            mle_params_estimated = estimator.margin_function_from_fit.get_gev_params(coordinate).to_dict()
             for key in ref.keys():
                 self.assertAlmostEqual(ref[key], mle_params_estimated[key], places=3)
 
@@ -42,32 +42,32 @@ class TestMarginTemporal(unittest.TestCase):
         margin_model = LinearNonStationaryLocationMarginModel(self.coordinates)
         estimator = LinearMarginEstimator(self.dataset, margin_model)
         estimator.fit()
-        self.assertNotEqual(estimator.margin_function_fitted.mu1_temporal_trend, 0.0)
+        self.assertNotEqual(estimator.margin_function_from_fit.mu1_temporal_trend, 0.0)
         # Checks that parameters returned are indeed different
         coordinate1 = np.array([0.0, 0.0, 1])
-        mle_params_estimated_year1 = estimator.margin_function_fitted.get_gev_params(coordinate1).to_dict()
+        mle_params_estimated_year1 = estimator.margin_function_from_fit.get_gev_params(coordinate1).to_dict()
         coordinate3 = np.array([0.0, 0.0, 3])
-        mle_params_estimated_year3 = estimator.margin_function_fitted.get_gev_params(coordinate3).to_dict()
+        mle_params_estimated_year3 = estimator.margin_function_from_fit.get_gev_params(coordinate3).to_dict()
         self.assertNotEqual(mle_params_estimated_year1, mle_params_estimated_year3)
 
     def test_margin_fit_nonstationary_with_start_point(self):
         # Create estimator
         estimator = self.fit_non_stationary_estimator(starting_point=2)
         # By default, estimator find the good margin
-        self.assertNotEqual(estimator.margin_function_fitted.mu1_temporal_trend, 0.0)
-        self.assertAlmostEqual(estimator.margin_function_fitted.mu1_temporal_trend,
+        self.assertNotEqual(estimator.margin_function_from_fit.mu1_temporal_trend, 0.0)
+        self.assertAlmostEqual(estimator.margin_function_from_fit.mu1_temporal_trend,
                                self.smooth_margin_model.margin_function_sample.mu1_temporal_trend,
                                places=3)
         # Checks starting point parameter are well passed
-        self.assertEqual(2, estimator.margin_function_fitted.starting_point)
+        self.assertEqual(2, estimator.margin_function_from_fit.starting_point)
         # Checks that parameters returned are indeed different
         coordinate1 = np.array([0.0, 0.0, 1])
-        mle_params_estimated_year1 = estimator.margin_function_fitted.get_gev_params(coordinate1).to_dict()
+        mle_params_estimated_year1 = estimator.margin_function_from_fit.get_gev_params(coordinate1).to_dict()
         coordinate2 = np.array([0.0, 0.0, 2])
-        mle_params_estimated_year2 = estimator.margin_function_fitted.get_gev_params(coordinate2).to_dict()
+        mle_params_estimated_year2 = estimator.margin_function_from_fit.get_gev_params(coordinate2).to_dict()
         self.assertEqual(mle_params_estimated_year1, mle_params_estimated_year2)
         coordinate5 = np.array([0.0, 0.0, 5])
-        mle_params_estimated_year5 = estimator.margin_function_fitted.get_gev_params(coordinate5).to_dict()
+        mle_params_estimated_year5 = estimator.margin_function_from_fit.get_gev_params(coordinate5).to_dict()
         self.assertNotEqual(mle_params_estimated_year5, mle_params_estimated_year2)
 
     def fit_non_stationary_estimator(self, starting_point):
@@ -80,8 +80,8 @@ class TestMarginTemporal(unittest.TestCase):
         # Create two different estimators
         estimator1 = self.fit_non_stationary_estimator(starting_point=3)
         estimator2 = self.fit_non_stationary_estimator(starting_point=20)
-        mu1_estimator1 = estimator1.margin_function_fitted.mu1_temporal_trend
-        mu1_estimator2 = estimator2.margin_function_fitted.mu1_temporal_trend
+        mu1_estimator1 = estimator1.margin_function_from_fit.mu1_temporal_trend
+        mu1_estimator2 = estimator2.margin_function_from_fit.mu1_temporal_trend
         self.assertNotEqual(mu1_estimator1, mu1_estimator2)
 
 
diff --git a/test/test_extreme_estimator/test_extreme_models/test_margin_temporal_transformed.py b/test/test_extreme_estimator/test_extreme_models/test_margin_temporal_transformed.py
index d3a9ca087df1121d101eda3065b802af6ef053ef..975c7f288d31a1492be3975c9e352ee08e63899b 100644
--- a/test/test_extreme_estimator/test_extreme_models/test_margin_temporal_transformed.py
+++ b/test/test_extreme_estimator/test_extreme_models/test_margin_temporal_transformed.py
@@ -41,8 +41,8 @@ class TestMarginTemporalTransformed(unittest.TestCase):
                'shape': 0.7289248773961512}
         for year in range(1, 3):
             coordinate = np.array([0.0, 0.0, year])
-            mle_params_estimated = estimator.margin_function_fitted.get_gev_params(coordinate,
-                                                                                   is_transformed=False).to_dict()
+            mle_params_estimated = estimator.margin_function_from_fit.get_gev_params(coordinate,
+                                                                                     is_transformed=False).to_dict()
             self.assertEqual(mle_params_estimated, ref)
 
     def test_margin_fit_nonstationary(self):
@@ -50,32 +50,32 @@ class TestMarginTemporalTransformed(unittest.TestCase):
         margin_model = LinearNonStationaryLocationMarginModel(self.coordinates)
         estimator = LinearMarginEstimator(self.dataset, margin_model)
         estimator.fit()
-        self.assertNotEqual(estimator.margin_function_fitted.mu1_temporal_trend, 0.0)
+        self.assertNotEqual(estimator.margin_function_from_fit.mu1_temporal_trend, 0.0)
         # Checks that parameters returned are indeed different
         coordinate1 = np.array([0.0, 0.0, 1])
-        mle_params_estimated_year1 = estimator.margin_function_fitted.get_gev_params(coordinate1,
-                                                                                     is_transformed=False).to_dict()
+        mle_params_estimated_year1 = estimator.margin_function_from_fit.get_gev_params(coordinate1,
+                                                                                       is_transformed=False).to_dict()
         coordinate3 = np.array([0.0, 0.0, 3])
-        mle_params_estimated_year3 = estimator.margin_function_fitted.get_gev_params(coordinate3,
-                                                                                     is_transformed=False).to_dict()
+        mle_params_estimated_year3 = estimator.margin_function_from_fit.get_gev_params(coordinate3,
+                                                                                       is_transformed=False).to_dict()
         self.assertNotEqual(mle_params_estimated_year1, mle_params_estimated_year3)
 
     def test_margin_fit_nonstationary_with_start_point(self):
         # Create estimator
         estimator = self.fit_non_stationary_estimator(starting_point=2)
         # By default, estimator find the good margin
-        self.assertNotEqual(estimator.margin_function_fitted.mu1_temporal_trend, 0.0)
+        self.assertNotEqual(estimator.margin_function_from_fit.mu1_temporal_trend, 0.0)
         # Checks that parameters returned are indeed different
         coordinate1 = np.array([0.0, 0.0, 1])
-        mle_params_estimated_year1 = estimator.margin_function_fitted.get_gev_params(coordinate1,
-                                                                                     is_transformed=False).to_dict()
+        mle_params_estimated_year1 = estimator.margin_function_from_fit.get_gev_params(coordinate1,
+                                                                                       is_transformed=False).to_dict()
         coordinate2 = np.array([0.0, 0.0, 2])
-        mle_params_estimated_year2 = estimator.margin_function_fitted.get_gev_params(coordinate2,
-                                                                                     is_transformed=False).to_dict()
+        mle_params_estimated_year2 = estimator.margin_function_from_fit.get_gev_params(coordinate2,
+                                                                                       is_transformed=False).to_dict()
         self.assertEqual(mle_params_estimated_year1, mle_params_estimated_year2)
         coordinate5 = np.array([0.0, 0.0, 5])
-        mle_params_estimated_year5 = estimator.margin_function_fitted.get_gev_params(coordinate5,
-                                                                                     is_transformed=False).to_dict()
+        mle_params_estimated_year5 = estimator.margin_function_from_fit.get_gev_params(coordinate5,
+                                                                                       is_transformed=False).to_dict()
         self.assertNotEqual(mle_params_estimated_year5, mle_params_estimated_year2)
 
     def fit_non_stationary_estimator(self, starting_point):
@@ -88,8 +88,8 @@ class TestMarginTemporalTransformed(unittest.TestCase):
         # Create two different estimators
         estimator1 = self.fit_non_stationary_estimator(starting_point=3)
         estimator2 = self.fit_non_stationary_estimator(starting_point=20)
-        mu1_estimator1 = estimator1.margin_function_fitted.mu1_temporal_trend
-        mu1_estimator2 = estimator2.margin_function_fitted.mu1_temporal_trend
+        mu1_estimator1 = estimator1.margin_function_from_fit.mu1_temporal_trend
+        mu1_estimator2 = estimator2.margin_function_from_fit.mu1_temporal_trend
         self.assertNotEqual(mu1_estimator1, mu1_estimator2)
 
 
diff --git a/test/test_extreme_estimator/test_extreme_models/test_max_stable_temporal.py b/test/test_extreme_estimator/test_extreme_models/test_max_stable_temporal.py
index c48667bd491c7d3ac62843ccac3a91696671072b..07570582b6c40c751639eb3b64ddb106b8a07ce2 100644
--- a/test/test_extreme_estimator/test_extreme_models/test_max_stable_temporal.py
+++ b/test/test_extreme_estimator/test_extreme_models/test_max_stable_temporal.py
@@ -37,7 +37,7 @@ class TestMaxStableTemporal(unittest.TestCase):
         ref = {'loc': 1.2091156634312243, 'scale': 1.1210085591373455, 'shape': 0.9831957705294134}
         for year in range(1, 3):
             coordinate = np.array([0.0, 0.0, year])
-            mle_params_estimated = estimator.margin_function_fitted.get_gev_params(coordinate).to_dict()
+            mle_params_estimated = estimator.margin_function_from_fit.get_gev_params(coordinate).to_dict()
             for key in ref.keys():
                 self.assertAlmostEqual(ref[key], mle_params_estimated[key], places=3)
 
@@ -47,32 +47,32 @@ class TestMaxStableTemporal(unittest.TestCase):
         estimator = FullEstimatorInASingleStepWithSmoothMargin(self.dataset, margin_model,
                                                                self.max_stable_model)
         estimator.fit()
-        self.assertNotEqual(estimator.margin_function_fitted.mu1_temporal_trend, 0.0)
+        self.assertNotEqual(estimator.margin_function_from_fit.mu1_temporal_trend, 0.0)
         # Checks that parameters returned are indeed different
         coordinate1 = np.array([0.0, 0.0, 1])
-        mle_params_estimated_year1 = estimator.margin_function_fitted.get_gev_params(coordinate1).to_dict()
+        mle_params_estimated_year1 = estimator.margin_function_from_fit.get_gev_params(coordinate1).to_dict()
         coordinate3 = np.array([0.0, 0.0, 3])
-        mle_params_estimated_year3 = estimator.margin_function_fitted.get_gev_params(coordinate3).to_dict()
+        mle_params_estimated_year3 = estimator.margin_function_from_fit.get_gev_params(coordinate3).to_dict()
         self.assertNotEqual(mle_params_estimated_year1, mle_params_estimated_year3)
 
     def test_margin_fit_nonstationary_with_start_point(self):
         # Create estimator
         estimator = self.fit_non_stationary_estimator(starting_point=2)
         # By default, estimator find the good margin
-        self.assertNotEqual(estimator.margin_function_fitted.mu1_temporal_trend, 0.0)
-        self.assertAlmostEqual(estimator.margin_function_fitted.mu1_temporal_trend,
+        self.assertNotEqual(estimator.margin_function_from_fit.mu1_temporal_trend, 0.0)
+        self.assertAlmostEqual(estimator.margin_function_from_fit.mu1_temporal_trend,
                                self.smooth_margin_model.margin_function_sample.mu1_temporal_trend,
                                places=2)
         # Checks starting point parameter are well passed
-        self.assertEqual(2, estimator.margin_function_fitted.starting_point)
+        self.assertEqual(2, estimator.margin_function_from_fit.starting_point)
         # Checks that parameters returned are indeed different
         coordinate1 = np.array([0.0, 0.0, 1])
-        mle_params_estimated_year1 = estimator.margin_function_fitted.get_gev_params(coordinate1).to_dict()
+        mle_params_estimated_year1 = estimator.margin_function_from_fit.get_gev_params(coordinate1).to_dict()
         coordinate2 = np.array([0.0, 0.0, 2])
-        mle_params_estimated_year2 = estimator.margin_function_fitted.get_gev_params(coordinate2).to_dict()
+        mle_params_estimated_year2 = estimator.margin_function_from_fit.get_gev_params(coordinate2).to_dict()
         self.assertEqual(mle_params_estimated_year1, mle_params_estimated_year2)
         coordinate5 = np.array([0.0, 0.0, 5])
-        mle_params_estimated_year5 = estimator.margin_function_fitted.get_gev_params(coordinate5).to_dict()
+        mle_params_estimated_year5 = estimator.margin_function_from_fit.get_gev_params(coordinate5).to_dict()
         self.assertNotEqual(mle_params_estimated_year5, mle_params_estimated_year2)
 
     def fit_non_stationary_estimator(self, starting_point):
@@ -86,8 +86,8 @@ class TestMaxStableTemporal(unittest.TestCase):
         # Create two different estimators
         estimator1 = self.fit_non_stationary_estimator(starting_point=3)
         estimator2 = self.fit_non_stationary_estimator(starting_point=20)
-        mu1_estimator1 = estimator1.margin_function_fitted.mu1_temporal_trend
-        mu1_estimator2 = estimator2.margin_function_fitted.mu1_temporal_trend
+        mu1_estimator1 = estimator1.margin_function_from_fit.mu1_temporal_trend
+        mu1_estimator2 = estimator2.margin_function_from_fit.mu1_temporal_trend
         self.assertNotEqual(mu1_estimator1, mu1_estimator2)
 
 
diff --git a/test/test_extreme_estimator/test_margin_fits/test_gev/test_gev_temporal.py b/test/test_extreme_estimator/test_margin_fits/test_gev/test_gev_temporal.py
index 64403ba019972cf7dc9d8b902235cde74acd6c50..3b7113fab76777e75e36c1c26a3db9d3837930d7 100644
--- a/test/test_extreme_estimator/test_margin_fits/test_gev/test_gev_temporal.py
+++ b/test/test_extreme_estimator/test_margin_fits/test_gev/test_gev_temporal.py
@@ -41,7 +41,7 @@ class TestGevTemporal(unittest.TestCase):
         estimator.fit()
         ref = {'loc': 0.0219, 'scale': 1.0347, 'shape': 0.8295}
         for year in range(1, 3):
-            mle_params_estimated = estimator.margin_function_fitted.get_gev_params(np.array([year])).to_dict()
+            mle_params_estimated = estimator.margin_function_from_fit.get_gev_params(np.array([year])).to_dict()
             for key in ref.keys():
                 self.assertAlmostEqual(ref[key], mle_params_estimated[key], places=3)
 
@@ -53,21 +53,21 @@ class TestGevTemporal(unittest.TestCase):
             estimator = LinearMarginEstimator(self.dataset, margin_model)
             estimator.fit()
             # Checks that parameters returned are indeed different
-            mle_params_estimated_year1 = estimator.margin_function_fitted.get_gev_params(np.array([1])).to_dict()
-            mle_params_estimated_year3 = estimator.margin_function_fitted.get_gev_params(np.array([3])).to_dict()
+            mle_params_estimated_year1 = estimator.margin_function_from_fit.get_gev_params(np.array([1])).to_dict()
+            mle_params_estimated_year3 = estimator.margin_function_from_fit.get_gev_params(np.array([3])).to_dict()
             self.assertNotEqual(mle_params_estimated_year1, mle_params_estimated_year3)
 
     def test_gev_temporal_margin_fit_nonstationary_with_start_point(self):
         # Create estimator
         estimator = self.fit_non_stationary_estimator(starting_point=3)
-        self.assertNotEqual(estimator.margin_function_fitted.mu1_temporal_trend, 0.0)
+        self.assertNotEqual(estimator.margin_function_from_fit.mu1_temporal_trend, 0.0)
         # Checks starting point parameter are well passed
-        self.assertEqual(3, estimator.margin_function_fitted.starting_point)
+        self.assertEqual(3, estimator.margin_function_from_fit.starting_point)
         # Checks that parameters returned are indeed different
-        mle_params_estimated_year1 = estimator.margin_function_fitted.get_gev_params(np.array([1])).to_dict()
-        mle_params_estimated_year3 = estimator.margin_function_fitted.get_gev_params(np.array([3])).to_dict()
+        mle_params_estimated_year1 = estimator.margin_function_from_fit.get_gev_params(np.array([1])).to_dict()
+        mle_params_estimated_year3 = estimator.margin_function_from_fit.get_gev_params(np.array([3])).to_dict()
         self.assertEqual(mle_params_estimated_year1, mle_params_estimated_year3)
-        mle_params_estimated_year5 = estimator.margin_function_fitted.get_gev_params(np.array([5])).to_dict()
+        mle_params_estimated_year5 = estimator.margin_function_from_fit.get_gev_params(np.array([5])).to_dict()
         self.assertNotEqual(mle_params_estimated_year5, mle_params_estimated_year3)
 
     def fit_non_stationary_estimator(self, starting_point):
@@ -80,8 +80,8 @@ class TestGevTemporal(unittest.TestCase):
         # Create two different estimators
         estimator1 = self.fit_non_stationary_estimator(starting_point=3)
         estimator2 = self.fit_non_stationary_estimator(starting_point=28)
-        mu1_estimator1 = estimator1.margin_function_fitted.mu1_temporal_trend
-        mu1_estimator2 = estimator2.margin_function_fitted.mu1_temporal_trend
+        mu1_estimator1 = estimator1.margin_function_from_fit.mu1_temporal_trend
+        mu1_estimator2 = estimator2.margin_function_from_fit.mu1_temporal_trend
         self.assertNotEqual(mu1_estimator1, mu1_estimator2)
 
 
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 86d74df487f43e3bdba8506d44e2d29576f692e2..c4c160527a2dc6bcfc194b8f9a314921f6dc7c00 100644
--- a/test/test_unitary/test_fitmaxstab/test_fitmaxstab_with_margin.py
+++ b/test/test_unitary/test_fitmaxstab/test_fitmaxstab_with_margin.py
@@ -32,7 +32,7 @@ class TestMaxStableFitWithConstantMargin(TestUnitaryAbstract):
         full_estimator = FullEstimatorInASingleStepWithSmoothMargin(dataset, margin_model,
                                                                     max_stable_model)
         full_estimator.fit()
-        return full_estimator.result_from_fit.all_parameters
+        return full_estimator.result_from_model_fit.all_parameters
 
     def test_max_stable_fit_with_constant_margin(self):
         self.compare()
@@ -59,7 +59,7 @@ class TestMaxStableFitWithLinearMargin(TestUnitaryAbstract):
         full_estimator = FullEstimatorInASingleStepWithSmoothMargin(dataset, margin_model,
                                                                     max_stable_model)
         full_estimator.fit()
-        return full_estimator.result_from_fit.all_parameters
+        return full_estimator.result_from_model_fit.all_parameters
 
     def test_max_stable_fit_with_linear_margin(self):
         self.compare()