diff --git a/extreme_estimator/R_fit/abstract_model.py b/extreme_estimator/R_fit/abstract_model.py
new file mode 100644
index 0000000000000000000000000000000000000000..8c6c4d86733088a70ce2687826be66c53fb714a3
--- /dev/null
+++ b/extreme_estimator/R_fit/abstract_model.py
@@ -0,0 +1,11 @@
+from extreme_estimator.R_fit.utils import get_loaded_r
+
+
+class AbstractModel(object):
+
+    def __init__(self, params_start_fit=None, params_sample=None):
+        self.default_params_start_fit = None
+        self.default_params_sample = None
+        self.user_params_start_fit = params_start_fit
+        self.user_params_sample = params_sample
+        self.r = get_loaded_r()
\ No newline at end of file
diff --git a/extreme_estimator/R_fit/gev_fit/abstract_margin_model.py b/extreme_estimator/R_fit/gev_fit/abstract_margin_model.py
index dc3e41a7ef16a8d90576d8cd54671a41b6b87765..e2c6aa8a59b7696d7e7ec1c185572ee0954fcfcc 100644
--- a/extreme_estimator/R_fit/gev_fit/abstract_margin_model.py
+++ b/extreme_estimator/R_fit/gev_fit/abstract_margin_model.py
@@ -1,44 +1,18 @@
-from extreme_estimator.R_fit.gev_fit.gev_mle_fit import GevMleFit, mle_gev
-from extreme_estimator.R_fit.utils import get_loaded_r
-
-
-def frechet_unitary_transformation(data, location, scale, shape):
-    """
-    Compute the unitary Frechet transformed data
-    (1 + \zeta \frac{z - \mu}{\sigma}) ^ {\frac{1}{\zeta}}
-    """
-    assert False
-    # todo: there is already a function doing that in R
-    return (1 + shape * (data - location) / scale) ** (1 / shape)
-
-
-class GevParameters(object):
-
-    def __init__(self, location, scale, shape):
-        self.location = location
-        self.scale = scale
-        self.shape = shape
+import numpy as np
+import pandas as pd
 
-
-def frechet_unitary_transformation_from_gev_parameters(data, gev_parameters: GevParameters):
-    return frechet_unitary_transformation(data, gev_parameters.location)
+from extreme_estimator.R_fit.abstract_model import AbstractModel
+from extreme_estimator.R_fit.gev_fit.gev_mle_fit import GevMleFit, mle_gev
 
 
-class AbstractMarginModel(object):
+class AbstractMarginModel(AbstractModel):
     GEV_SCALE = GevMleFit.GEV_SCALE
     GEV_LOCATION = GevMleFit.GEV_LOCATION
     GEV_SHAPE = GevMleFit.GEV_SHAPE
     GEV_PARAMETERS = [GEV_LOCATION, GEV_SCALE, GEV_SHAPE]
 
-    def __init__(self):
-        """
-        Class to fit a GEV a list of data. Compute also the corresponding unitary data
-
-        :param coordinate: Represents the spatio-temporal spatial_coordinates of the marginals
-        :param data:  array of data corresponding to this position (and potentially its neighborhood)
-        """
-        self.default_params = {gev_param: 1.0 for gev_param in self.GEV_PARAMETERS}
-        self.r = get_loaded_r()
+    def __init__(self, params_start_fit=None, params_sample=None):
+        super().__init__(params_start_fit, params_sample)
 
     # Define the method to sample/fit a single gev
 
@@ -53,20 +27,37 @@ class AbstractMarginModel(object):
     def fitgev(self, x_gev, estimator=GevMleFit):
         mle_params = mle_gev(x_gev=x_gev)
 
+    def gev_params_sample(self, coordinate) -> dict:
+        pass
+
     # Define the method to sample/fit all marginals globally in the child classes
 
     def fitmargin(self, maxima, coord):
-        df_gev_params = None
-        return df_gev_params
+        df_fit_gev_params = None
+        return df_fit_gev_params
 
     def rmargin(self, nb_obs, coord):
-        pass
-
-    def get_maxima(self, maxima_normalized, coord):
-        pass
-
-    def get_maxima_normalized(self, maxima, df_gev_params):
-        pass
+        maxima_gev = None
+        return maxima_gev
+
+    def frech2gev(self, maxima_frech: np.ndarray, coordinates: np.ndarray):
+        assert len(maxima_frech) == len(coordinates)
+        maxima_gev = []
+        for x_frech, coordinate in zip(maxima_frech, coordinates):
+            gev_params = self.gev_params_sample(coordinate)
+            x_gev = self.r.frech2gev(x_frech, **gev_params)
+            maxima_gev.append(x_gev)
+        return np.array(maxima_gev)
+
+    def gev2frech(self, maxima_gev: np.ndarray, df_gev_params: pd.DataFrame):
+        assert len(maxima_gev) == len(df_gev_params)
+        maxima_frech = []
+        for x_gev, (_, s_gev_params) in zip(maxima_gev, df_gev_params.iterrows()):
+            gev_params = dict(s_gev_params)
+            gev2frech_param = {'emp': False}
+            x_frech = self.r.gev2frech(x_gev, **gev_params, **gev2frech_param)
+            maxima_frech.append(x_frech)
+        return np.array(maxima_frech)
 
 
 class SmoothMarginModel(AbstractMarginModel):
@@ -74,4 +65,17 @@ class SmoothMarginModel(AbstractMarginModel):
 
 
 class ConstantMarginModel(SmoothMarginModel):
-    pass
+    def __init__(self, params_start_fit=None, params_sample=None):
+        super().__init__(params_start_fit, params_sample)
+        self.default_params_sample = {gev_param: 1.0 for gev_param in self.GEV_PARAMETERS}
+        self.default_params_start_fit = {gev_param: 1.0 for gev_param in self.GEV_PARAMETERS}
+
+    def gev_params_sample(self, coordinate):
+        return self.default_params_sample
+
+    def fitmargin(self, maxima, coord):
+        return pd.DataFrame([pd.Series(self.default_params_start_fit) for _ in maxima])
+
+
+
+
diff --git a/extreme_estimator/R_fit/max_stable_fit/abstract_max_stable_model.py b/extreme_estimator/R_fit/max_stable_fit/abstract_max_stable_model.py
index 8009570bb0a018a5ea11f11a2c5c56d397f891bc..67e4ab5b21ca3a5c93035abe0fbeb3976996207f 100644
--- a/extreme_estimator/R_fit/max_stable_fit/abstract_max_stable_model.py
+++ b/extreme_estimator/R_fit/max_stable_fit/abstract_max_stable_model.py
@@ -1,23 +1,20 @@
-import rpy2
 from enum import Enum
 
-from rpy2.robjects import ListVector
-from extreme_estimator.R_fit.utils import get_loaded_r
 import numpy as np
+import rpy2
+from rpy2.robjects import ListVector
 
+from extreme_estimator.R_fit.abstract_model import AbstractModel
 
-class AbstractMaxStableModel(object):
+
+class AbstractMaxStableModel(AbstractModel):
 
     def __init__(self, params_start_fit=None, params_sample=None):
+        super().__init__(params_start_fit, params_sample)
         self.cov_mod = None
-        self.default_params_start_fit = None
-        self.default_params_sample = None
-        self.user_params_start_fit = params_start_fit
-        self.user_params_sample = params_sample
-        self.r = get_loaded_r()
-
-    def fitmaxstab(self, maxima_normalized: np.ndarray, coord: np.ndarray, fit_marge=False):
-        assert all([isinstance(arr, np.ndarray) for arr in [maxima_normalized, coord]])
+
+    def fitmaxstab(self, maxima_frech: np.ndarray, coord: np.ndarray, fit_marge=False):
+        assert all([isinstance(arr, np.ndarray) for arr in [maxima_frech, coord]])
         #  Specify the fit params
         fit_params = {
             'fit.marge': fit_marge,
@@ -26,7 +23,8 @@ class AbstractMaxStableModel(object):
         # Run the fitmaxstab in R
         # todo: find how to specify the optim function to use
         try:
-            res = self.r.fitmaxstab(np.transpose(maxima_normalized), coord, **self.cov_mod_param, **fit_params)  # type: ListVector
+            res = self.r.fitmaxstab(np.transpose(maxima_frech), coord, **self.cov_mod_param,
+                                    **fit_params)  # type: ListVector
         except rpy2.rinterface.RRuntimeError as error:
             raise Exception('Some R exception have been launched at RunTime: {}'.format(error.__repr__()))
         # todo: maybe if the convergence was not successful I could try other starting point several times
@@ -35,12 +33,13 @@ class AbstractMaxStableModel(object):
         fitted_values = {key: fitted_values.rx2(key)[0] for key in fitted_values.names}
         return fitted_values
 
-    def rmaxstab(self, nb_obs: int, coord: np.ndarray, ) -> np.ndarray:
+    def rmaxstab(self, nb_obs: int, coord: np.ndarray) -> np.ndarray:
         """
         Return an numpy of maxima. With rows being the stations and columns being the years of maxima
         """
-        maxima = np.array(self.r.rmaxstab(nb_obs, coord, *list(self.cov_mod_param.values()), **self.params_sample))
-        return np.transpose(maxima)
+        maxima_frech = np.array(
+            self.r.rmaxstab(nb_obs, coord, *list(self.cov_mod_param.values()), **self.params_sample))
+        return np.transpose(maxima_frech)
 
     @property
     def cov_mod_param(self):
@@ -82,4 +81,4 @@ class AbstractMaxStableModelWithCovarianceFunction(AbstractMaxStableModel):
             'range': 3,
             'smooth': 0.5,
             'nugget': 0.5
-        }
\ No newline at end of file
+        }
diff --git a/extreme_estimator/R_fit/max_stable_fit/max_stable_fit.R b/extreme_estimator/R_fit/max_stable_fit/max_stable_fit.R
index 5a3754e9a9cad8b83c1ba62256d693f83e04a789..5c0491e8b4b948f6df4f378097f7c7ab4fc91c2f 100644
--- a/extreme_estimator/R_fit/max_stable_fit/max_stable_fit.R
+++ b/extreme_estimator/R_fit/max_stable_fit/max_stable_fit.R
@@ -24,6 +24,8 @@ if (call_main) {
     # res = fitmaxstab(data, coord, "whitmat", par())
     # print(res)
     # print(class(res))
+    gev2frech(x, loc, scale, shape, emp=FALSE)
+    frech2gev(x, loc, scale, shape)
     # print(names(res))
     for (name in names(res)){
         print(name)
diff --git a/extreme_estimator/estimator/full_estimator.py b/extreme_estimator/estimator/full_estimator.py
index 33fe652827cef315620179248d0faaf71e62163b..907ba60f6b726ef73579a26e5f40ef8a74884829 100644
--- a/extreme_estimator/estimator/full_estimator.py
+++ b/extreme_estimator/estimator/full_estimator.py
@@ -22,12 +22,13 @@ class SmoothMarginalsThenUnitaryMsp(AbstractFullEstimator):
     def _fit(self):
         # Estimate the margin parameters
         self.margin_estimator.fit()
-        # Compute the maxima_normalized
-        maxima_normalized = self.margin_estimator.margin_model.get_maxima_normalized(maxima=self.dataset.maxima,
-                                                                                     df_gev_params=self.margin_estimator.df_gev_params)
-        # Update maxima normalized field through the dataset object
-        print(maxima_normalized)
-        self.dataset.maxima_normalized = maxima_normalized
+        # Compute the maxima_frech
+        print(self.dataset.maxima_gev)
+        maxima_frech = self.margin_estimator.margin_model.gev2frech(maxima_gev=self.dataset.maxima_gev,
+                                                                    df_gev_params=self.margin_estimator.df_gev_params)
+        print(maxima_frech)
+        # Update maxima frech field through the dataset object
+        self.dataset.maxima_frech = maxima_frech
         # Estimate the max stable parameters
         self.max_stable_estimator.fit()
 
diff --git a/extreme_estimator/estimator/margin_estimator.py b/extreme_estimator/estimator/margin_estimator.py
index 2f78390baeeb8796b1ce697c145884a5e55419f4..0c997032c94b85e94e6da5856edbb0c0c74c18a2 100644
--- a/extreme_estimator/estimator/margin_estimator.py
+++ b/extreme_estimator/estimator/margin_estimator.py
@@ -7,7 +7,7 @@ class AbstractMarginEstimator(AbstractEstimator):
 
     def __init__(self, dataset: AbstractDataset):
         super().__init__(dataset)
-        assert dataset.temporal_observations.df_maxima is not None
+        assert self.dataset.maxima_gev is not None
 
 
 class PointWiseMarginEstimator(AbstractMarginEstimator):
@@ -23,5 +23,5 @@ class SmoothMarginEstimator(AbstractMarginEstimator):
         self.df_gev_params = None
 
     def _fit(self):
-        self.df_gev_params = self.margin_model.fitmargin(maxima=self.dataset.maxima,
-                                                         coord=self.dataset.coord)
+        self.df_gev_params = self.margin_model.fitmargin(maxima=self.dataset.maxima_gev,
+                                                         coord=self.dataset.coordinates)
diff --git a/extreme_estimator/estimator/max_stable_estimator.py b/extreme_estimator/estimator/max_stable_estimator.py
index 22d1fbaeca4a2b8916b9e804724af88af664f37f..96322572419ebcaf79ebee15c9f620df3b2e6a98 100644
--- a/extreme_estimator/estimator/max_stable_estimator.py
+++ b/extreme_estimator/estimator/max_stable_estimator.py
@@ -16,10 +16,10 @@ class AbstractMaxStableEstimator(AbstractEstimator):
 class MaxStableEstimator(AbstractMaxStableEstimator):
 
     def _fit(self):
-        assert self.dataset.maxima_normalized is not None
+        assert self.dataset.maxima_frech is not None
         self.max_stable_params_fitted = self.max_stable_model.fitmaxstab(
-            maxima_normalized=self.dataset.maxima_normalized,
-            coord=self.dataset.coord)
+            maxima_frech=self.dataset.maxima_frech,
+            coord=self.dataset.coordinates)
 
     def _error(self, true_max_stable_params: dict):
         absolute_errors = {param_name: np.abs(param_true_value - self.max_stable_params_fitted[param_name])
diff --git a/spatio_temporal_dataset/dataset/abstract_dataset.py b/spatio_temporal_dataset/dataset/abstract_dataset.py
index 855389c6a886dfe6fb341dbdef57681ea8d409cc..33faaa6d5771fe11fa0618bd25b4b1be1b580440 100644
--- a/spatio_temporal_dataset/dataset/abstract_dataset.py
+++ b/spatio_temporal_dataset/dataset/abstract_dataset.py
@@ -9,7 +9,6 @@ from spatio_temporal_dataset.spatial_coordinates.abstract_spatial_coordinates im
 class AbstractDataset(object):
 
     def __init__(self, temporal_observations: AbstractTemporalObservations, spatial_coordinates: AbstractSpatialCoordinates):
-        # assert
         # is_same_index = temporal_observations.index == spatial_coordinates.index  # type: pd.Series
         # assert is_same_index.all()
         self.temporal_observations = temporal_observations
@@ -32,21 +31,21 @@ class AbstractDataset(object):
     @property
     def df_dataset(self) -> pd.DataFrame:
         # Merge dataframes with the maxima and with the coordinates
-        return self.temporal_observations.df_maxima.join(self.spatial_coordinates.df_coord)
+        return self.temporal_observations.df_maxima_gev.join(self.spatial_coordinates.df_coordinates)
 
     @property
-    def coord(self):
-        return self.spatial_coordinates.coord
+    def coordinates(self):
+        return self.spatial_coordinates.coordinates
 
     @property
-    def maxima(self) -> np.ndarray:
-        return self.temporal_observations.maxima
+    def maxima_gev(self) -> np.ndarray:
+        return self.temporal_observations.maxima_gev
 
     @property
-    def maxima_normalized(self):
-        return self.temporal_observations.maxima_normalized
+    def maxima_frech(self):
+        return self.temporal_observations.maxima_frech
 
-    @maxima_normalized.setter
-    def maxima_normalized(self, maxima_normalized_to_set):
-        self.temporal_observations.maxima_normalized = maxima_normalized_to_set
+    @maxima_frech.setter
+    def maxima_frech(self, maxima_frech_to_set):
+        self.temporal_observations.maxima_frech = maxima_frech_to_set
 
diff --git a/spatio_temporal_dataset/dataset/simulation_dataset.py b/spatio_temporal_dataset/dataset/simulation_dataset.py
index e76668b282d94ad2a775a73d310ced359bde5f3f..6a571271da98c652878457fdc71ffd6c305cd99f 100644
--- a/spatio_temporal_dataset/dataset/simulation_dataset.py
+++ b/spatio_temporal_dataset/dataset/simulation_dataset.py
@@ -48,8 +48,9 @@ class MarginDataset(SimulatedDataset):
 class FullSimulatedDataset(SimulatedDataset):
 
     @classmethod
-    def from_sampling(cls, nb_obs: int, max_stable_model: AbstractMaxStableModel,
-                      spatial_coordinates: AbstractSpatialCoordinates):
+    def from_double_sampling(cls, nb_obs: int, max_stable_model: AbstractMaxStableModel,
+                             spatial_coordinates: AbstractSpatialCoordinates,
+                             margin_model: AbstractMarginModel):
         temporal_obs = FullAnnualMaxima.from_double_sampling(nb_obs, max_stable_model,
-                                                             spatial_coordinates)
+                                                             spatial_coordinates, margin_model)
         return cls(temporal_obs, spatial_coordinates, max_stable_model)
diff --git a/spatio_temporal_dataset/spatial_coordinates/abstract_spatial_coordinates.py b/spatio_temporal_dataset/spatial_coordinates/abstract_spatial_coordinates.py
index 187e0635b326d60afe999aa8f7cbbe70ef9012e5..2def2b96ece46ca9a6a071fa0882dd02f7023774 100644
--- a/spatio_temporal_dataset/spatial_coordinates/abstract_spatial_coordinates.py
+++ b/spatio_temporal_dataset/spatial_coordinates/abstract_spatial_coordinates.py
@@ -17,21 +17,21 @@ class AbstractSpatialCoordinates(object):
     TRAIN_SPLIT_STR = 'train_split'
     TEST_SPLIT_STR = 'test_split'
 
-    def __init__(self, df_coord: pd.DataFrame, s_split: pd.Series = None):
-        self.df_coord = df_coord
+    def __init__(self, df_coordinates: pd.DataFrame, s_split: pd.Series = None):
+        self.df_coordinates = df_coordinates
         self.s_split = s_split
 
     @classmethod
     def from_df(cls, df: pd.DataFrame):
         #  X and Y coordinates must be defined
         assert cls.COORD_X in df.columns and cls.COORD_Y in df.columns
-        df_coord = df.loc[:, cls.coord_columns(df)]
+        df_coordinates = df.loc[:, cls.coordinates_columns(df)]
         # Potentially, a split column can be specified
         s_split = df[cls.COORD_SPLIT] if cls.COORD_SPLIT in df.columns else None
-        return cls(df_coord=df_coord, s_split=s_split)
+        return cls(df_coordinates=df_coordinates, s_split=s_split)
 
     @classmethod
-    def coord_columns(cls, df_coord: pd.DataFrame) -> List[str]:
+    def coordinates_columns(cls, df_coord: pd.DataFrame) -> List[str]:
         # If a Z coordinate is in the DataFrame, then
         coord_columns = [cls.COORD_X, cls.COORD_Y]
         if cls.COORD_Z in df_coord.columns:
@@ -40,12 +40,12 @@ class AbstractSpatialCoordinates(object):
 
     @property
     def columns(self):
-        return self.coord_columns(df_coord=self.df_coord)
+        return self.coordinates_columns(df_coord=self.df_coordinates)
 
     @property
     def df(self) -> pd.DataFrame:
         # Merged DataFrame of df_coord and s_split
-        return self.df_coord if self.s_split is None else self.df_coord.join(self.s_split)
+        return self.df_coordinates if self.s_split is None else self.df_coordinates.join(self.s_split)
 
     @classmethod
     def from_csv(cls, csv_path: str = None):
@@ -66,52 +66,52 @@ class AbstractSpatialCoordinates(object):
             df_sample = pd.DataFrame.sample(coordinates.df, n=nb_points)
             return cls.from_df(df=df_sample)
 
-    def df_coord_split(self, split_str: str) -> pd.DataFrame:
+    def df_coordinates_split(self, split_str: str) -> pd.DataFrame:
         assert self.s_split is not None
         ind = self.s_split == split_str
-        return self.df_coord.loc[ind]
+        return self.df_coordinates.loc[ind]
 
-    def coord_values(self, df_coord: pd.DataFrame) -> np.ndarray:
-        return df_coord.loc[:, self.coord_columns(df_coord)].values
+    def coordinates_values(self, df_coordinates: pd.DataFrame) -> np.ndarray:
+        return df_coordinates.loc[:, self.coordinates_columns(df_coordinates)].values
 
     @property
-    def coord(self) -> np.ndarray:
-        return self.coord_values(df_coord=self.df_coord)
+    def coordinates(self) -> np.ndarray:
+        return self.coordinates_values(df_coordinates=self.df_coordinates)
 
     @property
-    def coord_train(self) -> np.ndarray:
-        return self.coord_values(df_coord=self.df_coord_split(self.TRAIN_SPLIT_STR))
+    def coordinates_train(self) -> np.ndarray:
+        return self.coordinates_values(df_coordinates=self.df_coordinates_split(self.TRAIN_SPLIT_STR))
 
     @property
-    def coord_test(self) -> np.ndarray:
-        return self.coord_values(df_coord=self.df_coord_split(self.TEST_SPLIT_STR))
+    def coordinates_test(self) -> np.ndarray:
+        return self.coordinates_values(df_coordinates=self.df_coordinates_split(self.TEST_SPLIT_STR))
 
     @property
     def index(self):
-        return self.df_coord.index
+        return self.df_coordinates.index
 
     #  Visualization
 
     def visualization_2D(self):
-        x, y = self.coord[:, 0], self.coord[:, 1]
+        x, y = self.coordinates[:, 0], self.coordinates[:, 1]
         plt.scatter(x, y)
         plt.show()
 
     def visualization_3D(self):
-        assert len(self.coord_columns(self.df_coord)) == 3
+        assert len(self.coordinates_columns(self.df_coordinates)) == 3
         fig = plt.figure()
         ax = fig.add_subplot(111, projection='3d')  # type: Axes3D
-        x, y, z = self.coord[:, 0], self.coord[:, 1], self.coord[:, 2]
+        x, y, z = self.coordinates[:, 0], self.coordinates[:, 1], self.coordinates[:, 2]
         ax.scatter(x, y, z, marker='^')
         plt.show()
 
     #  Magic Methods
 
     def __len__(self):
-        return len(self.df_coord)
+        return len(self.df_coordinates)
 
     def __mul__(self, other: float):
-        self.df_coord *= other
+        self.df_coordinates *= other
         return self
 
     def __rmul__(self, other):
diff --git a/spatio_temporal_dataset/spatial_coordinates/transformed_coordinates.py b/spatio_temporal_dataset/spatial_coordinates/transformed_coordinates.py
index 6b19e87a4cc344cc5d5f4182373e7eade916261a..8f38066ce46e19ec229282b8071868c1c821689d 100644
--- a/spatio_temporal_dataset/spatial_coordinates/transformed_coordinates.py
+++ b/spatio_temporal_dataset/spatial_coordinates/transformed_coordinates.py
@@ -7,7 +7,7 @@ class TransformedCoordinates(AbstractSpatialCoordinates):
     @classmethod
     def from_coordinates(cls, spatial_coordinates: AbstractSpatialCoordinates,
                          transformation_function: AbstractTransformation):
-        df_coord_transformed = spatial_coordinates.df_coord.copy()
+        df_coord_transformed = spatial_coordinates.df_coordinates.copy()
         df_coord_transformed = transformation_function.transform(df_coord=df_coord_transformed)
         return cls(df_coord=df_coord_transformed, s_split=spatial_coordinates.s_split)
 
diff --git a/spatio_temporal_dataset/temporal_observations/abstract_temporal_observations.py b/spatio_temporal_dataset/temporal_observations/abstract_temporal_observations.py
index c9d7c9e210a032ebf545355dfa43efa6d560b05d..34e1796ab7b7925d9275cbd98369a9b942ce2ab8 100644
--- a/spatio_temporal_dataset/temporal_observations/abstract_temporal_observations.py
+++ b/spatio_temporal_dataset/temporal_observations/abstract_temporal_observations.py
@@ -3,34 +3,34 @@ import pandas as pd
 
 class AbstractTemporalObservations(object):
 
-    def __init__(self, df_maxima_normalized: pd.DataFrame = None, df_maxima: pd.DataFrame = None):
+    def __init__(self, df_maxima_frech: pd.DataFrame = None, df_maxima_gev: pd.DataFrame = None):
         """
         Main attribute of the class is the DataFrame df_maxima
         Index are stations index
         Columns are the temporal moment of the maxima
         """
-        self.df_maxima_normalized = df_maxima_normalized
-        self.df_maxima = df_maxima
+        self.df_maxima_frech = df_maxima_frech
+        self.df_maxima_gev = df_maxima_gev
 
     @classmethod
     def from_df(cls, df):
         pass
 
     @property
-    def maxima(self):
-        return self.df_maxima.values
+    def maxima_gev(self):
+        return self.df_maxima_gev.values
 
     @property
-    def maxima_normalized(self):
-        return self.df_maxima_normalized.values
+    def maxima_frech(self):
+        return self.df_maxima_frech.values
 
-    @maxima_normalized.setter
-    def maxima_normalized(self, maxima_normalized_to_set):
-        assert self.df_maxima_normalized is None
-        assert maxima_normalized_to_set is not None
-        assert maxima_normalized_to_set.shape == self.maxima.shape
-        self.df_maxima_normalized = pd.DataFrame(data=maxima_normalized_to_set, index=self.df_maxima.index,
-                                                 columns=self.df_maxima.columns)
+    @maxima_frech.setter
+    def maxima_frech(self, maxima_frech_to_set):
+        assert maxima_frech_to_set is not None
+        assert maxima_frech_to_set.shape == self.maxima_gev.shape
+        self.df_maxima_frech = pd.DataFrame(data=maxima_frech_to_set,
+                                            index=self.df_maxima_gev.index,
+                                            columns=self.df_maxima_gev.columns)
 
     @property
     def column_to_time_index(self):
@@ -38,14 +38,7 @@ class AbstractTemporalObservations(object):
 
     @property
     def index(self):
-        return self.df_maxima.index
+        return self.df_maxima_gev.index
 
 
-class RealTemporalObservations(object):
 
-    def __init__(self):
-        pass
-
-
-class NormalizedTemporalObservations(object):
-    pass
diff --git a/spatio_temporal_dataset/temporal_observations/annual_maxima_observations.py b/spatio_temporal_dataset/temporal_observations/annual_maxima_observations.py
index 7964e829fa2b9ede4e5cf7071a6c42400563d560..6488ebb1e6e7e02ee36cdafe4a88546d98dc8928 100644
--- a/spatio_temporal_dataset/temporal_observations/annual_maxima_observations.py
+++ b/spatio_temporal_dataset/temporal_observations/annual_maxima_observations.py
@@ -19,9 +19,9 @@ class MarginAnnualMaxima(AnnualMaxima):
     @classmethod
     def from_sampling(cls, nb_obs: int, spatial_coordinates: AbstractSpatialCoordinates,
                       margin_model: AbstractMarginModel):
-        maxima = margin_model.rmargin(nb_obs=nb_obs, coord=spatial_coordinates.coord)
-        df_maxima = pd.DataFrame(data=maxima, index=spatial_coordinates.index)
-        return cls(df_maxima=df_maxima)
+        maxima_gev = margin_model.rmargin(nb_obs=nb_obs, coord=spatial_coordinates.coordinates)
+        df_maxima_gev = pd.DataFrame(data=maxima_gev, index=spatial_coordinates.index)
+        return cls(df_maxima_gev=df_maxima_gev)
 
 
 class MaxStableAnnualMaxima(AbstractTemporalObservations):
@@ -29,9 +29,9 @@ class MaxStableAnnualMaxima(AbstractTemporalObservations):
     @classmethod
     def from_sampling(cls, nb_obs: int, max_stable_model: AbstractMaxStableModel,
                       spatial_coordinates: AbstractSpatialCoordinates):
-        maxima_normalized = max_stable_model.rmaxstab(nb_obs=nb_obs, coord=spatial_coordinates.coord)
-        df_maxima_normalized = pd.DataFrame(data=maxima_normalized, index=spatial_coordinates.index)
-        return cls(df_maxima_normalized=df_maxima_normalized)
+        maxima_frech = max_stable_model.rmaxstab(nb_obs=nb_obs, coord=spatial_coordinates.coordinates)
+        df_maxima_frech = pd.DataFrame(data=maxima_frech, index=spatial_coordinates.index)
+        return cls(df_maxima_frech=df_maxima_frech)
 
 
 class FullAnnualMaxima(MaxStableAnnualMaxima):
@@ -41,7 +41,7 @@ class FullAnnualMaxima(MaxStableAnnualMaxima):
                              spatial_coordinates: AbstractSpatialCoordinates,
                              margin_model: AbstractMarginModel):
         max_stable_annual_maxima = super().from_sampling(nb_obs, max_stable_model, spatial_coordinates)
-        #  Compute df_maxima from df_maxima_normalized
-        maxima = margin_model.get_maxima(max_stable_annual_maxima.maxima_normalized, spatial_coordinates.coord)
-        max_stable_annual_maxima.df_maxima = pd.DataFrame(data=maxima, index=spatial_coordinates.index)
+        #  Compute df_maxima_gev from df_maxima_frech
+        maxima_gev = margin_model.frech2gev(max_stable_annual_maxima.maxima_frech, spatial_coordinates.coordinates)
+        max_stable_annual_maxima.df_maxima_gev = pd.DataFrame(data=maxima_gev, index=spatial_coordinates.index)
         return max_stable_annual_maxima
diff --git a/test/extreme_estimator/test_full_estimators.py b/test/extreme_estimator/test_full_estimators.py
index 075e56bba37df0a1fe1b193a9eb0aa014a7a3096..9ec15296c44efccdd75fa816d365f2a7df7d1bf8 100644
--- a/test/extreme_estimator/test_full_estimators.py
+++ b/test/extreme_estimator/test_full_estimators.py
@@ -2,7 +2,7 @@ import unittest
 
 from extreme_estimator.estimator.full_estimator import FullEstimatorInASingleStep, \
     FullEstimatorInASingleStepWithSmoothMarginals, SmoothMarginalsThenUnitaryMsp
-from spatio_temporal_dataset.dataset.simulation_dataset import MarginDataset
+from spatio_temporal_dataset.dataset.simulation_dataset import MarginDataset, FullSimulatedDataset
 from spatio_temporal_dataset.spatial_coordinates.generated_coordinates import CircleCoordinatesRadius1
 from test.extreme_estimator.test_margin_estimators import TestMarginEstimators
 from test.extreme_estimator.test_max_stable_estimators import TestMaxStableEstimators
@@ -22,8 +22,9 @@ class TestFullEstimators(unittest.TestCase):
     def test_full_estimators(self):
         print(self.margin_models, self.max_stable_models)
         for margin_model, max_stable_model in product(self.margin_models, self.max_stable_models):
-            dataset = MarginDataset.from_sampling(nb_obs=10, margin_model=margin_model,
-                                                  spatial_coordinates=self.spatial_coord)
+            dataset = FullSimulatedDataset.from_double_sampling(nb_obs=10, margin_model=margin_model,
+                                                                spatial_coordinates=self.spatial_coord,
+                                                                max_stable_model=max_stable_model)
 
             for estimator_class in self.FULL_ESTIMATORS:
                 estimator = estimator_class(dataset=dataset, margin_model=margin_model,