diff --git a/experiment/paper_past_snow_loads/result_trends_and_return_levels/eurocode_visualizer.py b/experiment/paper_past_snow_loads/result_trends_and_return_levels/eurocode_visualizer.py
index 6b887c9c8f872cffd49882f66c2e30d8d9614289..78bd644ddaf0967e516c14dec9b6c0da12f2783f 100644
--- a/experiment/paper_past_snow_loads/result_trends_and_return_levels/eurocode_visualizer.py
+++ b/experiment/paper_past_snow_loads/result_trends_and_return_levels/eurocode_visualizer.py
@@ -37,8 +37,8 @@ def plot_uncertainty_massifs(altitude_to_visualizer: Dict[int, StudyVisualizerFo
     :return:
     """
     visualizer = list(altitude_to_visualizer.values())[0]
-    # Subdivide massif names in group of 5
-    m = 4
+    # Subdivide massif names in group of 3
+    m = 3
     uncertainty_massif_names = visualizer.uncertainty_massif_names
     n = (len(uncertainty_massif_names) // m) + 1
     for i in list(range(n))[:]:
@@ -140,6 +140,7 @@ def plot_single_uncertainty_massif_and_non_stationary_context(ax, massif_name, n
     title = '{} ({} Eurocodes area) with a {} model'.format(massif_name_str, eurocode_region_str,
                                                             non_stationary_context)
     ax.set_title(title)
+    ax.set_xticks(altitudes)
     ax.set_ylabel('50-year return level of SL (kN $m^-2$)')
     ax.set_xlabel('Altitude (m)')
     ax.grid()
diff --git a/experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py b/experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py
index 172b03961f99fa5da71179fbea1d566a81082e2d..8554af5c72d3dacc19da037db7b6875ea8484e4d 100644
--- a/experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py
+++ b/experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py
@@ -67,9 +67,16 @@ def major_result():
 
 if __name__ == '__main__':
     # minor_result(altitude=1800)
-    intermediate_result(altitudes=[1500, 1800], massif_names=['Chartreuse'],
-                        uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.ci_mle],
-                        non_stationary_uncertainty=[False])
+    # intermediate_result(altitudes=[1500, 1800], massif_names=['Chartreuse'],
+    #                     uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.ci_mle],
+    #                     non_stationary_uncertainty=[False])
+    # intermediate_result(altitudes=[1500, 1800], massif_names=None,
+    #                     uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.ci_mle],
+    #                     non_stationary_uncertainty=[False])
+    # intermediate_result(altitudes=[300, 600, 900, 1200, 1500, 1800], massif_names=None,
+    #                     uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.ci_mle],
+    #                     non_stationary_uncertainty=[False])
     intermediate_result(altitudes=[300, 600, 900, 1200, 1500, 1800], massif_names=None,
-                        uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.ci_mle],
-                        non_stationary_uncertainty=[False])
+                        uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.ci_mle,
+                                             ConfidenceIntervalMethodFromExtremes.ci_bayes],
+                        non_stationary_uncertainty=[False, True])
diff --git a/experiment/paper_past_snow_loads/result_trends_and_return_levels/study_visualizer_for_non_stationary_trends.py b/experiment/paper_past_snow_loads/result_trends_and_return_levels/study_visualizer_for_non_stationary_trends.py
index e4f210545285ebbbbe98240c045beb00ed245269..b80f46acc75ac728883d4f41cb8edf96c7ab85fe 100644
--- a/experiment/paper_past_snow_loads/result_trends_and_return_levels/study_visualizer_for_non_stationary_trends.py
+++ b/experiment/paper_past_snow_loads/result_trends_and_return_levels/study_visualizer_for_non_stationary_trends.py
@@ -6,6 +6,7 @@ from typing import Dict, Tuple
 import numpy as np
 from cached_property import cached_property
 
+from experiment.eurocode_data.utils import EUROCODE_QUANTILE
 from experiment.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
 from experiment.meteo_france_data.scm_models_data.visualization.study_visualization.study_visualizer import \
     StudyVisualizer
@@ -70,6 +71,12 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
         return {m: np.count_nonzero(maxima) / len(maxima) for m, (_, maxima) in
                 self.massif_name_to_years_and_maxima.items()}
 
+    @cached_property
+    def massif_name_to_eurocode_quantile_in_practice(self):
+        """Due to missing data, the the eurocode quantile which 0.98 if we have all the data
+        correspond in practice to the quantile psnow x 0.98 of the data where there is snow"""
+        return {m: p * EUROCODE_QUANTILE for m, p in self.massif_name_to_psnow.items()}
+
     @cached_property
     def massif_name_to_non_null_years_and_maxima(self):
         d = {}
@@ -80,7 +87,7 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
 
     @cached_property
     def massif_name_to_minimized_aic_non_stationary_trend_test(self) -> Dict[str, AbstractGevTrendTest]:
-        starting_year = 1958
+        starting_year = None
         massif_name_to_trend_test_that_minimized_aic = {}
         for massif_name, (x, y) in self.massif_name_to_non_null_years_and_maxima.items():
             non_stationary_trend_test = [t(x, y, starting_year) for t in self.non_stationary_trend_test]
diff --git a/extreme_fit/model/result_from_model_fit/result_from_extremes/eurocode_return_level_uncertainties.py b/extreme_fit/model/result_from_model_fit/result_from_extremes/eurocode_return_level_uncertainties.py
index 99599dcad12e0d6a85111a6f40a6344069cd2452..71d7e18340600b67e1caf02e61602851e2158d07 100644
--- a/extreme_fit/model/result_from_model_fit/result_from_extremes/eurocode_return_level_uncertainties.py
+++ b/extreme_fit/model/result_from_model_fit/result_from_extremes/eurocode_return_level_uncertainties.py
@@ -49,7 +49,7 @@ class EurocodeConfidenceIntervalFromExtremes(object):
         else:
             fit_method = TemporalMarginFitMethod.extremes_fevd_mle
         # Fitted estimator
-        fitted_estimator = fitted_linear_margin_estimator(model_class, coordinates, dataset, starting_year=1958,
+        fitted_estimator = fitted_linear_margin_estimator(model_class, coordinates, dataset, starting_year=None,
                                                           fit_method=fit_method)
         # Load object from result from extremes
         return cls.from_estimator_extremes(fitted_estimator, ci_method, temporal_covariate)