diff --git a/extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py b/extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py
index 50a12b74a5d4413f40bac2fdee4a2f5fed1c5cd4..32270574da66638ecd3f14b4042725bb2bdf0e4f 100644
--- a/extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py
+++ b/extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py
@@ -490,7 +490,7 @@ class StudyVisualizer(VisualizationParameters):
         estimator.fit()
 
         # Set visualization attributes for margin_fct
-        margin_fct = estimator.function_from_fit
+        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/extreme_fit/estimator/abstract_estimator.py b/extreme_fit/estimator/abstract_estimator.py
index 68c3bbf9b3f5fc8b291d3a5a8d772aeb8516dcd6..a7057f22c22ccc62a2efe7c30c9505945460d072 100644
--- a/extreme_fit/estimator/abstract_estimator.py
+++ b/extreme_fit/estimator/abstract_estimator.py
@@ -36,7 +36,7 @@ class AbstractEstimator(object):
         return self._result_from_fit
 
     @cached_property
-    def function_from_fit(self) -> AbstractFunction:
+    def margin_function_from_fit(self) -> AbstractFunction:
         raise NotImplementedError
 
 
diff --git a/extreme_fit/estimator/full_estimator/abstract_full_estimator.py b/extreme_fit/estimator/full_estimator/abstract_full_estimator.py
index be66fb1b9a70cb23608c9bb92cc0fa9e62d43d8c..6841026d9948a25910f7f55a64e511178844231d 100644
--- a/extreme_fit/estimator/full_estimator/abstract_full_estimator.py
+++ b/extreme_fit/estimator/full_estimator/abstract_full_estimator.py
@@ -33,7 +33,7 @@ class SmoothMarginalsThenUnitaryMsp(AbstractFullEstimator):
         maxima_gev_train = self.dataset.maxima_gev
         maxima_frech = AbstractMarginModel.gev2frech(maxima_gev=maxima_gev_train,
                                                      coordinates_values=self.dataset.coordinates_values(),
-                                                     margin_function=self.margin_estimator.function_from_fit)
+                                                     margin_function=self.margin_estimator.margin_function_from_fit)
         # Update maxima frech field through the dataset object
         self.dataset.set_maxima_frech(maxima_frech)
         # Estimate the max stable parameters
@@ -84,7 +84,7 @@ class FullEstimatorInASingleStepWithSmoothMargin(AbstractFullEstimator):
         )
 
     @cached_property
-    def function_from_fit(self) -> LinearMarginFunction:
+    def margin_function_from_fit(self) -> LinearMarginFunction:
         return load_margin_function(self, self.linear_margin_model)
 
 
diff --git a/extreme_fit/estimator/margin_estimator/abstract_margin_estimator.py b/extreme_fit/estimator/margin_estimator/abstract_margin_estimator.py
index 97a64880c93cfc371fee1d54050712d4f0b8b548..f6455310da7afc09b4d448aa4f2f7ca74692b7ca 100644
--- a/extreme_fit/estimator/margin_estimator/abstract_margin_estimator.py
+++ b/extreme_fit/estimator/margin_estimator/abstract_margin_estimator.py
@@ -59,7 +59,7 @@ class LinearMarginEstimator(AbstractMarginEstimator):
             climate_coordinates_with_effects=self.margin_model.climate_coordinates_with_effects)
 
     @cached_property
-    def function_from_fit(self) -> LinearMarginFunction:
+    def margin_function_from_fit(self) -> LinearMarginFunction:
         return load_margin_function(self, self.margin_model)
 
     @property
@@ -70,7 +70,7 @@ class LinearMarginEstimator(AbstractMarginEstimator):
     def nllh(self):
         maxima_values = self.dataset.maxima_gev
         coordinate_values = self.coordinates_for_nllh
-        return compute_nllh(coordinate_values, maxima_values, self.function_from_fit)
+        return compute_nllh(coordinate_values, maxima_values, self.margin_function_from_fit)
 
     def sorted_empirical_standard_gumbel_quantiles(self, coordinate_for_filter=None):
         sorted_empirical_quantiles = []
@@ -82,7 +82,7 @@ class LinearMarginEstimator(AbstractMarginEstimator):
                 keep = any([(f is not None) and (c == f) for c, f in zip(coordinate, coordinate_for_filter)])
                 if not keep:
                     continue
-            gev_param = self.function_from_fit.get_params(
+            gev_param = self.margin_function_from_fit.get_params(
                 coordinate=coordinate,
                 is_transformed=False)
             maximum_standardized = gev_param.gumbel_standardization(maximum[0])
@@ -95,7 +95,7 @@ class LinearMarginEstimator(AbstractMarginEstimator):
         coordinate_values = self.dataset.df_coordinates.values
         assert len(standard_gumbel_quantiles) == len(coordinate_values)
         for quantile, coordinate in zip(standard_gumbel_quantiles, coordinate_values):
-            gev_param = self.function_from_fit.get_params(
+            gev_param = self.margin_function_from_fit.get_params(
                 coordinate=coordinate,
                 is_transformed=False)
             maximum = gev_param.gumbel_inverse_standardization(quantile)
@@ -118,8 +118,8 @@ class LinearMarginEstimator(AbstractMarginEstimator):
 
     @property
     def nb_params(self):
-        nb_params = self.function_from_fit.nb_params
-        nb_params += self.function_from_fit.nb_params_for_climate_effects
+        nb_params = self.margin_function_from_fit.nb_params
+        nb_params += self.margin_function_from_fit.nb_params_for_climate_effects
         if isinstance(self.margin_model, AbstractTemporalLinearMarginModel) and self.margin_model.is_gumbel_model:
             nb_params -= 1
         return nb_params
diff --git a/extreme_fit/estimator/margin_estimator/utils.py b/extreme_fit/estimator/margin_estimator/utils.py
index 084f670375f2eb87ac9db7428dad006a3e3c4b6a..24975769d4696a71b931dca2bd163e4ec3079488 100644
--- a/extreme_fit/estimator/margin_estimator/utils.py
+++ b/extreme_fit/estimator/margin_estimator/utils.py
@@ -46,7 +46,7 @@ def _fitted_stationary_gev(fit_method, model_class, starting_year, x_gev):
     dataset = AbstractDataset(observations=observations, coordinates=coordinates)
     estimator = fitted_linear_margin_estimator(model_class, coordinates, dataset, starting_year, fit_method)
     first_coordinate = coordinates.coordinates_values()[0]
-    gev_param = estimator.function_from_fit.get_params(first_coordinate)
+    gev_param = estimator.margin_function_from_fit.get_params(first_coordinate)
     # Warning
     if not -0.5 < gev_param.shape < 0.5:
         pass
diff --git a/extreme_fit/estimator/quantile_estimator/abstract_quantile_estimator.py b/extreme_fit/estimator/quantile_estimator/abstract_quantile_estimator.py
index c16b79c355aac1aac7ed72a55722d3e21f3c1ede..102da1394e1b0fc7fa0bbb4627a6514f9b5ab6c2 100644
--- a/extreme_fit/estimator/quantile_estimator/abstract_quantile_estimator.py
+++ b/extreme_fit/estimator/quantile_estimator/abstract_quantile_estimator.py
@@ -15,7 +15,7 @@ class AbstractQuantileEstimator(AbstractEstimator, ABC):
         self.quantile = quantile
 
     @cached_property
-    def function_from_fit(self) -> AbstractQuantileFunction:
+    def margin_function_from_fit(self) -> AbstractQuantileFunction:
         raise NotImplementedError
 
 
diff --git a/extreme_fit/estimator/quantile_estimator/quantile_estimator_from_margin.py b/extreme_fit/estimator/quantile_estimator/quantile_estimator_from_margin.py
index 68e8518cbb5c92c7a461fdc447b56f45b4592b95..343da2158c5183cfd12e26c421ad42c5fd7e8b79 100644
--- a/extreme_fit/estimator/quantile_estimator/quantile_estimator_from_margin.py
+++ b/extreme_fit/estimator/quantile_estimator/quantile_estimator_from_margin.py
@@ -17,6 +17,6 @@ class QuantileEstimatorFromMargin(LinearMarginEstimator, AbstractQuantileEstimat
         super().__init__(dataset=dataset, quantile=quantile, margin_model=margin_model)
 
     @cached_property
-    def function_from_fit(self) -> AbstractQuantileFunction:
-        linear_margin_function = super().function_from_fit  # type: AbstractMarginFunction
+    def margin_function_from_fit(self) -> AbstractQuantileFunction:
+        linear_margin_function = super().margin_function_from_fit  # type: AbstractMarginFunction
         return QuantileFunctionFromMarginFunction(self.dataset.coordinates, linear_margin_function, self.quantile)
diff --git a/extreme_fit/estimator/quantile_estimator/quantile_estimator_from_regression.py b/extreme_fit/estimator/quantile_estimator/quantile_estimator_from_regression.py
index acd1e56a4e136574ec96222848001d09533e6710..c7fea54c9b166491418699268f4612f2efc30d3a 100644
--- a/extreme_fit/estimator/quantile_estimator/quantile_estimator_from_regression.py
+++ b/extreme_fit/estimator/quantile_estimator/quantile_estimator_from_regression.py
@@ -20,7 +20,7 @@ class QuantileRegressionEstimator(AbstractQuantileEstimator):
         return self.quantile_regression_model.fit()
 
     @cached_property
-    def function_from_fit(self) -> AbstractQuantileFunction:
+    def margin_function_from_fit(self) -> AbstractQuantileFunction:
         result_from_model_fit = self.result_from_model_fit  # type: ResultFromQuantreg
         coefs = result_from_model_fit.coefficients
         nb_coefs = len(coefs)
diff --git a/extreme_trend/one_fold_fit/altitudes_studies_visualizer_for_non_stationary_models.py b/extreme_trend/one_fold_fit/altitudes_studies_visualizer_for_non_stationary_models.py
index c6ed249c62855e4f81daa2a29b81f0e06f435a9a..6b41b09d4fb2cbb69ca4f5b4199a563a9486c62f 100644
--- a/extreme_trend/one_fold_fit/altitudes_studies_visualizer_for_non_stationary_models.py
+++ b/extreme_trend/one_fold_fit/altitudes_studies_visualizer_for_non_stationary_models.py
@@ -328,7 +328,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
                     y_list = []
                     for t in t_list:
                         coordinate = np.array([altitude, t])
-                        gev_params = one_fold_fit.best_function_from_fit.get_params(coordinate, is_transformed=False)
+                        gev_params = one_fold_fit.best_margin_function_from_fit.get_params(coordinate, is_transformed=False)
                         if plot_mean:
                             y = gev_params.mean
                         else:
@@ -378,7 +378,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
         step = 10
         altitudes = list(np.arange(min(altitudes), max(altitudes) + step, step))
         # Get all the correspond peak years
-        margin_function = self.massif_name_to_one_fold_fit[massif_name].best_function_from_fit
+        margin_function = self.massif_name_to_one_fold_fit[massif_name].best_margin_function_from_fit
         peak_years = []
         year_left = 1900
         switch_altitudes = []
diff --git a/extreme_trend/one_fold_fit/one_fold_fit.py b/extreme_trend/one_fold_fit/one_fold_fit.py
index 0ae4806a7a6f24eadf696a861b7db87aed182027..b061dd407db8ce462298135af694ba68c571615c 100644
--- a/extreme_trend/one_fold_fit/one_fold_fit.py
+++ b/extreme_trend/one_fold_fit/one_fold_fit.py
@@ -110,7 +110,7 @@ class OneFoldFit(object):
 
     def get_gev_params(self, altitude, year):
         coordinate = self.get_coordinate(altitude, year)
-        gev_params = self.best_function_from_fit.get_params(coordinate, is_transformed=False)
+        gev_params = self.best_margin_function_from_fit.get_params(coordinate, is_transformed=False)
         return gev_params
 
     def moment(self, altitudes, order=1):
@@ -199,7 +199,7 @@ class OneFoldFit(object):
                 coordinate_values_to_check = list(coordinate_values_for_the_fit) + coordinate_values_for_the_result
                 has_undefined_parameters = False
                 for coordinate in coordinate_values_to_check:
-                    gev_params = e.function_from_fit.get_params(coordinate)
+                    gev_params = e.margin_function_from_fit.get_params(coordinate)
                     if gev_params.has_undefined_parameters:
                         has_undefined_parameters = True
                         break
@@ -218,6 +218,12 @@ class OneFoldFit(object):
             print('Error for')
             print(self.massif_name, self.altitude_group)
             raise
+
+        #  Apply the goodness of fit
+        if self.only_models_that_pass_goodness_of_fit_test:
+            return [e for e in sorted_estimators if self.goodness_of_fit_test(e)]
+        if not (self.remove_physically_implausible_models or self.only_models_that_pass_goodness_of_fit_test):
+            assert len(sorted_estimators) == len(self.models_classes)
         return sorted_estimators
 
     def get_coordinate(self, altitude, year):
@@ -232,24 +238,13 @@ class OneFoldFit(object):
 
     def _compute_shape_for_reference_altitude(self, estimator):
         coordinate = self.get_coordinate(self.altitude_plot, self.covariate_after)
-        gev_params = estimator.function_from_fit.get_params(coordinate, is_transformed=False)
+        gev_params = estimator.margin_function_from_fit.get_params(coordinate, is_transformed=False)
         shape = gev_params.shape
         return shape
 
-    @cached_property
-    def sorted_estimators_with_stationary(self):
-        if self.only_models_that_pass_goodness_of_fit_test:
-            return [e for e in self.sorted_estimators
-                    if self.goodness_of_fit_test(e)
-                    ]
-        else:
-            if not self.remove_physically_implausible_models:
-                assert len(self.sorted_estimators) == len(self.models_classes)
-            return self.sorted_estimators
-
     @property
     def has_at_least_one_valid_model(self):
-        return len(self.sorted_estimators_with_stationary) > 0
+        return len(self.sorted_estimators) > 0
 
     @property
     def model_class_to_estimator_with_finite_aic(self):
@@ -258,7 +253,7 @@ class OneFoldFit(object):
     @property
     def best_estimator(self):
         if self.has_at_least_one_valid_model:
-            best_estimator = self.sorted_estimators_with_stationary[0]
+            best_estimator = self.sorted_estimators[0]
             return best_estimator
         else:
             raise ValueError('This object should not have been called because '
@@ -269,8 +264,8 @@ class OneFoldFit(object):
         return self.best_estimator.margin_model
 
     @property
-    def best_function_from_fit(self):
-        return self.best_estimator.function_from_fit
+    def best_margin_function_from_fit(self):
+        return self.best_estimator.margin_function_from_fit
 
     @property
     def best_shape(self):
@@ -282,7 +277,7 @@ class OneFoldFit(object):
 
     def best_coef(self, param_name, dim, degree):
         try:
-            coef = self.best_function_from_fit.param_name_to_coef[param_name]  # type: PolynomialAllCoef
+            coef = self.best_margin_function_from_fit.param_name_to_coef[param_name]  # type: PolynomialAllCoef
             if coef.dim_to_polynomial_coef is None:
                 return coef.intercept
             else:
@@ -403,7 +398,7 @@ class OneFoldFit(object):
 
         # First result - Compute the significance
         sign_of_changes = [self.sign_of_change(f) for f in bootstrap_fitted_functions]
-        if self.sign_of_change(self.best_function_from_fit) > 0:
+        if self.sign_of_change(self.best_margin_function_from_fit) > 0:
             is_significant = np.quantile(sign_of_changes, self.SIGNIFICANCE_LEVEL) > 0
         else:
             is_significant = np.quantile(sign_of_changes, 1 - self.SIGNIFICANCE_LEVEL) < 0
@@ -417,7 +412,7 @@ class OneFoldFit(object):
             for altitude in altitudes:
                 key = (altitude, year)
                 coordinate = self.get_coordinate(altitude, year)
-                mean_estimate = self.get_return_level(self.best_function_from_fit, coordinate)
+                mean_estimate = self.get_return_level(self.best_margin_function_from_fit, coordinate)
                 bootstrap_return_levels = [self.get_return_level(f, coordinate) for f in
                                            bootstrap_fitted_functions]
                 confidence_interval = tuple([np.quantile(bootstrap_return_levels, q)
@@ -437,7 +432,7 @@ class OneFoldFit(object):
         print('last temporal coordinate', last_temporal_coordinate)
         altitude = self.altitude_group.reference_altitude
         coordinate = self.get_coordinate(altitude, last_temporal_coordinate)
-        return self.get_return_level(self.best_function_from_fit, coordinate)
+        return self.get_return_level(self.best_margin_function_from_fit, coordinate)
 
     @property
     def bootstrap_fitted_functions_from_fit(self):
@@ -487,6 +482,6 @@ class OneFoldFit(object):
         dataset = AbstractDataset(observations=observations, coordinates=coordinates)
         model_class = type(self.best_margin_model)
 
-        function_from_fit = self.fitted_linear_margin_estimator(model_class, dataset).function_from_fit
+        function_from_fit = self.fitted_linear_margin_estimator(model_class, dataset).margin_function_from_fit
 
         return function_from_fit
diff --git a/extreme_trend/trend_test/abstract_gev_trend_test.py b/extreme_trend/trend_test/abstract_gev_trend_test.py
index e8c147a57acede8f6dd8740b62873265ef1e5cc9..b60b2d8c7cc8e37d9b4fcbc56eee442bca47f2b2 100644
--- a/extreme_trend/trend_test/abstract_gev_trend_test.py
+++ b/extreme_trend/trend_test/abstract_gev_trend_test.py
@@ -123,8 +123,8 @@ class AbstractGevTrendTest(object):
     # Evolution of the GEV parameters and corresponding quantiles
 
     def get_non_stationary_linear_coef(self, param_name: str):
-        return self.unconstrained_estimator.function_from_fit.get_coef(param_name,
-                                                                       AbstractCoordinates.COORDINATE_T)
+        return self.unconstrained_estimator.margin_function_from_fit.get_coef(param_name,
+                                                                              AbstractCoordinates.COORDINATE_T)
 
     @cached_property
     def unconstrained_estimator_gev_params_last_year(self) -> GevParams:
@@ -317,13 +317,13 @@ class AbstractGevTrendTest(object):
         ax.plot(extended_quantiles, extended_maxima, linestyle='-', label=label, color=color, linewidth=5)
 
     def get_unconstrained_gev_params(self, year):
-        gev_params_year = self.unconstrained_estimator.function_from_fit.get_params(
+        gev_params_year = self.unconstrained_estimator.margin_function_from_fit.get_params(
             coordinate=np.array([year]),
             is_transformed=False)
         return gev_params_year
 
     def get_unconstrained_time_derivative_gev_params(self, year):
-        gev_params_year = self.unconstrained_estimator.function_from_fit.get_first_derivative_param(
+        gev_params_year = self.unconstrained_estimator.margin_function_from_fit.get_first_derivative_param(
             coordinate=np.array([year]),
             is_transformed=False,
             dim=0)
@@ -399,7 +399,7 @@ class AbstractGevTrendTest(object):
         empirical_quantiles = []
         # for year, maximum in sorted(zip(self.years, self.maxima), key=lambda t: t[1]):
         for year, maximum in zip(self.years, self.maxima):
-            gev_param = estimator.function_from_fit.get_params(
+            gev_param = estimator.margin_function_from_fit.get_params(
                 coordinate=np.array([year]),
                 is_transformed=False)
             maximum_standardized = gev_param.gumbel_standardization(maximum)
@@ -446,7 +446,7 @@ class AbstractGevTrendTest(object):
         plt.gca().set_ylim(bottom=0)
 
     def get_gev_params_with_big_shape_and_correct_shape(self):
-        gev_params = self.unconstrained_estimator.function_from_fit.get_params(
+        gev_params = self.unconstrained_estimator.margin_function_from_fit.get_params(
             coordinate=np.array([YEAR_OF_INTEREST_FOR_RETURN_LEVEL]),
             is_transformed=False)  # type: GevParams
         gev_params_with_corrected_shape = GevParams(loc=gev_params.location,
diff --git a/projects/archive/projected_swe/old_weight_computer/non_stationary_weight_computer.py b/projects/archive/projected_swe/old_weight_computer/non_stationary_weight_computer.py
index 13844a737d2659eb139429c92d7e7a4ad2fe41b1..9cdc2709427c5dde455fca6cb8e78f11088532c7 100644
--- a/projects/archive/projected_swe/old_weight_computer/non_stationary_weight_computer.py
+++ b/projects/archive/projected_swe/old_weight_computer/non_stationary_weight_computer.py
@@ -25,4 +25,4 @@ class NllhWeightComputer(AbstractWeightComputer):
     def get_function_from_fit(self, ensemble_fit, massif_name, gcm_rcm_couple):
         visualizer = ensemble_fit.gcm_rcm_couple_to_visualizer[gcm_rcm_couple]
         one_fold_fit = visualizer.massif_name_to_one_fold_fit[massif_name]
-        return one_fold_fit.best_function_from_fit
+        return one_fold_fit.best_margin_function_from_fit
diff --git a/projects/archive/quantile_regression_vs_evt/annual_maxima_simulation/abstract_annual_maxima_simulation.py b/projects/archive/quantile_regression_vs_evt/annual_maxima_simulation/abstract_annual_maxima_simulation.py
index a5fe0bd353d999b9885c9c840e9e629142f19193..f2a8c9858c1f0b7af8ce6a4682daafd8030e7b4a 100644
--- a/projects/archive/quantile_regression_vs_evt/annual_maxima_simulation/abstract_annual_maxima_simulation.py
+++ b/projects/archive/quantile_regression_vs_evt/annual_maxima_simulation/abstract_annual_maxima_simulation.py
@@ -50,5 +50,5 @@ class AnnualMaximaSimulation(AbstractSimulation):
         true_gev_params = margin_model.margin_function.get_params(last_coordinate)
         true_quantile = true_gev_params.quantile(self.quantile_data)
         # Compute estimated values
-        estimated_quantiles = [estimator.function_from_fit.get_quantile(last_coordinate) for estimator in estimators]
+        estimated_quantiles = [estimator.margin_function_from_fit.get_quantile(last_coordinate) for estimator in estimators]
         return 100 * np.abs(np.array(estimated_quantiles) - true_quantile) / true_quantile
diff --git a/projects/projected_extreme_snowfall/results/part_3/plot_gcm_rcm_effects.py b/projects/projected_extreme_snowfall/results/part_3/plot_gcm_rcm_effects.py
index f54934f73dd899e5443c6a76bac459bc980de15d..49cf5db5a4ef4e29700d8dd1e5607e77a74e3672 100644
--- a/projects/projected_extreme_snowfall/results/part_3/plot_gcm_rcm_effects.py
+++ b/projects/projected_extreme_snowfall/results/part_3/plot_gcm_rcm_effects.py
@@ -60,7 +60,7 @@ def plot_curve_gcm_rcm_effect(ax, massif_name, visualizer_list: List[AltitudesSt
         indices = one_fold_fit.dataset.coordinates.get_indices_for_effects(climate_coordinates_with_effects,
                                                                            gcm_rcm_couple)
         assert len(indices) <= 2, indices
-        ordered_climate_effects = one_fold_fit.best_function_from_fit.param_name_to_ordered_climate_effects[param_name]
+        ordered_climate_effects = one_fold_fit.best_margin_function_from_fit.param_name_to_ordered_climate_effects[param_name]
         sum_effects = sum([ordered_climate_effects[i] for i in indices])
         effects.append(sum_effects)
     if len(gcm_rcm_couple) == 2:
diff --git a/projects/projected_extreme_snowfall/simulation/spline_simulation.py b/projects/projected_extreme_snowfall/simulation/spline_simulation.py
index 1b7c6540161b2a51caec474a89eab9e236b4b0b7..a96ee784c466f29988fc2d315e9bc8c167d28202 100644
--- a/projects/projected_extreme_snowfall/simulation/spline_simulation.py
+++ b/projects/projected_extreme_snowfall/simulation/spline_simulation.py
@@ -49,7 +49,7 @@ def get_margin_function_from_fit(model_class, gev_parameter):
                                                coordinates, dataset,
                                                starting_year=None,
                                                fit_method=MarginFitMethod.evgam)
-    return estimator.function_from_fit
+    return estimator.margin_function_from_fit
 
 
 def get_params_from_margin(margin_function, gev_parameter, x):
diff --git a/test/test_extreme_fit/test_estimator/test_margin_estimators.py b/test/test_extreme_fit/test_estimator/test_margin_estimators.py
index 2d998f99c5df07b82f29f34d23d3cc06c7f5eccf..c9e546fa47417062720418fa8597c68e70f52909 100644
--- a/test/test_extreme_fit/test_estimator/test_margin_estimators.py
+++ b/test/test_extreme_fit/test_estimator/test_margin_estimators.py
@@ -34,7 +34,7 @@ class TestSmoothMarginEstimator(unittest.TestCase):
                 # Plot
                 if self.DISPLAY:
                     margin_model.margin_function.visualize_function(show=True)
-                    estimator.function_from_fit.visualize_function(show=True)
+                    estimator.margin_function_from_fit.visualize_function(show=True)
         self.assertTrue(True)
 
 
diff --git a/test/test_extreme_fit/test_estimator/test_quantile_estimator.py b/test/test_extreme_fit/test_estimator/test_quantile_estimator.py
index c3fc11ae74def25ed6fff4e910d36faf938a25ee..6f60b89791f13a4fd33ea0e4d3d52e8a7c2f8ca9 100644
--- a/test/test_extreme_fit/test_estimator/test_quantile_estimator.py
+++ b/test/test_extreme_fit/test_estimator/test_quantile_estimator.py
@@ -33,7 +33,7 @@ class TestQuantileEstimator(unittest.TestCase):
             # Fit quantile estimators
             for quantile_estimator in quantile_estimators:
                 quantile_estimator.fit()
-                quantile_estimator.function_from_fit.visualize(show=self.DISPLAY)
+                quantile_estimator.margin_function_from_fit.visualize(show=self.DISPLAY)
 
         self.assertTrue(True)
 
diff --git a/test/test_extreme_fit/test_estimator/test_spatio_temporal_estimator/test_gev_spatio_temporal_polynomial_extremes_mle.py b/test/test_extreme_fit/test_estimator/test_spatio_temporal_estimator/test_gev_spatio_temporal_polynomial_extremes_mle.py
index 1725c32d124c7499bb9fc442c630252ceb427822..34b6e79386e675f5e0634679e9924c64b9ff14e1 100644
--- a/test/test_extreme_fit/test_estimator/test_spatio_temporal_estimator/test_gev_spatio_temporal_polynomial_extremes_mle.py
+++ b/test/test_extreme_fit/test_estimator/test_spatio_temporal_estimator/test_gev_spatio_temporal_polynomial_extremes_mle.py
@@ -73,7 +73,7 @@ class TestGevTemporalQuadraticExtremesMle(unittest.TestCase):
             confidence_interval = EurocodeConfidenceIntervalFromExtremes.from_estimator_extremes(estimator,
                                                                                                  ci_method=ConfidenceIntervalMethodFromExtremes.ci_mle,
                                                                                                  coordinate=coordinate)
-            gev_params = estimator.function_from_fit.get_params(coordinate)
+            gev_params = estimator.margin_function_from_fit.get_params(coordinate)
             return_level = gev_params.return_level(return_period=50)
             if np.isnan(return_level) or np.isnan(confidence_interval.mean_estimate):
                 self.assertTrue(np.isnan(return_level) and np.isnan(confidence_interval.mean_estimate))
diff --git a/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_evgam.py b/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_evgam.py
index 9952646a5d9f48cc08dd5c782e09964e6e9cc673..0d5d287cee25d1bc70347f8bf9a4fa599e41c705 100644
--- a/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_evgam.py
+++ b/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_evgam.py
@@ -45,7 +45,7 @@ class TestGevTemporalEvGam(unittest.TestCase):
                                                    fit_method=self.fit_method)
         ref = {'loc': 0.043852626609682636, 'scale': 2.0696383929640807, 'shape': 0.8290699088428524}
         for year in range(1, 3):
-            mle_params_estimated = estimator.function_from_fit.get_params(np.array([year])).to_dict()
+            mle_params_estimated = estimator.margin_function_from_fit.get_params(np.array([year])).to_dict()
             for key in ref.keys():
                 self.assertAlmostEqual(ref[key], mle_params_estimated[key], places=3)
 
@@ -56,21 +56,21 @@ class TestGevTemporalEvGam(unittest.TestCase):
             estimator = LinearMarginEstimator(self.dataset, margin_model)
             estimator.fit()
             # Checks that parameters returned are indeed different
-            mle_params_estimated_year1 = estimator.function_from_fit.get_params(np.array([1])).to_dict()
-            mle_params_estimated_year3 = estimator.function_from_fit.get_params(np.array([3])).to_dict()
+            mle_params_estimated_year1 = estimator.margin_function_from_fit.get_params(np.array([1])).to_dict()
+            mle_params_estimated_year3 = estimator.margin_function_from_fit.get_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.function_from_fit.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.function_from_fit.starting_point)
+        self.assertEqual(3, estimator.margin_function_from_fit.starting_point)
         # Checks that parameters returned are indeed different
-        mle_params_estimated_year1 = estimator.function_from_fit.get_params(np.array([1])).to_dict()
-        mle_params_estimated_year3 = estimator.function_from_fit.get_params(np.array([3])).to_dict()
+        mle_params_estimated_year1 = estimator.margin_function_from_fit.get_params(np.array([1])).to_dict()
+        mle_params_estimated_year3 = estimator.margin_function_from_fit.get_params(np.array([3])).to_dict()
         self.assertEqual(mle_params_estimated_year1, mle_params_estimated_year3)
-        mle_params_estimated_year5 = estimator.function_from_fit.get_params(np.array([5])).to_dict()
+        mle_params_estimated_year5 = estimator.margin_function_from_fit.get_params(np.array([5])).to_dict()
         self.assertNotEqual(mle_params_estimated_year5, mle_params_estimated_year3)
 
     def fit_non_stationary_estimator(self, starting_point):
@@ -84,8 +84,8 @@ class TestGevTemporalEvGam(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.function_from_fit.mu1_temporal_trend
-        mu1_estimator2 = estimator2.function_from_fit.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_fit/test_estimator/test_temporal_estimator/test_gev_temporal_extremes_bayesian.py b/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_extremes_bayesian.py
deleted file mode 100644
index 73023c28726c57676f038261a479a11ca0987d79..0000000000000000000000000000000000000000
--- a/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_extremes_bayesian.py
+++ /dev/null
@@ -1,76 +0,0 @@
-import unittest
-
-import numpy as np
-import pandas as pd
-
-from extreme_trend.trend_test.abstract_gev_trend_test import fitted_linear_margin_estimator
-from extreme_fit.model.margin_model.utils import \
-    MarginFitMethod
-from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import StationaryTemporalModel, \
-    NonStationaryLocationTemporalModel, NonStationaryLocationAndScaleTemporalModel
-from extreme_fit.model.utils import r, set_seed_r
-from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
-from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_coordinates import \
-    AbstractTemporalCoordinates
-from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
-from spatio_temporal_dataset.spatio_temporal_observations.abstract_spatio_temporal_observations import \
-    AbstractSpatioTemporalObservations
-
-
-# class TestGevTemporalExtremesBayesian(unittest.TestCase):
-#
-#     def setUp(self) -> None:
-#         set_seed_r()
-#         r("""
-#         N <- 50
-#         loc = 0; scale = 1; shape <- 1
-#         x_gev <- rgev(N, loc = loc, scale = scale, shape = shape)
-#         start_loc = 0; start_scale = 1; start_shape = 1
-#         """)
-#         # Compute the stationary temporal margin with isMev
-#         self.start_year = 0
-#         df = pd.DataFrame({AbstractCoordinates.COORDINATE_T: range(self.start_year, self.start_year + 50)})
-#         self.coordinates = AbstractTemporalCoordinates.from_df(df)
-#         df2 = pd.DataFrame(data=np.array(r['x_gev']), index=df.index)
-#         observations = AbstractSpatioTemporalObservations(df_maxima_gev=df2)
-#         self.dataset = AbstractDataset(observations=observations, coordinates=self.coordinates)
-#         self.fit_method = MarginFitMethod.extremes_fevd_bayesian
-#
-#     def test_gev_temporal_margin_fit_stationary(self):
-#         # Create estimator
-#         estimator = fitted_linear_margin_estimator(StationaryTemporalModel, self.coordinates, self.dataset,
-#                                                           starting_year=0,
-#                                                           fit_method=self.fit_method)
-#         ref = {'loc': 0.34272436381693616, 'scale': 1.3222588712831973, 'shape': 0.30491484962825105}
-#         for year in range(1, 3):
-#             mle_params_estimated = estimator.function_from_fit.get_params(np.array([year])).to_dict()
-#             for key in ref.keys():
-#                 self.assertAlmostEqual(ref[key], mle_params_estimated[key], places=3)
-#
-#     def test_gev_temporal_margin_fit_non_stationary_location(self):
-#         # Create estimator
-#         estimator = fitted_linear_margin_estimator(NonStationaryLocationTemporalModel, self.coordinates, self.dataset,
-#                                                    starting_year=0,
-#                                                    fit_method=self.fit_method)
-#         mu1_values = estimator.result_from_model_fit.df_posterior_samples.iloc[:, 1]
-#         self.assertTrue((mu1_values != 0).any())
-#         # Checks that parameters returned are indeed different
-#         mle_params_estimated_year1 = estimator.function_from_fit.get_params(np.array([1])).to_dict()
-#         mle_params_estimated_year3 = estimator.function_from_fit.get_params(np.array([3])).to_dict()
-#         self.assertNotEqual(mle_params_estimated_year1, mle_params_estimated_year3)
-#
-#     def test_gev_temporal_margin_fit_non_stationary_location_and_scale(self):
-#         # Create estimator
-#         estimator = fitted_linear_margin_estimator(NonStationaryLocationAndScaleTemporalModel, self.coordinates, self.dataset,
-#                                                    starting_year=0,
-#                                                    fit_method=self.fit_method)
-#         mu1_values = estimator.result_from_model_fit.df_posterior_samples.iloc[:, 1]
-#         self.assertTrue((mu1_values != 0).any())
-#         # Checks that parameters returned are indeed different
-#         mle_params_estimated_year1 = estimator.function_from_fit.get_params(np.array([1])).to_dict()
-#         mle_params_estimated_year3 = estimator.function_from_fit.get_params(np.array([3])).to_dict()
-#         self.assertNotEqual(mle_params_estimated_year1, mle_params_estimated_year3)
-
-
-if __name__ == '__main__':
-    unittest.main()
diff --git a/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_extremes_gumbel.py b/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_extremes_gumbel.py
index 7a2863f31ec06b0f897afb5f3a8d70df6b2f7136..6ea39cfc6aa30be183d8e63676fadbf5dac03a57 100644
--- a/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_extremes_gumbel.py
+++ b/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_extremes_gumbel.py
@@ -42,7 +42,7 @@ class TestGevTemporalExtremesGumbel(unittest.TestCase):
                                                    fit_method=MarginFitMethod.extremes_fevd_mle)
         ref = {'loc': -0.0862185692806497, 'scale': 1.0818465357627252, 'shape': 0}
         for year in range(1, 3):
-            mle_params_estimated = estimator.function_from_fit.get_params(np.array([year])).to_dict()
+            mle_params_estimated = estimator.margin_function_from_fit.get_params(np.array([year])).to_dict()
             for key in ref.keys():
                 self.assertAlmostEqual(ref[key], mle_params_estimated[key], places=3)
 
diff --git a/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_extremes_l_moments.py b/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_extremes_l_moments.py
index 153043706530dbe6d5f6982c0c3ff55574576340..6f78c1fbb081c8b9b6c9b111958bcbc892e1989d 100644
--- a/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_extremes_l_moments.py
+++ b/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_extremes_l_moments.py
@@ -45,7 +45,7 @@ class TestGevTemporalExtremesLMoments(unittest.TestCase):
                                                    fit_method=self.fit_method)
         ref = {'loc': 0.0813843045950251, 'scale': 1.1791830110181365, 'shape': 0.6610403806908737}
         for year in range(1, 3):
-            mle_params_estimated = estimator.function_from_fit.get_params(np.array([year])).to_dict()
+            mle_params_estimated = estimator.margin_function_from_fit.get_params(np.array([year])).to_dict()
             for key in ref.keys():
                 self.assertAlmostEqual(ref[key], mle_params_estimated[key], places=3)
 
diff --git a/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_extremes_mle.py b/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_extremes_mle.py
index ebade9d06d248f048e1878fa4528fa78b78d1721..4f7e3230f94460dabc3913ff2690031dd9a20e90 100644
--- a/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_extremes_mle.py
+++ b/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_extremes_mle.py
@@ -43,7 +43,7 @@ class TestGevTemporalExtremesMle(unittest.TestCase):
                                                    fit_method=self.fit_method)
         ref = {'loc': 0.02191974259369493, 'scale': 1.0347946062900268, 'shape': 0.829052520147379}
         for year in range(1, 3):
-            mle_params_estimated = estimator.function_from_fit.get_params(np.array([year])).to_dict()
+            mle_params_estimated = estimator.margin_function_from_fit.get_params(np.array([year])).to_dict()
             for key in ref.keys():
                 self.assertAlmostEqual(ref[key], mle_params_estimated[key], places=3)
             self.assertAlmostEqual(estimator.result_from_model_fit.nllh, estimator.nllh)
@@ -56,8 +56,8 @@ class TestGevTemporalExtremesMle(unittest.TestCase):
                                                    starting_year=0,
                                                    fit_method=self.fit_method)
         # Checks that parameters returned are indeed different
-        mle_params_estimated_year1 = estimator.function_from_fit.get_params(np.array([1])).to_dict()
-        mle_params_estimated_year3 = estimator.function_from_fit.get_params(np.array([3])).to_dict()
+        mle_params_estimated_year1 = estimator.margin_function_from_fit.get_params(np.array([1])).to_dict()
+        mle_params_estimated_year3 = estimator.margin_function_from_fit.get_params(np.array([3])).to_dict()
         self.assertNotEqual(mle_params_estimated_year1, mle_params_estimated_year3)
         self.assertAlmostEqual(estimator.result_from_model_fit.nllh, estimator.nllh)
         self.assertAlmostEqual(estimator.result_from_model_fit.aic, estimator.aic)
@@ -70,8 +70,8 @@ class TestGevTemporalExtremesMle(unittest.TestCase):
                                                    starting_year=0,
                                                    fit_method=self.fit_method)
         # Checks that parameters returned are indeed different
-        mle_params_estimated_year1 = estimator.function_from_fit.get_params(np.array([1])).to_dict()
-        mle_params_estimated_year3 = estimator.function_from_fit.get_params(np.array([3])).to_dict()
+        mle_params_estimated_year1 = estimator.margin_function_from_fit.get_params(np.array([1])).to_dict()
+        mle_params_estimated_year3 = estimator.margin_function_from_fit.get_params(np.array([3])).to_dict()
         self.assertNotEqual(mle_params_estimated_year1, mle_params_estimated_year3)
         self.assertAlmostEqual(estimator.result_from_model_fit.nllh, estimator.nllh)
         self.assertAlmostEqual(estimator.result_from_model_fit.aic, estimator.aic)
diff --git a/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_is_mev.py b/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_is_mev.py
index 4431b8ee9b1110ef671d77df581c53df75696f48..920807d4f32de91afc39ed958c3dca2bd168d513 100644
--- a/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_is_mev.py
+++ b/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_is_mev.py
@@ -45,7 +45,7 @@ class TestGevTemporal(unittest.TestCase):
                                                    fit_method=self.fit_method)
         ref = {'loc': 0.04309190816463247, 'scale': 2.0688696961628437, 'shape': 0.8291528207825063}
         for year in range(1, 3):
-            mle_params_estimated = estimator.function_from_fit.get_params(np.array([year])).to_dict()
+            mle_params_estimated = estimator.margin_function_from_fit.get_params(np.array([year])).to_dict()
             for key in ref.keys():
                 self.assertAlmostEqual(ref[key], mle_params_estimated[key], places=3)
 
@@ -56,21 +56,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.function_from_fit.get_params(np.array([1])).to_dict()
-            mle_params_estimated_year3 = estimator.function_from_fit.get_params(np.array([3])).to_dict()
+            mle_params_estimated_year1 = estimator.margin_function_from_fit.get_params(np.array([1])).to_dict()
+            mle_params_estimated_year3 = estimator.margin_function_from_fit.get_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.function_from_fit.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.function_from_fit.starting_point)
+        self.assertEqual(3, estimator.margin_function_from_fit.starting_point)
         # Checks that parameters returned are indeed different
-        mle_params_estimated_year1 = estimator.function_from_fit.get_params(np.array([1])).to_dict()
-        mle_params_estimated_year3 = estimator.function_from_fit.get_params(np.array([3])).to_dict()
+        mle_params_estimated_year1 = estimator.margin_function_from_fit.get_params(np.array([1])).to_dict()
+        mle_params_estimated_year3 = estimator.margin_function_from_fit.get_params(np.array([3])).to_dict()
         self.assertEqual(mle_params_estimated_year1, mle_params_estimated_year3)
-        mle_params_estimated_year5 = estimator.function_from_fit.get_params(np.array([5])).to_dict()
+        mle_params_estimated_year5 = estimator.margin_function_from_fit.get_params(np.array([5])).to_dict()
         self.assertNotEqual(mle_params_estimated_year5, mle_params_estimated_year3)
 
     def fit_non_stationary_estimator(self, starting_point):
@@ -84,8 +84,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.function_from_fit.mu1_temporal_trend
-        mu1_estimator2 = estimator2.function_from_fit.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_fit/test_estimator/test_temporal_estimator/test_gev_temporal_polynomial_evgam.py b/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_polynomial_evgam.py
index 2dd680eb50bfe079282ecd708f8eadf4a2cc012d..d7bc6981df15019d76c5c1ea8645ae7ff2879f15 100644
--- a/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_polynomial_evgam.py
+++ b/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_polynomial_evgam.py
@@ -44,9 +44,9 @@ class TestGevTemporalPolynomialEvgam(unittest.TestCase):
                                                    starting_year=0,
                                                    fit_method=self.fit_method)
         # Checks that parameters returned are indeed different
-        mle_params_estimated_year1 = estimator.function_from_fit.get_params(np.array([1])).to_dict()
-        mle_params_estimated_year3 = estimator.function_from_fit.get_params(np.array([21])).to_dict()
-        mle_params_estimated_year5 = estimator.function_from_fit.get_params(np.array([41])).to_dict()
+        mle_params_estimated_year1 = estimator.margin_function_from_fit.get_params(np.array([1])).to_dict()
+        mle_params_estimated_year3 = estimator.margin_function_from_fit.get_params(np.array([21])).to_dict()
+        mle_params_estimated_year5 = estimator.margin_function_from_fit.get_params(np.array([41])).to_dict()
         self.assertNotEqual(mle_params_estimated_year1, mle_params_estimated_year3)
         self.assertNotEqual(mle_params_estimated_year3, mle_params_estimated_year5)
         # Assert the relationship for the location is indeed quadratic
diff --git a/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_polynomial_extremes_mle.py b/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_polynomial_extremes_mle.py
index 82386487164a0517a1bfeab1078e3c78b7ce4c43..3446f7df2cfbc5e779aecf73bef9444d341b006e 100644
--- a/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_polynomial_extremes_mle.py
+++ b/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_polynomial_extremes_mle.py
@@ -44,9 +44,9 @@ class TestGevTemporalQuadraticExtremesMle(unittest.TestCase):
                                                    starting_year=0,
                                                    fit_method=self.fit_method)
         # Checks that parameters returned are indeed different
-        mle_params_estimated_year1 = estimator.function_from_fit.get_params(np.array([1])).to_dict()
-        mle_params_estimated_year3 = estimator.function_from_fit.get_params(np.array([21])).to_dict()
-        mle_params_estimated_year5 = estimator.function_from_fit.get_params(np.array([41])).to_dict()
+        mle_params_estimated_year1 = estimator.margin_function_from_fit.get_params(np.array([1])).to_dict()
+        mle_params_estimated_year3 = estimator.margin_function_from_fit.get_params(np.array([21])).to_dict()
+        mle_params_estimated_year5 = estimator.margin_function_from_fit.get_params(np.array([41])).to_dict()
         self.assertNotEqual(mle_params_estimated_year1, mle_params_estimated_year3)
         self.assertNotEqual(mle_params_estimated_year3, mle_params_estimated_year5)
         # Assert the relationship for the location is indeed quadratic
diff --git a/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_spline.py b/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_spline.py
index 149bc07efd1184f2ddbb85a8d2482da5b7cb118f..3aebd07dd4e0f5cfc9a9117c27447f7e0a458a9d 100644
--- a/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_spline.py
+++ b/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_spline.py
@@ -53,9 +53,9 @@ class TestGevTemporalSpline(unittest.TestCase):
                                                    starting_year=None,
                                                    fit_method=self.fit_method)
         # Checks that parameters returned are indeed different
-        mle_params_estimated_year_first = estimator.function_from_fit.get_params(np.array([self.start_year])).to_dict()
-        mle_params_estimated_year_middle = estimator.function_from_fit.get_params(np.array([self.middle_year])).to_dict()
-        mle_params_estimated_year_last = estimator.function_from_fit.get_params(np.array([self.last_year])).to_dict()
+        mle_params_estimated_year_first = estimator.margin_function_from_fit.get_params(np.array([self.start_year])).to_dict()
+        mle_params_estimated_year_middle = estimator.margin_function_from_fit.get_params(np.array([self.middle_year])).to_dict()
+        mle_params_estimated_year_last = estimator.margin_function_from_fit.get_params(np.array([self.last_year])).to_dict()
         self.assertNotEqual(mle_params_estimated_year_first, mle_params_estimated_year_middle)
         self.assertNotEqual(mle_params_estimated_year_middle, mle_params_estimated_year_last)
         # # Assert the relationship for the location is different between the beginning and the end
@@ -71,7 +71,7 @@ class TestGevTemporalSpline(unittest.TestCase):
         for idx, nb_year in enumerate(range(5)):
             year = self.start_year + nb_year
             gev_params_from_result = estimator.result_from_model_fit.get_gev_params_from_result(idx).to_dict()
-            my_gev_params = estimator.function_from_fit.get_params(np.array([year])).to_dict()
+            my_gev_params = estimator.margin_function_from_fit.get_params(np.array([year])).to_dict()
             for param_name in GevParams.PARAM_NAMES:
                 self.assertAlmostEqual(gev_params_from_result[param_name], my_gev_params[param_name],
                                        msg='for the {} parameter at year={}'.format(param_name, year),
diff --git a/test/test_extreme_fit/test_model/test_margin_temporal.py b/test/test_extreme_fit/test_model/test_margin_temporal.py
index dc1f602b5287583546dd01ecb0de01c7d5ca480e..156c286792d24bfbf329ddef509f7d3927a0e3af 100644
--- a/test/test_extreme_fit/test_model/test_margin_temporal.py
+++ b/test/test_extreme_fit/test_model/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.function_from_fit.get_params(coordinate).to_dict()
+            mle_params_estimated = estimator.margin_function_from_fit.get_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.function_from_fit.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.function_from_fit.get_params(coordinate1).to_dict()
+        mle_params_estimated_year1 = estimator.margin_function_from_fit.get_params(coordinate1).to_dict()
         coordinate3 = np.array([0.0, 0.0, 3])
-        mle_params_estimated_year3 = estimator.function_from_fit.get_params(coordinate3).to_dict()
+        mle_params_estimated_year3 = estimator.margin_function_from_fit.get_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.function_from_fit.mu1_temporal_trend, 0.0)
-        self.assertAlmostEqual(estimator.function_from_fit.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.mu1_temporal_trend,
                                places=3)
         # Checks starting point parameter are well passed
-        self.assertEqual(2, estimator.function_from_fit.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.function_from_fit.get_params(coordinate1).to_dict()
+        mle_params_estimated_year1 = estimator.margin_function_from_fit.get_params(coordinate1).to_dict()
         coordinate2 = np.array([0.0, 0.0, 2])
-        mle_params_estimated_year2 = estimator.function_from_fit.get_params(coordinate2).to_dict()
+        mle_params_estimated_year2 = estimator.margin_function_from_fit.get_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.function_from_fit.get_params(coordinate5).to_dict()
+        mle_params_estimated_year5 = estimator.margin_function_from_fit.get_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.function_from_fit.mu1_temporal_trend
-        mu1_estimator2 = estimator2.function_from_fit.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_fit/test_model/test_margin_temporal_transformed.py b/test/test_extreme_fit/test_model/test_margin_temporal_transformed.py
index c03c419338ce7a0fb0e6eb11db420314ebbe737e..f782b58ab0e15f23ecdfa27184a4a35c2789c36a 100644
--- a/test/test_extreme_fit/test_model/test_margin_temporal_transformed.py
+++ b/test/test_extreme_fit/test_model/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.function_from_fit.get_params(coordinate,
-                                                                          is_transformed=False).to_dict()
+            mle_params_estimated = estimator.margin_function_from_fit.get_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.function_from_fit.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.function_from_fit.get_params(coordinate1,
-                                                                            is_transformed=False).to_dict()
+        mle_params_estimated_year1 = estimator.margin_function_from_fit.get_params(coordinate1,
+                                                                                   is_transformed=False).to_dict()
         coordinate3 = np.array([0.0, 0.0, 3])
-        mle_params_estimated_year3 = estimator.function_from_fit.get_params(coordinate3,
-                                                                            is_transformed=False).to_dict()
+        mle_params_estimated_year3 = estimator.margin_function_from_fit.get_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.function_from_fit.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.function_from_fit.get_params(coordinate1,
-                                                                            is_transformed=False).to_dict()
+        mle_params_estimated_year1 = estimator.margin_function_from_fit.get_params(coordinate1,
+                                                                                   is_transformed=False).to_dict()
         coordinate2 = np.array([0.0, 0.0, 2])
-        mle_params_estimated_year2 = estimator.function_from_fit.get_params(coordinate2,
-                                                                            is_transformed=False).to_dict()
+        mle_params_estimated_year2 = estimator.margin_function_from_fit.get_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.function_from_fit.get_params(coordinate5,
-                                                                            is_transformed=False).to_dict()
+        mle_params_estimated_year5 = estimator.margin_function_from_fit.get_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.function_from_fit.mu1_temporal_trend
-        mu1_estimator2 = estimator2.function_from_fit.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_fit/test_model/test_max_stable_temporal.py b/test/test_extreme_fit/test_model/test_max_stable_temporal.py
index 38591b03839d97c3fc80dd2c176a2dbeb369772f..3e2bbe5129964379a306d37ab8afd3aa5d2fc298 100644
--- a/test/test_extreme_fit/test_model/test_max_stable_temporal.py
+++ b/test/test_extreme_fit/test_model/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.function_from_fit.get_params(coordinate).to_dict()
+            mle_params_estimated = estimator.margin_function_from_fit.get_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.function_from_fit.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.function_from_fit.get_params(coordinate1).to_dict()
+        mle_params_estimated_year1 = estimator.margin_function_from_fit.get_params(coordinate1).to_dict()
         coordinate3 = np.array([0.0, 0.0, 3])
-        mle_params_estimated_year3 = estimator.function_from_fit.get_params(coordinate3).to_dict()
+        mle_params_estimated_year3 = estimator.margin_function_from_fit.get_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.function_from_fit.mu1_temporal_trend, 0.0)
-        self.assertAlmostEqual(estimator.function_from_fit.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.mu1_temporal_trend,
                                places=2)
         # Checks starting point parameter are well passed
-        self.assertEqual(2, estimator.function_from_fit.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.function_from_fit.get_params(coordinate1).to_dict()
+        mle_params_estimated_year1 = estimator.margin_function_from_fit.get_params(coordinate1).to_dict()
         coordinate2 = np.array([0.0, 0.0, 2])
-        mle_params_estimated_year2 = estimator.function_from_fit.get_params(coordinate2).to_dict()
+        mle_params_estimated_year2 = estimator.margin_function_from_fit.get_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.function_from_fit.get_params(coordinate5).to_dict()
+        mle_params_estimated_year5 = estimator.margin_function_from_fit.get_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.function_from_fit.mu1_temporal_trend
-        mu1_estimator2 = estimator2.function_from_fit.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_trend/test_ensemble_fit_multiplie_altitudes.py b/test/test_extreme_trend/test_ensemble_fit_multiplie_altitudes.py
index 6f1c92814fe004986fce64e730c138256fef2d25..318a8c5ae9f44e47cf81edd5fb5191c7fc3787e2 100644
--- a/test/test_extreme_trend/test_ensemble_fit_multiplie_altitudes.py
+++ b/test/test_extreme_trend/test_ensemble_fit_multiplie_altitudes.py
@@ -44,7 +44,7 @@ class TestEnsembleFit(unittest.TestCase):
                                                temporal_covariate_for_fit=temporal_covariate,
                                                only_models_that_pass_goodness_of_fit_test=False)
 
-            _ = ensemble_fit.visualizer.massif_name_to_one_fold_fit[self.massif_names[0]].best_function_from_fit
+            _ = ensemble_fit.visualizer.massif_name_to_one_fold_fit[self.massif_names[0]].best_margin_function_from_fit
         self.assertTrue(True)
 
     # def test_ensembe_fit_with_effect(self):