diff --git a/experiment/meteo_france_data/scm_models_data/abstract_study.py b/experiment/meteo_france_data/scm_models_data/abstract_study.py
index b16ead03bcd1eb286b72bab364937aef29fd470d..8c4071292247cd07419fa85b4445f37479cf6d79 100644
--- a/experiment/meteo_france_data/scm_models_data/abstract_study.py
+++ b/experiment/meteo_france_data/scm_models_data/abstract_study.py
@@ -102,6 +102,10 @@ class AbstractStudy(object):
     def observations_annual_maxima(self) -> AnnualMaxima:
         return AnnualMaxima(df_maxima_gev=pd.DataFrame(self.year_to_annual_maxima, index=self.study_massif_names))
 
+    def annual_maxima_and_years(self, massif_name) -> Tuple[np.ndarray, np.ndarray]:
+        df = self.observations_annual_maxima.df_maxima_gev
+        return df.loc[massif_name].values, np.array(df.columns)
+
     @cached_property
     def year_to_annual_maxima(self) -> OrderedDict:
         # Map each year to an array of size nb_massif
diff --git a/experiment/paper1_steps/1 - non stationary model choice/__init__.py b/experiment/paper1/1 - non stationary model choice/__init__.py
similarity index 100%
rename from experiment/paper1_steps/1 - non stationary model choice/__init__.py
rename to experiment/paper1/1 - non stationary model choice/__init__.py
diff --git a/experiment/paper1_steps/1 - non stationary model choice/main2_choice_to_not_use_starting_years.py b/experiment/paper1/1 - non stationary model choice/main2_choice_to_not_use_starting_years.py
similarity index 94%
rename from experiment/paper1_steps/1 - non stationary model choice/main2_choice_to_not_use_starting_years.py
rename to experiment/paper1/1 - non stationary model choice/main2_choice_to_not_use_starting_years.py
index c43a29e7bf0b3b6c28b316174350081b0de4be05..b5a2f7c5975764a55f7c4d2ec5533ba67690abf3 100644
--- a/experiment/paper1_steps/1 - non stationary model choice/main2_choice_to_not_use_starting_years.py	
+++ b/experiment/paper1/1 - non stationary model choice/main2_choice_to_not_use_starting_years.py	
@@ -2,7 +2,7 @@ import time
 
 from experiment.meteo_france_data.scm_models_data.visualization.hypercube_visualization.altitude_year_hypercube_visualizer import \
     Altitude_Hypercube_Year_Visualizer
-from experiment.paper1_steps.utils import get_full_altitude_visualizer, FULL_ALTITUDES
+from experiment.paper1.utils import get_full_altitude_visualizer, FULL_ALTITUDES
 
 
 def main_fast_spatial_repartition():
diff --git a/experiment/paper1_steps/1 - non stationary model choice/main3_non_stationary_strength_evolution.py b/experiment/paper1/1 - non stationary model choice/main3_non_stationary_strength_evolution.py
similarity index 95%
rename from experiment/paper1_steps/1 - non stationary model choice/main3_non_stationary_strength_evolution.py
rename to experiment/paper1/1 - non stationary model choice/main3_non_stationary_strength_evolution.py
index 29713ba1b35b6ada4e64ae01fb2e2f864a3e3ff5..6b174865b1510aa59cc9b5959600e49243d85efd 100644
--- a/experiment/paper1_steps/1 - non stationary model choice/main3_non_stationary_strength_evolution.py	
+++ b/experiment/paper1/1 - non stationary model choice/main3_non_stationary_strength_evolution.py	
@@ -9,7 +9,7 @@ from experiment.trend_analysis.univariate_test.gev_trend_test_two_parameters imp
 """
 Visualize the 0.99 quantile initial value and its evolution
 """
-from experiment.paper1_steps.utils import get_full_altitude_visualizer, FULL_ALTITUDES
+from experiment.paper1.utils import get_full_altitude_visualizer, FULL_ALTITUDES
 
 
 def main_fast_spatial_risk_evolution():
diff --git a/experiment/paper1_steps/1 - non stationary model choice/main4_common_spatial_altitude_starting_years_impact.py b/experiment/paper1/1 - non stationary model choice/main4_common_spatial_altitude_starting_years_impact.py
similarity index 97%
rename from experiment/paper1_steps/1 - non stationary model choice/main4_common_spatial_altitude_starting_years_impact.py
rename to experiment/paper1/1 - non stationary model choice/main4_common_spatial_altitude_starting_years_impact.py
index ac8a63e3512088748050aa26901277f730633ae5..10d4fe7bd18bd83ea5821ab4d3557aea7f2f7245 100644
--- a/experiment/paper1_steps/1 - non stationary model choice/main4_common_spatial_altitude_starting_years_impact.py	
+++ b/experiment/paper1/1 - non stationary model choice/main4_common_spatial_altitude_starting_years_impact.py	
@@ -9,7 +9,7 @@ from experiment.trend_analysis.univariate_test.gev_trend_test_two_parameters imp
 """
 Visualize the 0.99 quantile initial value and its evolution
 """
-from experiment.paper1_steps.utils import get_full_altitude_visualizer, FULL_ALTITUDES
+from experiment.paper1.utils import get_full_altitude_visualizer, FULL_ALTITUDES
 
 
 def main_fast_spatial_risk_evolution():
diff --git a/experiment/paper1_steps/1 - non stationary model choice/main4_common_spatial_starting_years_impact.py b/experiment/paper1/1 - non stationary model choice/main4_common_spatial_starting_years_impact.py
similarity index 97%
rename from experiment/paper1_steps/1 - non stationary model choice/main4_common_spatial_starting_years_impact.py
rename to experiment/paper1/1 - non stationary model choice/main4_common_spatial_starting_years_impact.py
index e4e0248f537a53b98170222a743d10795139cef1..2973ad0f869ac3fa045b7dbb77244873a9af8fe1 100644
--- a/experiment/paper1_steps/1 - non stationary model choice/main4_common_spatial_starting_years_impact.py	
+++ b/experiment/paper1/1 - non stationary model choice/main4_common_spatial_starting_years_impact.py	
@@ -9,7 +9,7 @@ from experiment.trend_analysis.univariate_test.gev_trend_test_two_parameters imp
 """
 Visualize the 0.99 quantile initial value and its evolution
 """
-from experiment.paper1_steps.utils import get_full_altitude_visualizer, FULL_ALTITUDES
+from experiment.paper1.utils import get_full_altitude_visualizer, FULL_ALTITUDES
 
 
 def main_fast_spatial_risk_evolution():
diff --git a/experiment/paper1_steps/1 - non stationary model choice/main4_individual_starting_years_impact.py b/experiment/paper1/1 - non stationary model choice/main4_individual_starting_years_impact.py
similarity index 96%
rename from experiment/paper1_steps/1 - non stationary model choice/main4_individual_starting_years_impact.py
rename to experiment/paper1/1 - non stationary model choice/main4_individual_starting_years_impact.py
index c8a3fbc8bca4a4d1dca10f84975b4f22247ca9d6..ada1e0731e2342f1d21ef4ad73ecdfc96e81aa5a 100644
--- a/experiment/paper1_steps/1 - non stationary model choice/main4_individual_starting_years_impact.py	
+++ b/experiment/paper1/1 - non stationary model choice/main4_individual_starting_years_impact.py	
@@ -9,7 +9,7 @@ from experiment.trend_analysis.univariate_test.gev_trend_test_two_parameters imp
 """
 Visualize the 0.99 quantile initial value and its evolution
 """
-from experiment.paper1_steps.utils import get_full_altitude_visualizer, FULL_ALTITUDES
+from experiment.paper1.utils import get_full_altitude_visualizer, FULL_ALTITUDES
 
 
 def main_fast_spatial_risk_evolution():
diff --git a/experiment/paper1_steps/__init__.py b/experiment/paper1/__init__.py
similarity index 100%
rename from experiment/paper1_steps/__init__.py
rename to experiment/paper1/__init__.py
diff --git a/experiment/paper1/main_drawing_bayesian.py b/experiment/paper1/main_drawing_bayesian.py
new file mode 100644
index 0000000000000000000000000000000000000000..713147f277497804c630d25d3c2a6928f75db35e
--- /dev/null
+++ b/experiment/paper1/main_drawing_bayesian.py
@@ -0,0 +1,94 @@
+import seaborn as sns
+
+import matplotlib.pyplot as plt
+from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSwe3Days
+from experiment.meteo_france_data.scm_models_data.visualization.utils import create_adjusted_axes
+from experiment.trend_analysis.univariate_test.utils import load_temporal_coordinates_and_dataset, \
+    fitted_linear_margin_estimator
+from extreme_fit.distribution.gev.gev_params import GevParams
+from extreme_fit.model.margin_model.linear_margin_model.abstract_temporal_linear_margin_model import \
+    TemporalMarginFitMethod
+from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import StationaryTemporalModel
+from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_extract_eurocode_return_level import \
+    ExtractEurocodeReturnLevelFromMyBayesianExtremes
+from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import \
+    ConfidenceIntervalMethodFromExtremes
+
+
+def main_drawing_bayesian():
+    return_level_bayesian = get_return_level_bayesian_example()
+    print(return_level_bayesian.result_from_fit.df_all_samples)
+    print(return_level_bayesian.result_from_fit.df_posterior_samples)
+    print(return_level_bayesian.posterior_eurocode_return_level_samples_for_temporal_covariate)
+    axes = create_adjusted_axes(1, 3)
+    # Plot the trajectories on the first axis
+    ax_trajectories = axes[0]
+    ax_trajectories_inverted = ax_trajectories.twinx()
+    df_all_samples = return_level_bayesian.result_from_fit.df_all_samples
+    iteration_step = [i + 1 for i in df_all_samples.index]
+    lns = []
+    # Last color is for the return level
+    colors = ['r', 'g', 'b', 'tab:orange']
+    gev_param_name_to_color = dict(zip(GevParams.PARAM_NAMES, colors))
+    gev_param_name_to_ax = dict(
+        zip(GevParams.PARAM_NAMES, [ax_trajectories, ax_trajectories, ax_trajectories_inverted]))
+    gev_param_name_to_label = {n: GevParams.greek_letter_from_gev_param_name(n) for n in GevParams.PARAM_NAMES}
+    for j, gev_param_name in enumerate(GevParams.PARAM_NAMES[:]):
+        label = gev_param_name_to_label[gev_param_name]
+        ax = gev_param_name_to_ax[gev_param_name]
+        color = gev_param_name_to_color[gev_param_name]
+        ln = ax.plot(iteration_step, df_all_samples.iloc[:, j].values, label=label, color=color)
+        lns.extend(ln)
+    ax_trajectories_inverted.set_ylim(-0.3, 10)
+    ax_trajectories.set_ylim(-0.3, 90)
+    for ax in [ax_trajectories, ax_trajectories_inverted]:
+        ax.set_xlim(min(iteration_step), max(iteration_step))
+    # labs = [l.get_label() for l in lns]
+    # ax_trajectories.legend(lns, labs, loc=0)
+    ax.axvline(x=return_level_bayesian.result_from_fit.burn_in_nb, color='k', linestyle='--')
+    ax_trajectories_inverted.legend(loc=1)
+    ax_trajectories.legend(loc=2)
+    ax_trajectories.set_xlabel("MCMC iterations")
+
+    # Plot the parameter posterior on axes 1
+    ax_parameter_posterior = axes[1]
+    ax_parameter_posterior_flip = ax_parameter_posterior.twiny()
+    gev_param_name_to_ax = dict(
+        zip(GevParams.PARAM_NAMES, [ax_parameter_posterior, ax_parameter_posterior.twiny(), ax_parameter_posterior_flip]))
+    df_posterior_samples = return_level_bayesian.result_from_fit.df_posterior_samples
+    lns = []
+    for j, gev_param_name in enumerate(GevParams.PARAM_NAMES[:]):
+        label = gev_param_name_to_label[gev_param_name]
+        color = gev_param_name_to_color[gev_param_name]
+        ax = gev_param_name_to_ax[gev_param_name]
+        ln = sns.kdeplot(df_posterior_samples.iloc[:, j], ax=ax, label=label, color=color)
+        lns.append(ln)
+    labs = [l.get_label() for l in lns]
+    ax_parameter_posterior.legend(lns, labs, loc=0)
+
+    # Plot the return level posterior on the last axes
+    ax_return_level_posterior = axes[2]
+    sns.kdeplot(return_level_bayesian.posterior_eurocode_return_level_samples_for_temporal_covariate,
+                ax=ax_return_level_posterior, color=colors[-1])
+    ax_return_level_posterior.set_xlabel("$z_p(\\theta)")
+    ax_return_level_posterior.set_ylabel("$p(z_p(\\theta)|y)")
+
+    plt.show()
+
+
+def get_return_level_bayesian_example():
+    maxima, years = CrocusSwe3Days(altitude=1800).annual_maxima_and_years('Vercors')
+    # plt.plot(years, maxima)
+    # plt.show()
+    model_class = StationaryTemporalModel
+    coordinates, dataset = load_temporal_coordinates_and_dataset(maxima, years)
+    fitted_estimator = fitted_linear_margin_estimator(model_class, coordinates, dataset, starting_year=1958,
+                                                      fit_method=TemporalMarginFitMethod.extremes_fevd_bayesian)
+    return_level_bayesian = ExtractEurocodeReturnLevelFromMyBayesianExtremes(estimator=fitted_estimator,
+                                                                             ci_method=ConfidenceIntervalMethodFromExtremes.my_bayes,
+                                                                             temporal_covariate=2017)
+    return return_level_bayesian
+
+
+if __name__ == '__main__':
+    main_drawing_bayesian()
diff --git a/experiment/paper1_steps/poster_EVAN2019/__init__.py b/experiment/paper1/poster_EVAN2019/__init__.py
similarity index 100%
rename from experiment/paper1_steps/poster_EVAN2019/__init__.py
rename to experiment/paper1/poster_EVAN2019/__init__.py
diff --git a/experiment/paper1_steps/poster_EVAN2019/main_poster_EVAN2019.py b/experiment/paper1/poster_EVAN2019/main_poster_EVAN2019.py
similarity index 99%
rename from experiment/paper1_steps/poster_EVAN2019/main_poster_EVAN2019.py
rename to experiment/paper1/poster_EVAN2019/main_poster_EVAN2019.py
index 063783cc2e2ae3f825cf4c08897483565bccdebb..ac9ddb195f475dfdb7e079769cde6ed1662503ec 100644
--- a/experiment/paper1_steps/poster_EVAN2019/main_poster_EVAN2019.py
+++ b/experiment/paper1/poster_EVAN2019/main_poster_EVAN2019.py
@@ -7,7 +7,7 @@ from experiment.meteo_france_data.scm_models_data.visualization.study_visualizat
 from experiment.trend_analysis.univariate_test.gev_trend_test_one_parameter import GevScaleTrendTest, \
     GevLocationTrendTest
 from experiment.trend_analysis.univariate_test.gev_trend_test_two_parameters import GevLocationAndScaleTrendTest
-from experiment.paper1_steps.utils import get_full_altitude_visualizer
+from experiment.paper1.utils import get_full_altitude_visualizer
 
 POSTER_ALTITUDES = [900, 1800, 2700]
 import matplotlib as mpl
diff --git a/experiment/paper1_steps/poster_EVAN2019/shape_prior_check/__init__.py b/experiment/paper1/poster_EVAN2019/shape_prior_check/__init__.py
similarity index 100%
rename from experiment/paper1_steps/poster_EVAN2019/shape_prior_check/__init__.py
rename to experiment/paper1/poster_EVAN2019/shape_prior_check/__init__.py
diff --git a/experiment/paper1_steps/poster_EVAN2019/shape_prior_check/analyse_shape_from_some_experiment.py b/experiment/paper1/poster_EVAN2019/shape_prior_check/analyse_shape_from_some_experiment.py
similarity index 100%
rename from experiment/paper1_steps/poster_EVAN2019/shape_prior_check/analyse_shape_from_some_experiment.py
rename to experiment/paper1/poster_EVAN2019/shape_prior_check/analyse_shape_from_some_experiment.py
diff --git a/experiment/paper1_steps/poster_EVAN2019/shape_prior_check/hist_values_shape.png b/experiment/paper1/poster_EVAN2019/shape_prior_check/hist_values_shape.png
similarity index 100%
rename from experiment/paper1_steps/poster_EVAN2019/shape_prior_check/hist_values_shape.png
rename to experiment/paper1/poster_EVAN2019/shape_prior_check/hist_values_shape.png
diff --git a/experiment/paper1_steps/poster_EVAN2019/shape_prior_check/shape_from_some_experiment.txt b/experiment/paper1/poster_EVAN2019/shape_prior_check/shape_from_some_experiment.txt
similarity index 100%
rename from experiment/paper1_steps/poster_EVAN2019/shape_prior_check/shape_from_some_experiment.txt
rename to experiment/paper1/poster_EVAN2019/shape_prior_check/shape_from_some_experiment.txt
diff --git a/experiment/paper1_steps/poster_EVAN2019/shape_prior_check/some_experiment_EVAN.py b/experiment/paper1/poster_EVAN2019/shape_prior_check/some_experiment_EVAN.py
similarity index 96%
rename from experiment/paper1_steps/poster_EVAN2019/shape_prior_check/some_experiment_EVAN.py
rename to experiment/paper1/poster_EVAN2019/shape_prior_check/some_experiment_EVAN.py
index d93577b43e06cd15fa86ec350ca77e58a2b7e174..9d0175bfb4872dbef79b61d4ec0b5fb9fb72cf87 100644
--- a/experiment/paper1_steps/poster_EVAN2019/shape_prior_check/some_experiment_EVAN.py
+++ b/experiment/paper1/poster_EVAN2019/shape_prior_check/some_experiment_EVAN.py
@@ -9,7 +9,7 @@ from experiment.trend_analysis.univariate_test.abstract_comparison_non_stationar
 from experiment.trend_analysis.univariate_test.gev_trend_test_one_parameter import GevScaleTrendTest, \
     GevLocationTrendTest
 from experiment.trend_analysis.univariate_test.gev_trend_test_two_parameters import GevLocationAndScaleTrendTest
-from experiment.paper1_steps.utils import get_full_altitude_visualizer
+from experiment.paper1.utils import get_full_altitude_visualizer
 
 POSTER_ALTITUDES = [900, 1800, 2700]
 import matplotlib as mpl
diff --git a/experiment/paper1_steps/snow_load_trends_percentage_significant.py b/experiment/paper1/snow_load_trends_percentage_significant.py
similarity index 100%
rename from experiment/paper1_steps/snow_load_trends_percentage_significant.py
rename to experiment/paper1/snow_load_trends_percentage_significant.py
diff --git a/experiment/paper1_steps/utils.py b/experiment/paper1/utils.py
similarity index 100%
rename from experiment/paper1_steps/utils.py
rename to experiment/paper1/utils.py
diff --git a/experiment/paper1_steps/validations/__init__.py b/experiment/paper1/validations/__init__.py
similarity index 100%
rename from experiment/paper1_steps/validations/__init__.py
rename to experiment/paper1/validations/__init__.py
diff --git a/experiment/paper1_steps/validations/main0_comparison_with_observations.py b/experiment/paper1/validations/main0_comparison_with_observations.py
similarity index 100%
rename from experiment/paper1_steps/validations/main0_comparison_with_observations.py
rename to experiment/paper1/validations/main0_comparison_with_observations.py
diff --git a/experiment/paper1_steps/validations/main1_good_stationary_gev_fit.py b/experiment/paper1/validations/main1_good_stationary_gev_fit.py
similarity index 100%
rename from experiment/paper1_steps/validations/main1_good_stationary_gev_fit.py
rename to experiment/paper1/validations/main1_good_stationary_gev_fit.py