diff --git a/extreme_fit/model/abstract_model.py b/extreme_fit/model/abstract_model.py
index bfb98d73a76ed3f7c7d609b29a7c099b4fdff10c..71fe7319db2cd738920f0e4ef626d3ff80f422aa 100644
--- a/extreme_fit/model/abstract_model.py
+++ b/extreme_fit/model/abstract_model.py
@@ -1,13 +1,13 @@
 class AbstractModel(object):
 
-    def __init__(self, params_start_fit=None, params_sample=None):
+    def __init__(self, params_sample=None):
         self.default_params = None
-        self.user_params_start_fit = params_start_fit
         self.user_params_sample = params_sample
 
     @property
     def params_start_fit(self) -> dict:
-        return self.default_params.copy()
+        # return self.default_params.copy()
+        return self.params_sample.copy()
 
     @property
     def params_sample(self) -> dict:
diff --git a/extreme_fit/model/margin_model/abstract_margin_model.py b/extreme_fit/model/margin_model/abstract_margin_model.py
index 4d657ef5eee475dfc78fa45db2caa002d451f27c..f3129472da53c7d81a98d07ffb17db0febe73470 100644
--- a/extreme_fit/model/margin_model/abstract_margin_model.py
+++ b/extreme_fit/model/margin_model/abstract_margin_model.py
@@ -19,10 +19,9 @@ class AbstractMarginModel(AbstractModel, ABC):
         -margin_function_start_fit for starting to fit
     """
 
-    def __init__(self, coordinates: AbstractCoordinates,
-                 params_start_fit=None, params_sample=None,
+    def __init__(self, coordinates: AbstractCoordinates, params_sample=None,
                  params_class=GevParams):
-        super().__init__(params_start_fit, params_sample)
+        super().__init__(params_sample)
         assert isinstance(coordinates, AbstractCoordinates), type(coordinates)
         self.coordinates = coordinates
         self.margin_function_sample = None  # type: AbstractMarginFunction
diff --git a/extreme_fit/model/margin_model/linear_margin_model/abstract_temporal_linear_margin_model.py b/extreme_fit/model/margin_model/linear_margin_model/abstract_temporal_linear_margin_model.py
index 4f00dd918c04f7f91bb7fd35b531e3f633a128dc..1005c89de4ea3fb94b5b633c87bd6809a488c22f 100644
--- a/extreme_fit/model/margin_model/linear_margin_model/abstract_temporal_linear_margin_model.py
+++ b/extreme_fit/model/margin_model/linear_margin_model/abstract_temporal_linear_margin_model.py
@@ -20,17 +20,16 @@ class AbstractTemporalLinearMarginModel(LinearMarginModel):
     """Linearity only with respect to the temporal coordinates"""
 
     def __init__(self, coordinates: AbstractCoordinates,
-                 params_start_fit=None,
                  params_sample=None, starting_point=None,
                  fit_method=MarginFitMethod.is_mev_gev_fit,
                  nb_iterations_for_bayesian_fit=5000,
-                 params_start_fit_bayesian=None,
+                 params_initial_fit_bayesian=None,
                  type_for_MLE="GEV",
                  params_class=GevParams):
-        super().__init__(coordinates, params_start_fit, params_sample, starting_point,
+        super().__init__(coordinates, params_sample, starting_point,
                          params_class)
         self.type_for_mle = type_for_MLE
-        self.params_start_fit_bayesian = params_start_fit_bayesian
+        self.params_initial_fit_bayesian = params_initial_fit_bayesian
         self.nb_iterations_for_bayesian_fit = nb_iterations_for_bayesian_fit
         assert isinstance(fit_method, MarginFitMethod), fit_method
         self.fit_method = fit_method
@@ -92,8 +91,8 @@ class AbstractTemporalLinearMarginModel(LinearMarginModel):
 
     def extremes_fevd_bayesian_fit(self, x, df_coordinates_temp) -> AbstractResultFromExtremes:
         r_type_argument_kwargs, y = self.extreme_arguments(df_coordinates_temp)
-        params_start_fit = self.params_start_fit_bayesian if self.params_start_fit_bayesian is not None else {}
-        r_type_argument_kwargs['initial'] = r.list(**params_start_fit)
+        params_initial_fit = self.params_initial_fit_bayesian if self.params_initial_fit_bayesian is not None else {}
+        r_type_argument_kwargs['initial'] = r.list(**params_initial_fit)
         # Assert for any non-stationary model that the shape parameter is constant
         # (because the prior function considers that the last parameter should be the shape)
         assert GevParams.SHAPE not in self.margin_function_start_fit.param_name_to_dims \
diff --git a/extreme_fit/model/margin_model/linear_margin_model/linear_margin_model.py b/extreme_fit/model/margin_model/linear_margin_model/linear_margin_model.py
index bce0d46a162fa075ed85d223ff4bf3febf2bf894..23ecd372f38e214a152eaadd55d3437453f5a618 100644
--- a/extreme_fit/model/margin_model/linear_margin_model/linear_margin_model.py
+++ b/extreme_fit/model/margin_model/linear_margin_model/linear_margin_model.py
@@ -12,7 +12,7 @@ class LinearMarginModel(ParametricMarginModel):
         for param_name, coef_list in param_name_to_coef_list.items():
             for idx, coef in enumerate(coef_list, -1):
                 params[(param_name, idx)] = coef
-        return cls(coordinates, params_sample=params, params_start_fit=params, params_class=params_class, **kwargs)
+        return cls(coordinates, params_sample=params, params_class=params_class, **kwargs)
 
     def load_margin_functions(self, param_name_to_dims=None):
         assert param_name_to_dims is not None, 'LinearMarginModel cannot be used for sampling/fitting \n' \
diff --git a/extreme_fit/model/margin_model/linear_margin_model/temporal_linear_margin_models.py b/extreme_fit/model/margin_model/linear_margin_model/temporal_linear_margin_models.py
index 20a37d346879ba33a58a9d1f9819d49c6acddacb..1f8c1e85187be003fb95527694f9bff71fdc0afc 100644
--- a/extreme_fit/model/margin_model/linear_margin_model/temporal_linear_margin_models.py
+++ b/extreme_fit/model/margin_model/linear_margin_model/temporal_linear_margin_models.py
@@ -70,11 +70,11 @@ class NonStationaryLocationAndScaleTemporalModel(AbstractTemporalLinearMarginMod
 
 class GumbelTemporalModel(StationaryTemporalModel):
 
-    def __init__(self, coordinates: AbstractCoordinates, params_start_fit=None,
+    def __init__(self, coordinates: AbstractCoordinates,
                  params_sample=None, starting_point=None, fit_method=MarginFitMethod.is_mev_gev_fit,
-                 nb_iterations_for_bayesian_fit=5000, params_start_fit_bayesian=None):
-        super().__init__(coordinates, params_start_fit, params_sample, starting_point, fit_method,
-                         nb_iterations_for_bayesian_fit, params_start_fit_bayesian, type_for_MLE="Gumbel")
+                 nb_iterations_for_bayesian_fit=5000, params_initial_fit_bayesian=None):
+        super().__init__(coordinates, params_sample, starting_point, fit_method,
+                         nb_iterations_for_bayesian_fit, params_initial_fit_bayesian, type_for_MLE="Gumbel")
 
 
 class NonStationaryLocationGumbelModel(GumbelTemporalModel, NonStationaryLocationTemporalModel):
diff --git a/extreme_fit/model/margin_model/parametric_margin_model.py b/extreme_fit/model/margin_model/parametric_margin_model.py
index 973dcc807caf661234ee381db2d35e1079d03648..e9af27077a3b7bf2f4677d55d05061da0ff4cf03 100644
--- a/extreme_fit/model/margin_model/parametric_margin_model.py
+++ b/extreme_fit/model/margin_model/parametric_margin_model.py
@@ -16,7 +16,7 @@ from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoo
 
 class ParametricMarginModel(AbstractMarginModel, ABC):
 
-    def __init__(self, coordinates: AbstractCoordinates, params_start_fit=None,
+    def __init__(self, coordinates: AbstractCoordinates,
                  params_sample=None, starting_point=None, params_class=GevParams,
                  fit_method=MarginFitMethod.spatial_extremes_mle):
         """
@@ -27,7 +27,7 @@ class ParametricMarginModel(AbstractMarginModel, ABC):
         self.margin_function_sample = None  # type: ParametricMarginFunction
         self.margin_function_start_fit = None  # type: ParametricMarginFunction
         self.drop_duplicates = True
-        super().__init__(coordinates, params_start_fit, params_sample, params_class)
+        super().__init__(coordinates, params_sample, params_class)
 
     def fitmargin_from_maxima_gev(self, data: np.ndarray, df_coordinates_spat: pd.DataFrame,
                                   df_coordinates_temp: pd.DataFrame) -> ResultFromSpatialExtreme:
diff --git a/extreme_fit/model/margin_model/spline_margin_model.py b/extreme_fit/model/margin_model/spline_margin_model.py
index ec509617a2e2bdf2b3cd0e75857aca64743938a2..cfe0e1b04a3475579c53d77a7c0c0548dfb48867 100644
--- a/extreme_fit/model/margin_model/spline_margin_model.py
+++ b/extreme_fit/model/margin_model/spline_margin_model.py
@@ -11,9 +11,9 @@ from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoo
 
 class SplineMarginModel(ParametricMarginModel):
 
-    def __init__(self, coordinates: AbstractCoordinates, params_start_fit=None,
+    def __init__(self, coordinates: AbstractCoordinates,
                  params_sample=None):
-        super().__init__(coordinates, params_start_fit, params_sample)
+        super().__init__(coordinates, params_sample)
 
     def load_margin_functions(self, param_name_to_dims: Dict[str, List[int]] = None,
                               param_name_to_coef: Dict[str, AbstractCoef] = None,
diff --git a/extreme_fit/model/max_stable_model/abstract_max_stable_model.py b/extreme_fit/model/max_stable_model/abstract_max_stable_model.py
index 6776d8d086eecd3d797cf5fdc94a502ea07b3068..40634790f32cae5c0d16b8498dca4419cf8abf29 100644
--- a/extreme_fit/model/max_stable_model/abstract_max_stable_model.py
+++ b/extreme_fit/model/max_stable_model/abstract_max_stable_model.py
@@ -14,8 +14,8 @@ from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoo
 
 class AbstractMaxStableModel(AbstractModel):
 
-    def __init__(self, params_start_fit=None, params_sample=None):
-        super().__init__(params_start_fit, params_sample)
+    def __init__(self, params_sample=None):
+        super().__init__(params_sample)
         self.cov_mod = None
 
     @property
@@ -105,9 +105,9 @@ class CovarianceFunction(Enum):
 
 class AbstractMaxStableModelWithCovarianceFunction(AbstractMaxStableModel):
 
-    def __init__(self, params_start_fit=None, params_sample=None,
+    def __init__(self, params_sample=None,
                  covariance_function: CovarianceFunction = None):
-        super().__init__(params_start_fit, params_sample)
+        super().__init__(params_sample)
         assert covariance_function is not None
         self.covariance_function = covariance_function
         self.default_params = {
diff --git a/projects/contrasting_trends_in_snow_loads/check_mcmc_convergence_for_return_levels/gelman_convergence_test.py b/projects/contrasting_trends_in_snow_loads/check_mcmc_convergence_for_return_levels/gelman_convergence_test.py
index 5d8b829d47aead3e8c186bdc42161de777baf197..e5f6b7a5575abd544f01cb9cbf21d5693d261266 100644
--- a/projects/contrasting_trends_in_snow_loads/check_mcmc_convergence_for_return_levels/gelman_convergence_test.py
+++ b/projects/contrasting_trends_in_snow_loads/check_mcmc_convergence_for_return_levels/gelman_convergence_test.py
@@ -31,16 +31,16 @@ def compute_refined_gelman_score(means, variances, N, M):
 
 def compute_gelman_convergence_value(non_null_years_and_maxima, mcmc_iterations, model_class, nb_chains):
     s_means, s_variances = [], []
-    df_params_start_fit = sample_starting_value(nb_chains)
-    for i, row in df_params_start_fit.iterrows():
+    df_params_initial_fit_bayesian = sample_starting_value(nb_chains)
+    for i, row in df_params_initial_fit_bayesian.iterrows():
         s_mean, s_variance = compute_mean_and_variance(mcmc_iterations, model_class, non_null_years_and_maxima,
-                                                       params_start_fit=row.to_dict())
+                                                       params_initial_fit_bayesian=row.to_dict())
         s_means.append(s_mean)
         s_variances.append(s_variance)
     df_mean = pd.concat(s_means, axis=1).transpose()
     df_variance = pd.concat(s_variances, axis=1).transpose()
     param_name_to_R = {param_name: compute_gelman_score(df_mean[param_name], df_variance[param_name], N=mcmc_iterations // 2, M=nb_chains)
-                       for param_name in df_params_start_fit.columns}
+                       for param_name in df_params_initial_fit_bayesian.columns}
     print(param_name_to_R)
     return max(param_name_to_R.values())
 
@@ -53,13 +53,13 @@ def sample_starting_value(nb_chains):
     })
 
 
-def compute_mean_and_variance(mcmc_iterations, model_class, non_null_years_and_maxima, params_start_fit):
+def compute_mean_and_variance(mcmc_iterations, model_class, non_null_years_and_maxima, params_initial_fit_bayesian):
     maxima, years = non_null_years_and_maxima
     coordinates, dataset = load_temporal_coordinates_and_dataset(maxima, years)
     fitted_estimator = fitted_linear_margin_estimator(model_class, coordinates, dataset, starting_year=None,
                                                       fit_method=MarginFitMethod.extremes_fevd_bayesian,
                                                       nb_iterations_for_bayesian_fit=mcmc_iterations,
-                                                      params_start_fit_bayesian=params_start_fit)
+                                                      params_initial_fit_bayesian=params_initial_fit_bayesian)
     res = fitted_estimator.result_from_model_fit  # type: ResultFromBayesianExtremes
     return res.mean_posterior_parameters, res.variance_posterior_parameters