diff --git a/extreme_fit/model/margin_model/polynomial_margin_model/gev_altitudinal_models.py b/extreme_fit/model/margin_model/polynomial_margin_model/gev_altitudinal_models.py
index dadc236f55c180066251165a99bfe93797bb3881..6007c798e76f5f4c4c381d081d178da9befbcabe 100644
--- a/extreme_fit/model/margin_model/polynomial_margin_model/gev_altitudinal_models.py
+++ b/extreme_fit/model/margin_model/polynomial_margin_model/gev_altitudinal_models.py
@@ -34,23 +34,28 @@ class AbstractAltitudinalModel(AbstractSpatioTemporalPolynomialModel):
 
     @property
     def name_str(self):
-        s = ''
+        l = []
         if self.dim_to_str_number(GevParams.SHAPE, self.coordinates.idx_x_coordinates) == '1':
-            s += '\\zeta_z'
+            l += ['\\zeta_z']
         if self.dim_to_str_number(GevParams.LOC, self.coordinates.idx_temporal_coordinates) in ['1', '2']:
-            s += '\\mu_t'
+            s = '\\mu_t'
             if self.dim_to_str_number(GevParams.LOC, self.coordinates.idx_temporal_coordinates) == '2':
                 s += '^2'
+            l += [s]
         if self.dim_to_str_number(GevParams.SCALE, self.coordinates.idx_temporal_coordinates) in ['1', '2']:
-            s += '\\sigma_t'
+            s = '\\sigma_t'
             if self.dim_to_str_number(GevParams.SCALE, self.coordinates.idx_temporal_coordinates) == '2':
                 s += '^2'
+            l += [s]
         if self.dim_to_str_number(GevParams.SHAPE, self.coordinates.idx_temporal_coordinates) in ['1', '2']:
-            s += '\\zeta_t'
+            s = '\\zeta_t'
             if self.dim_to_str_number(GevParams.SHAPE, self.coordinates.idx_temporal_coordinates) == '2':
                 s += '^2'
-        if len(s) == 0:
+            l += [s]
+        if len(l) == 0:
             s = '0'
+        else:
+            s = ','.join(l)
         return '$ \\boldsymbol{ \\mathcal{M}_{%s} }$' % s
 
 
diff --git a/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py b/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
index 64cb1a6098387332c91b71f39dfec7c6f91e779e..452965cc10d9b1a69938ecc4978fa3e644fe41de 100644
--- a/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
+++ b/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
@@ -31,7 +31,7 @@ def main():
         massif_names = None
         altitudes_list = altitudes_for_groups[:2]
     elif fast:
-        massif_names = ['Vanoise', 'Haute-Tarentaise', 'Vercors']
+        massif_names = ['Vanoise', 'Haute-Maurienne', 'Vercors']
         altitudes_list = altitudes_for_groups[:2]
     else:
         massif_names = None
@@ -57,15 +57,15 @@ def main_loop(altitudes_list, massif_names, seasons, study_classes):
 def plot_visualizers(massif_names, visualizer_list):
     # plot_histogram_all_trends_against_altitudes(massif_names, visualizer_list)
     # plot_histogram_all_models_against_altitudes(massif_names, visualizer_list)
-    for relative in [True, False]:
-        plot_shoe_plot_changes_against_altitude(massif_names, visualizer_list, relative=relative)
+    # for relative in [True, False]:
+    #     plot_shoe_plot_changes_against_altitude(massif_names, visualizer_list, relative=relative)
     # plot_coherence_curves(massif_names, visualizer_list)
     pass
 
 
 def plot_visualizer(massif_names, visualizer):
     # Plot time series
-    # visualizer.studies.plot_maxima_time_series(massif_names=massif_names)
+    visualizer.studies.plot_maxima_time_series(massif_names=massif_names)
     # Plot moments against altitude
     # for std in [True, False][:]:
     #     for change in [True, False, None]:
diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitude_group.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitude_group.py
index fecd357a13dddb89a3505d090c8e43a9bad54934..196c968472f509a04337f15ca0f92551ee2a0036 100644
--- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitude_group.py
+++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitude_group.py
@@ -25,6 +25,11 @@ class AbstractAltitudeGroup(object):
 
     @property
     def xlabel(self):
+        return 'Elevation = {} m. Models are estimated with\n' \
+               'maxima from group {}, i.e. {}'.format(self.reference_altitude, self.group_id, self.formula)
+
+    @property
+    def formula(self):
         # warning: the label could not correspond to all massifs, some might have been fitted with less data
         # idx = get_index_group_from_reference_altitude(reference_altitude)
         # min_altitude, *_, max_altitude = altitudes_for_groups[idx]
@@ -36,8 +41,12 @@ class AbstractAltitudeGroup(object):
         else:
             min_altitude, max_altitude = 1000 * i, 1000 * (i + 1)
             formula = 'between {} m and {} m'.format(min_altitude, max_altitude)
-        return 'Altitude = {} m. Selected models were estimated with\n' \
-               'altitude group {}, i.e. with maxima {}'.format(self.reference_altitude, self.group_id, formula)
+        return formula
+
+    @property
+    def formula_upper(self):
+        f = self.formula
+        return f[0].upper() + f[1:]
 
 
 class LowAltitudeGroup(AbstractAltitudeGroup):
diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
index 8825fb09ee8fc4a3781d935435f3e84282b39bc1..3ff15e6f2b56193c3ba8cef7ef22806406beea19 100644
--- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
+++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
@@ -429,6 +429,8 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
                 if one_fold.is_significant:
                     non_stationary_significant_changes.append(change)
 
+        moment = 'relative mean' if relative else 'Mean'
+        print('{} for {}m'.format(moment, self.altitude_group.reference_altitude), np.mean(changes))
         return changes, non_stationary_changes, non_stationary_significant_changes
 
     def get_valid_names(self, massif_names):
diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py
index 0f4a644ae0730cae9c57a6a871cf3c980a85ac24..b41207fc0d01d6f40e0867241b1b4892bdb4af93 100644
--- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py
+++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py
@@ -313,7 +313,7 @@ class OneFoldFit(object):
         return standard_gumbel_quantiles
 
     def best_confidence_interval(self, altitude, year) -> EurocodeConfidenceIntervalFromExtremes:
-        EurocodeConfidenceIntervalFromExtremes.quantile_level = 1 - (1 / self.return_period)
+        EurocodeConfidenceIntervalFromExtremes.quantile_level = self.quantile_level
         return EurocodeConfidenceIntervalFromExtremes.from_estimator_extremes(estimator_extremes=self.best_estimator,
                                                                               ci_method=ConfidenceIntervalMethodFromExtremes.ci_mle,
                                                                               coordinate=np.array([altitude, year]))
diff --git a/projects/altitude_spatial_model/altitudes_fit/plots/plot_histogram_altitude_studies.py b/projects/altitude_spatial_model/altitudes_fit/plots/plot_histogram_altitude_studies.py
index 8c04dd7457cba190e0573111d0e9c1b72f62bf9e..fcc2b133363bf81a944c2924a1464db51edb410f 100644
--- a/projects/altitude_spatial_model/altitudes_fit/plots/plot_histogram_altitude_studies.py
+++ b/projects/altitude_spatial_model/altitudes_fit/plots/plot_histogram_altitude_studies.py
@@ -100,11 +100,12 @@ def plot_histogram_all_trends_against_altitudes(massif_names, visualizer_list: L
                linewidth=linewidth)
     ax.legend(loc='upper left', prop={'size': size})
     ax.set_ylabel('Percentage of massifs (\%) ', fontsize=legend_fontsize)
-    ax.set_xlabel('Altitudes', fontsize=legend_fontsize)
+    ax.set_xlabel('Elevation', fontsize=legend_fontsize)
     ax.tick_params(axis='both', which='major', labelsize=labelsize)
     ax.set_xticks(x)
     ax.yaxis.grid()
-    ax.set_xticklabels([str(v.altitude_group.reference_altitude) for v in visualizer_list])
+    ax.set_ylim(bottom=0)
+    ax.set_xticklabels([v.altitude_group.formula_upper for v in visualizer_list])
 
     plot_nb_massif_on_upper_axis(ax, labelsize, legend_fontsize, nb_massifs, x)
 
@@ -135,23 +136,22 @@ def plot_shoe_plot_changes_against_altitude(massif_names, visualizer_list: List[
 
 
     x = np.array([4 * width * (i + 1) for i in range(len(nb_massifs))])
-
-    for changes in all_changes:
-        print(changes)
     for j, (changes, label, color) in enumerate(list(zip(all_changes, labels, colors)), -1):
-        bplot = ax.boxplot(list(changes), positions=x + j * width, widths=width, patch_artist=True)
+        positions = x + j * width
+        bplot = ax.boxplot(list(changes), positions=positions, widths=width, patch_artist=True, showmeans=True)
         for patch in bplot['boxes']:
             patch.set_facecolor(color)
 
     custom_lines = [Line2D([0], [0], color=color, lw=4) for color in colors]
-    ax.legend(custom_lines, labels, loc='lower right')
+    loc = 'lower right' if relative else 'upper left'
+    ax.legend(custom_lines, labels, loc=loc)
 
     start = 'Relative changes' if relative else 'Changes'
     unit = '\%' if relative else visualizer.study.variable_unit
-    ax.set_ylabel('{} between 1969 and 2019 in {}-year return levels ({})'.format(start, OneFoldFit.return_period,
+    ax.set_ylabel('{} of {}-year return levels between 1969 and 2019 ({})'.format(start, OneFoldFit.return_period,
                                                                                   unit),
                   fontsize=legend_fontsize)
-    ax.set_xlabel('Altitudes', fontsize=legend_fontsize)
+    ax.set_xlabel('Elevation', fontsize=legend_fontsize)
     ax.tick_params(axis='both', which='major', labelsize=labelsize)
     ax.set_xticks(x)
     ax.yaxis.grid()
@@ -179,4 +179,4 @@ def plot_nb_massif_on_upper_axis(ax, labelsize, legend_fontsize, nb_massifs, x):
     ax_twiny.tick_params(labelsize=labelsize)
     ax_twiny.set_xticklabels(nb_massifs)
     ax_twiny.set_xlim(ax.get_xlim())
-    ax_twiny.set_xlabel('Total number of massifs at each altitude (for the percentage)', fontsize=legend_fontsize)
+    ax_twiny.set_xlabel('Total number of massifs at each elevation (for the percentage)', fontsize=legend_fontsize)
diff --git a/projects/altitude_spatial_model/preliminary_analysis.py b/projects/altitude_spatial_model/preliminary_analysis.py
index b2942820e28021b4a0205dad43fb295a808ada2e..d7de53e3838dc9920bfe28af75978b9d31c6cb56 100644
--- a/projects/altitude_spatial_model/preliminary_analysis.py
+++ b/projects/altitude_spatial_model/preliminary_analysis.py
@@ -89,7 +89,7 @@ class PointwiseGevStudyVisualizer(AltitudesStudies):
             # Plot map of slope for each massif
             visualizer = StudyVisualizer(study=self.study, show=False, save_to_file=True)
             idx = 8 if param_name == GevParams.SHAPE else 1
-            label = 'Altitude gradient for the {}'.format(ylabel[:-idx] + '/100m)')
+            label = 'Elevation gradient for the {}'.format(ylabel[:-idx] + '/100m)')
             gev_param_name_to_graduation = {
                 GevParams.LOC: 0.5,
                 GevParams.SCALE: 0.1,