diff --git a/experiment/regression_margin/regression_margin.py b/experiment/regression_margin/regression_margin.py
index e7089feeadcf3771b4243ce6522c8e0f4369a46e..bb96a5dfbae56a181a5409a10092ef6bedd1233e 100644
--- a/experiment/regression_margin/regression_margin.py
+++ b/experiment/regression_margin/regression_margin.py
@@ -10,12 +10,12 @@ import matplotlib.pyplot as plt
 
 from spatio_temporal_dataset.dataset.simulation_dataset import FullSimulatedDataset
 
-nb_points = 50
-nb_obs = 100
+nb_points = 5
+nb_obs = 10
+show = False
 
 coordinates = LinSpaceCoordinates.from_nb_points(nb_points=nb_points)
 
-
 ########## GENERATING THE DATA #####################
 
 # MarginModel Linear with respect to the shape (from 0.01 to 0.02)
@@ -24,19 +24,20 @@ max_stable_model = Smith()
 dataset = FullSimulatedDataset.from_double_sampling(nb_obs=nb_obs, margin_model=margin_model,
                                                     coordinates=coordinates,
                                                     max_stable_model=max_stable_model)
-# Visualize the sampling margin
-dataset.margin_model.margin_function_sample.visualize_all()
-# Plot a realization from the maxima gev (i.e the maxima obtained just by simulating the marginal law)
-for maxima_gev in np.transpose(dataset.maxima_gev):
-    plt.plot(coordinates.coordinates_values, maxima_gev)
-plt.show()
+if show:
+    # Visualize the sampling margin
+    dataset.margin_model.margin_function_sample.visualize_all()
+    # Plot a realization from the maxima gev (i.e the maxima obtained just by simulating the marginal law)
+    for maxima_gev in np.transpose(dataset.maxima_gev):
+        plt.plot(coordinates.coordinates_values, maxima_gev)
+    plt.show()
 
 ######### FITTING A MODEL #################
 
 margin_model = LinearAllParametersAllAxisMarginModel(coordinates)
 full_estimator = FullEstimatorInASingleStepWithSmoothMargin(dataset, margin_model, max_stable_model)
 full_estimator.fit()
-print(full_estimator.full_params_fitted)
+full_estimator.margin_function_fitted.visualize_all()
 
 # Plot the margin functions
 # margin_model.margin_function_sample.visualize_2D()
diff --git a/extreme_estimator/estimator/full_estimator.py b/extreme_estimator/estimator/full_estimator.py
index 1fbad39888e2ce8cb6857e8943046d3bb7274b8e..35d9ca96edadbe49fcc6bf3123409e6050af714d 100644
--- a/extreme_estimator/estimator/full_estimator.py
+++ b/extreme_estimator/estimator/full_estimator.py
@@ -1,6 +1,7 @@
 from extreme_estimator.extreme_models.margin_model.abstract_margin_model import AbstractMarginModel
 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 extreme_estimator.extreme_models.margin_model.smooth_margin_model import LinearMarginModel
 from extreme_estimator.extreme_models.max_stable_model.abstract_max_stable_model import AbstractMaxStableModel
 from extreme_estimator.estimator.abstract_estimator import AbstractEstimator
@@ -60,7 +61,9 @@ class FullEstimatorInASingleStepWithSmoothMargin(AbstractFullEstimator):
                  max_stable_model: AbstractMaxStableModel):
         super().__init__(dataset)
         self.max_stable_model = max_stable_model
-        self.smooth_margin_function_to_fit = margin_model.margin_function_start_fit
+        self.linear_margin_model = margin_model
+        self.linear_margin_function_to_fit = self.linear_margin_model.margin_function_start_fit
+        assert isinstance(self.linear_margin_function_to_fit, LinearMarginFunction)
 
     def _fit(self):
         # Estimate both the margin and the max-stable structure
@@ -68,11 +71,14 @@ class FullEstimatorInASingleStepWithSmoothMargin(AbstractFullEstimator):
             maxima_frech=self.dataset.maxima_frech,
             df_coordinates=self.dataset.df_coordinates,
             fit_marge=True,
-            fit_marge_form_dict=self.smooth_margin_function_to_fit.form_dict,
-            margin_start_dict=self.smooth_margin_function_to_fit.coef_dict
+            fit_marge_form_dict=self.linear_margin_function_to_fit.form_dict,
+            margin_start_dict=self.linear_margin_function_to_fit.coef_dict
         )
-        # Initialize
-        # self._margin_function_fitted =
+        # 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)
 
 
 class PointwiseAndThenUnitaryMsp(AbstractFullEstimator):
diff --git a/extreme_estimator/extreme_models/margin_model/margin_function/linear_margin_function.py b/extreme_estimator/extreme_models/margin_model/margin_function/linear_margin_function.py
index 62377cea4185370ccbd4ea95277d969e4c8475bc..2d50202ea120a1e0c9aa7157a0b5af05dfe10a33 100644
--- a/extreme_estimator/extreme_models/margin_model/margin_function/linear_margin_function.py
+++ b/extreme_estimator/extreme_models/margin_model/margin_function/linear_margin_function.py
@@ -1,4 +1,4 @@
-from typing import Dict, List, Tuple
+from typing import Dict, List
 
 from extreme_estimator.extreme_models.margin_model.margin_function.independent_margin_function import \
     IndependentMarginFunction
@@ -29,8 +29,8 @@ class LinearMarginFunction(IndependentMarginFunction):
                  gev_param_name_to_linear_dims: Dict[str, List[int]],
                  gev_param_name_to_linear_coef: Dict[str, LinearCoef]):
         super().__init__(coordinates)
-        self.gev_param_name_to_linear_coef = gev_param_name_to_linear_coef
-        self.gev_param_name_to_linear_dims = gev_param_name_to_linear_dims
+        self.gev_param_name_to_linear_coef = gev_param_name_to_linear_coef  # type: Dict[str, LinearCoef]
+        self.gev_param_name_to_linear_dims = gev_param_name_to_linear_dims  # type: Dict[str, List[int]]
         # Build gev_parameter_to_param_function dictionary
         self.gev_param_name_to_param_function = {}  # type: Dict[str, ParamFunction]
 
@@ -54,11 +54,18 @@ class LinearMarginFunction(IndependentMarginFunction):
             self.gev_param_name_to_param_function[gev_param_name] = param_function
 
     @classmethod
-    def from_margin_coeff_dict(cls, coordinates, margin_coeff_dict):
-        pass
+    def from_coef_dict(cls, coordinates: AbstractCoordinates, gev_param_name_to_linear_dims: Dict[str, List[int]],
+                       coef_dict: Dict[str, float]):
+        gev_param_name_to_linear_coef = {}
+        for gev_param_name in GevParams.GEV_PARAM_NAMES:
+            linear_dims = gev_param_name_to_linear_dims.get(gev_param_name, [])
+            linear_coef = LinearCoef.from_coef_dict(coef_dict=coef_dict, gev_param_name=gev_param_name,
+                                                    linear_dims=linear_dims)
+            gev_param_name_to_linear_coef[gev_param_name] = linear_coef
+        return cls(coordinates, gev_param_name_to_linear_dims, gev_param_name_to_linear_coef)
 
     @property
-    def form_dict(self) -> dict:
+    def form_dict(self) -> Dict[str, str]:
         form_dict = {}
         for gev_param_name in GevParams.GEV_PARAM_NAMES:
             linear_dims = self.gev_param_name_to_linear_dims.get(gev_param_name, [])
@@ -66,7 +73,7 @@ class LinearMarginFunction(IndependentMarginFunction):
         return form_dict
 
     @property
-    def coef_dict(self) -> dict:
+    def coef_dict(self) -> Dict[str, float]:
         coef_dict = {}
         for gev_param_name in GevParams.GEV_PARAM_NAMES:
             linear_dims = self.gev_param_name_to_linear_dims.get(gev_param_name, [])
diff --git a/extreme_estimator/extreme_models/margin_model/param_function/linear_coef.py b/extreme_estimator/extreme_models/margin_model/param_function/linear_coef.py
index 8374737ddd8ddc99391c1bcdad04e45dc4cc2551..e6c825705d45ef9880420b659a6b398934dea6c4 100644
--- a/extreme_estimator/extreme_models/margin_model/param_function/linear_coef.py
+++ b/extreme_estimator/extreme_models/margin_model/param_function/linear_coef.py
@@ -12,7 +12,7 @@ class LinearCoef(object):
         dim = 3 correspond to the coordinate Z
     """
 
-    def __init__(self, gev_param_name: str, default_value: float = 0.0, dim_to_coef: Dict[int, float] = None):
+    def __init__(self, gev_param_name: str, dim_to_coef: Dict[int, float] = None, default_value: float = 0.0):
         self.gev_param_name = gev_param_name
         self.dim_to_coef = dim_to_coef
         self.default_value = default_value
@@ -27,21 +27,28 @@ class LinearCoef(object):
     def intercept(self):
         return self.get_coef(dim=0)
 
-    @classmethod
-    def from_dict(cls, coef_dict: Dict[int, float], gev_param_name: str, default_value: float = 0.0):
-        pass
+    @staticmethod
+    def coef_template_str(gev_param_name):
+        return gev_param_name + 'Coeff{}'
 
-    def coef_dict(self, linear_dims):
-        coef_dict = {}
-        coef_template_str = self.gev_param_name + 'Coeff{}'
+    @classmethod
+    def from_coef_dict(cls, coef_dict: Dict[str, float], gev_param_name: str, linear_dims):
+        dims = [0] + linear_dims
+        dim_to_coef = {}
+        for j, dim in enumerate(dims, 1):
+            coef = coef_dict[cls.coef_template_str(gev_param_name).format(j)]
+            dim_to_coef[dim] = coef
+        return cls(gev_param_name, dim_to_coef)
+
+    def coef_dict(self, linear_dims) -> Dict[str, float]:
         # Constant param must be specified for all the parameters
-        coef_dict[coef_template_str.format(1)] = self.intercept
+        coef_dict = {self.coef_template_str(self.gev_param_name).format(1): self.intercept}
         # Specify only the param that belongs to dim_to_coef
         for j, dim in enumerate(linear_dims, 2):
-            coef_dict[coef_template_str.format(j)] = self.dim_to_coef[dim]
+            coef_dict[self.coef_template_str(self.gev_param_name).format(j)] = self.dim_to_coef[dim]
         return coef_dict
 
-    def form_dict(self, linear_dims):
+    def form_dict(self, linear_dims) -> Dict[str, str]:
         """
         Example of formula that could be specified:
         loc.form = loc ~ coord_x
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 d8051e8676d3eca14276397ff2010d97996eb615..a166af4daf5818cbc782c27c287d688682026539 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
@@ -32,6 +32,7 @@ class AbstractMaxStableModel(AbstractModel):
         data = np.transpose(maxima_frech)
 
         # Prepare the coord
+        df_coordinates = df_coordinates.copy()
         # In the one dimensional case, fitmaxstab isn't working
         # therefore, we treat our 1D coordinate as 2D coordinate on the line y=x, and enforce iso=TRUE
         fitmaxstab_with_one_dimensional_data = len(df_coordinates.columns) == 1