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 cc7fe98e7d4548feb9c9f498a5572dccbf3a4fc4..9506fb209f6e9dad3bb0e80219430b7068aa7b05 100644
--- a/experiment/meteo_france_data/scm_models_data/abstract_study.py
+++ b/experiment/meteo_france_data/scm_models_data/abstract_study.py
@@ -406,7 +406,8 @@ class AbstractStudy(object):
             legend_elements = cls.get_legend_for_model_symbol(marker_style_to_label_name, markersize=8)
             ax.legend(handles=legend_elements[:], loc='upper right')
             # ax.legend(handles=legend_elements[4:], bbox_to_anchor=(0.01, 0.03), loc='lower left')
-            ax.annotate("Filled symbol = significant trend ", xy=(0.05, 0.015), xycoords='axes fraction', fontsize=7)
+            ax.annotate("Filled symbol = significant trend w.r.t $\mathcal{M}_0$",
+                        xy=(0.05, 0.015), xycoords='axes fraction', fontsize=7)
 
         if show:
             plt.show()
diff --git a/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/main_study_visualizer.py b/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/main_study_visualizer.py
index 1b9e076e8a86f1b36f6ec7c4381f4cb684dc8992..cacd15f398166f173ac705c7f29bfa464af84670 100644
--- a/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/main_study_visualizer.py
+++ b/experiment/meteo_france_data/scm_models_data/visualization/study_visualization/main_study_visualizer.py
@@ -35,7 +35,7 @@ SCM_STUDY_CLASS_TO_ABBREVIATION = {
     SafranSnowfall: 'SF3',
     CrocusSweTotal: 'SWE',
     CrocusSwe3Days: 'SWE3',
-    CrocusSnowLoadEurocode: 'GSL from annual maximum of HS and {}'.format(eurocode_snow_density),
+    CrocusSnowLoadEurocode: 'GSL from annual maximum of HS \nand {}'.format(eurocode_snow_density),
     CrocusDepth: 'SD',
     CrocusSnowLoadTotal: 'GSL',
     CrocusSnowLoad3Days: 'GSL3',
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 1ee9c7bb673c17648933fa7df935d865b5e8f446..858f8b338090ebd93283543176c6709a31f49b1d 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
@@ -63,12 +63,14 @@ def intermediate_result(altitudes, massif_names=None,
             _ = compute_minimized_aic(visualizer)
     # Compute common max value for the colorbar
     altitudes_for_plot_trend = [900, 1800, 2700]
+    altitudes_for_plot_trend = altitudes
     visualizers_for_altitudes = [visualizer
                                  for altitude, visualizer in altitude_to_visualizer.items()
                                  if altitude in altitudes_for_plot_trend]
     max_abs_tdrl = max([visualizer.max_abs_change for visualizer in visualizers_for_altitudes])
     for visualizer in visualizers_for_altitudes:
-        visualizer.plot_trends(max_abs_tdrl, add_colorbar=visualizer.study.altitude == 2700)
+        # visualizer.plot_trends(max_abs_tdrl, add_colorbar=visualizer.study.altitude == 2700)
+        visualizer.plot_trends(None, add_colorbar=True)
 
     # Plot graph
     plot_uncertainty_massifs(altitude_to_visualizer)
@@ -78,9 +80,9 @@ def intermediate_result(altitudes, massif_names=None,
 
 def major_result():
     uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.my_bayes,
-                           ConfidenceIntervalMethodFromExtremes.ci_mle][:]
+                           ConfidenceIntervalMethodFromExtremes.ci_mle][1:]
     massif_names = None
-    for study_class in paper_study_classes[:1]:
+    for study_class in paper_study_classes[:2]:
         if study_class == CrocusSnowLoadEurocode:
             non_stationary_uncertainty = [False]
         else:
@@ -89,21 +91,21 @@ def major_result():
 
 
 if __name__ == '__main__':
-    # major_result()
-    # intermediate_result(altitudes=paper_altitudes, massif_names=['Maurienne'],
+    major_result()
+    # intermediate_result(altitudes=[900, 1200], massif_names=['Belledonne'],
     #                     uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
-    #                        ConfidenceIntervalMethodFromExtremes.ci_mle][1:],
+    #                                          ConfidenceIntervalMethodFromExtremes.ci_mle][1:],
     #                     non_stationary_uncertainty=[False, True][1:],
     #                     multiprocessing=True)
-    # intermediate_result(altitudes=paper_altitudes, massif_names=['Maurienne'],
+    # intermediate_result(altitudes=[900, 1200], massif_names=['Maurienne'],
     #                     uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.my_bayes,
-    #                                          ConfidenceIntervalMethodFromExtremes.ci_mle][:],
+    #                                          ConfidenceIntervalMethodFromExtremes.ci_mle][1:],
     #                     non_stationary_uncertainty=[False, True][:],
     #                     multiprocessing=True)
     # intermediate_result(altitudes=[900, 1200], massif_names=None)
     # intermediate_result(ALL_ALTITUDES_WITHOUT_NAN)
     # intermediate_result(paper_altitudes)
-    minor_result(altitude=900)
+    # minor_result(altitude=900)
     # intermediate_result(altitudes=[1500, 1800], massif_names=['Chartreuse'],
     #                     uncertainty_methods=[ConfidenceIntervalMethodFromExtremes.ci_mle,
     #                                          ConfidenceIntervalMethodFromExtremes.ci_bayes],
diff --git a/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py b/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py
index 8c0ac530463856bde162a0e1c42aa24845ea9715..41406af0c0f250680125915ccdb111134a6bfe09 100644
--- a/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py
+++ b/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_curves.py
@@ -69,8 +69,8 @@ def plot_single_uncertainty_massif(altitude_to_visualizer: Dict[int, StudyVisual
 def get_label_name(non_stationary_context, ci_method_name):
     model_symbol = 'N' if non_stationary_context else '0'
     parameter = ', 2017' if non_stationary_context else ''
-    model_name = ' $ \widehat{z_p}(\\boldsymbol{\\theta_{\mathcal{M}_'
-    model_name += model_symbol
+    model_name = ' $ \widehat{z_p}(\\boldsymbol{\\theta_{\mathcal{M}'
+    # model_name += '_' + model_symbol
     model_name += '}}'
     model_name += parameter
     model_name += ')_{ \\textrm{' + ci_method_name.upper().split(' ')[1] + '}} $ '
@@ -100,6 +100,7 @@ def plot_single_uncertainty_massif_and_non_stationary_context(ax, massif_name, n
         color = ci_method_to_color[uncertainty_method]
         valid_altitudes = plot_valid_return_level_uncertainties(alpha, altitude_to_visualizer, altitudes, ax, color,
                                                                 massif_name, non_stationary_context, uncertainty_method)
+        # Plot some data for the non valid altitudes
 
         # Plot bars of TDRL only in the non stationary case
         if j == 0 and non_stationary_context:
@@ -133,8 +134,8 @@ def add_title(ax, eurocode_region, massif_name, non_stationary_context):
 
 
 def plot_tdrl_bars(altitude_to_visualizer, ax, massif_name, valid_altitudes, legend_size, fontsize):
-    visualizers = [v for a, v in altitude_to_visualizer.items() if
-                   a in valid_altitudes and massif_name in v.uncertainty_massif_names]
+    visualizers = [v for a, v in altitude_to_visualizer.items()
+                   if a in valid_altitudes and massif_name in v.uncertainty_massif_names]
     if len(visualizers) > 0:
         tdrl_values = [v.massif_name_to_tdrl_value[massif_name] for v in visualizers]
         # Plot bars
@@ -151,15 +152,17 @@ def plot_tdrl_bars(altitude_to_visualizer, ax, massif_name, valid_altitudes, leg
             # ax.plot([altitude], [value / 2], **marker_kwargs)
             # Better to plot all the markers on the same line
             ax.plot([altitude], 0, **marker_kwargs)
-    # Add a legend plot
-    visualizer = visualizers[0]
-    legend_elements = AbstractStudy.get_legend_for_model_symbol(visualizer.marker_to_label, markersize=9)
-    ax2 = ax.twinx()
-    # ax2.legend(handles=legend_elements, bbox_to_anchor=(0.93, 0.7), loc='upper right')
-    # ax2.annotate("Filled symbol = significant trend ", xy=(0.85, 0.5), xycoords='axes fraction', fontsize=7)
-    ax2.legend(handles=legend_elements, loc='upper right', prop={'size': legend_size})
-    ax2.annotate("Filled symbol =\n significant trend ", xy=(0.6, 0.85), xycoords='axes fraction', fontsize=fontsize)
-    ax2.set_yticks([])
+        # Add a legend plot
+        visualizer = visualizers[0]
+        markers = [v.massif_name_to_marker_style[massif_name]['marker'] for v in visualizers]
+        marker_to_label = {m: visualizer.all_marker_style_to_label_name[m] for m in markers}
+        legend_elements = AbstractStudy.get_legend_for_model_symbol(marker_to_label, markersize=9)
+        ax2 = ax.twinx()
+        # ax2.legend(handles=legend_elements, bbox_to_anchor=(0.93, 0.7), loc='upper right')
+        # ax2.annotate("Filled symbol = significant trend ", xy=(0.85, 0.5), xycoords='axes fraction', fontsize=7)
+        ax2.legend(handles=legend_elements, loc='upper right', prop={'size': legend_size})
+        ax2.annotate("Filled symbol =\nsignificant trend  \nw.r.t $\mathcal{M}_0$", xy=(0.6, 0.85), xycoords='axes fraction', fontsize=fontsize)
+        ax2.set_yticks([])
 
 
 def plot_valid_return_level_uncertainties(alpha, altitude_to_visualizer, altitudes, ax, color, massif_name,
diff --git a/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_histogram.py b/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_histogram.py
index 9fb61fd0b53d8d7fb5ca4614ac783d7cf1ab11b1..68e535ebb8d4856e78c65f783b38874dd8b7e84a 100644
--- a/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_histogram.py
+++ b/experiment/paper_past_snow_loads/result_trends_and_return_levels/plot_uncertainty_histogram.py
@@ -7,7 +7,7 @@ from experiment.paper_past_snow_loads.paper_utils import dpi_paper1_figure
 from experiment.paper_past_snow_loads.study_visualizer_for_non_stationary_trends import \
     StudyVisualizerForNonStationaryTrends
 from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import ci_method_to_color, \
-    ci_method_to_label
+    ci_method_to_label, ConfidenceIntervalMethodFromExtremes
 
 
 def plot_uncertainty_histogram(altitude_to_visualizer: Dict[int, StudyVisualizerForNonStationaryTrends]):
@@ -36,13 +36,17 @@ def plot_histogram(altitude_to_visualizer, non_stationary_context):
         if len(visualizer.uncertainty_methods) == 2:
             offset = -50 if j == 0 else 50
             bincenters = altitudes + offset
-        width = 100
+            width = 100
+        else:
+            width = 200
         plot_histogram_ci_method(visualizers, non_stationary_context, ci_method, ax, bincenters, width=width)
     fontsize_label = 15
     legend_size = 15
     ax.set_xticks(altitudes)
     ax.tick_params(labelsize=fontsize_label)
-    ax.legend(loc='upper left', prop={'size': legend_size})
+    if not (len(visualizer.uncertainty_methods) == 1
+            and visualizer.uncertainty_methods[0] == ConfidenceIntervalMethodFromExtremes.ci_mle):
+        ax.legend(loc='upper left', prop={'size': legend_size})
     ax.set_ylabel('Massifs exceeding French standards (\%)', fontsize=fontsize_label)
     ax.set_xlabel('Altitude (m)', fontsize=fontsize_label)
     ax.set_ylim([0, 100])
diff --git a/experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py b/experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py
index acb422a39f9814ffcbfb0f433cce062bc8b899c4..862ac213818193a02ce247c3548a09ea6cc808c5 100644
--- a/experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py
+++ b/experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py
@@ -27,7 +27,8 @@ from experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_two
     GevLocationAndScaleTrendTest, GevLocationAgainstGumbel, GevScaleAgainstGumbel
 from experiment.trend_analysis.univariate_test.extreme_trend_test.trend_test_two_parameters.gumbel_test_two_parameters import \
     GumbelLocationAndScaleTrendTest
-from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import StationaryTemporalModel
+from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import StationaryTemporalModel, \
+    GumbelTemporalModel
 from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import \
     ConfidenceIntervalMethodFromExtremes
 from extreme_fit.model.result_from_model_fit.result_from_extremes.eurocode_return_level_uncertainties import \
@@ -78,7 +79,6 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
                                                                 # ["v", "^", "D", "X", "x", 7, 6, "d"]))
         else:
             self.non_stationary_trend_test = list(self.non_stationary_trend_test_to_marker.keys())
-        self.marker_to_label = OrderedDict([(t.marker, t.label) for t in self.non_stationary_trend_test])
         self.global_max_abs_change = None
 
     # Utils
@@ -137,14 +137,12 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
     def plot_trends(self, max_abs_tdrl=None,  add_colorbar=True):
         if max_abs_tdrl is not None:
             self.global_max_abs_change = max_abs_tdrl
-        marker_style_selected = set([d['marker'] for d in self.massif_name_to_marker_style.values()])
-        marker_style_to_label_name = {m: l for m, l in self.marker_to_label.items() if m in marker_style_selected}
         ax = self.study.visualize_study(massif_name_to_value=self.massif_name_to_change_value,
                                         replace_blue_by_white=False,
                                         axis_off=False, show_label=False,
                                         add_colorbar=add_colorbar,
                                         massif_name_to_marker_style=self.massif_name_to_marker_style,
-                                        marker_style_to_label_name=marker_style_to_label_name,
+                                        marker_style_to_label_name=self.selected_marker_style_to_label_name,
                                         massif_name_to_color=self.massif_name_to_color,
                                         cmap=self.cmap,
                                         show=False,
@@ -158,6 +156,17 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
                                   dpi=500)
         plt.close()
 
+    @cached_property
+    def all_marker_style_to_label_name(self):
+        return OrderedDict([(t.marker, t.label) for t in self.non_stationary_trend_test])
+
+    @cached_property
+    def selected_marker_style_to_label_name(self):
+        marker_style_selected = set([d['marker'] for d in self.massif_name_to_marker_style.values()])
+        marker_style_to_label_name = {m: l for m, l in self.all_marker_style_to_label_name.items()
+                                      if m in marker_style_selected}
+        return marker_style_to_label_name
+
     @property
     def label(self):
         if self.relative_change_trend_plot:
@@ -238,7 +247,7 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
 
     def massif_name_and_non_stationary_context_to_model_class(self, massif_name, non_stationary_context):
         if not non_stationary_context:
-            return StationaryTemporalModel
+            return GumbelTemporalModel
         else:
             return self.massif_name_to_minimized_aic_non_stationary_trend_test[massif_name].unconstrained_model_class
 
diff --git a/experiment/trend_analysis/univariate_test/extreme_trend_test/trend_test_one_parameter/gumbel_trend_test_one_parameter.py b/experiment/trend_analysis/univariate_test/extreme_trend_test/trend_test_one_parameter/gumbel_trend_test_one_parameter.py
index 698a7063dc4d49b705e658a6e7a8e34dc0a99913..927b7371fd0822e1c002d66a69bc8b0aabe5e950 100644
--- a/experiment/trend_analysis/univariate_test/extreme_trend_test/trend_test_one_parameter/gumbel_trend_test_one_parameter.py
+++ b/experiment/trend_analysis/univariate_test/extreme_trend_test/trend_test_one_parameter/gumbel_trend_test_one_parameter.py
@@ -24,10 +24,6 @@ class GumbelVersusGumbel(GevTrendTestOneParameter):
     def total_number_of_parameters_for_unconstrained_model(self) -> int:
         return 2
 
-    @property
-    def time_derivative_of_return_level(self):
-        return super().time_derivative_of_return_level()
-
     @classproperty
     def label(self):
         return super().label % '0'
diff --git a/extreme_fit/distribution/gev/ci_fevd_fixed.R b/extreme_fit/distribution/gev/ci_fevd_fixed.R
index 3b5a7ba4e5bd1e313303f78fc7e1ca792a14090b..785b3ca3d8a7c347581a99aa853ca4eb2f1dd959 100644
--- a/extreme_fit/distribution/gev/ci_fevd_fixed.R
+++ b/extreme_fit/distribution/gev/ci_fevd_fixed.R
@@ -1,11 +1,15 @@
 
 # todo: report bug on exTremes, fix is around line 510 where I set theta[4] = 0 (and it was Nan before)
 # it was important to do that to make the code work in the Gumbel case
+# todo: also luine 11  x$results$num.pars$shape = 0
+# Also aroudn line 433,  grads = grads[1:dim(cov.theta)[1]]
 ci.fevd.mle_fixed <- function (x, alpha = 0.05, type = c("return.level", "parameter"),
     return.period = 100, which.par = 1, R = 502, method = c("normal",
         "boot", "proflik"), xrange = NULL, nint = 20, verbose = FALSE,
     tscale = FALSE, return.samples = FALSE, ...)
 {
+    if (x$type == "Gumbel")
+        x$results$num.pars$shape = 0
     if (missing(method))
         miss.meth <- TRUE
     else miss.meth <- FALSE
@@ -13,6 +17,7 @@ ci.fevd.mle_fixed <- function (x, alpha = 0.05, type = c("return.level", "parame
     method <- match.arg(method)
     type <- tolower(type)
     type <- match.arg(type)
+
     theta.hat <- x$results$par
     theta.names <- names(theta.hat)
     np <- length(theta.hat)
@@ -426,6 +431,7 @@ ci.rl.ns.fevd.mle_fixed <- function (x, alpha = 0.05, return.period = 100, metho
             stop("ci: negative Std. Err. estimates obtained.  Not trusting any of them.")
         grads <- t(rlgrad.fevd(x, period = return.period, qcov = qcov,
             qcov.base = qcov.base))
+        grads = grads[1:dim(cov.theta)[1]]
         se.theta <- sqrt(diag(t(grads) %*% cov.theta %*% grads))
         out <- cbind(c(res) - z.alpha * se.theta, c(res), c(res) +
             z.alpha * se.theta, se.theta)
diff --git a/extreme_fit/distribution/gev/main_fevd_mle.R b/extreme_fit/distribution/gev/main_fevd_mle.R
index cf67516c5419d3fb3a21261553dfb8f6ebf416a2..592dbaacb6b73d53308b5046c37c67a099c51a74 100644
--- a/extreme_fit/distribution/gev/main_fevd_mle.R
+++ b/extreme_fit/distribution/gev/main_fevd_mle.R
@@ -11,14 +11,19 @@ source('ci_fevd_fixed.R')
 # Sample from a GEV
 set.seed(42)
 N <- 50
-loc = 0; scale = 1; shape <- 0.1
+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
+# N <- 50
+# loc = 0; scale = 1; shape <- 0.1
+# x_gev <- rgev(N, loc = loc, scale = scale, shape = shape)
 coord <- matrix(ncol=1, nrow = N)
 coord[,1]=seq(0,N-1,1)
 colnames(coord) = c("T")
 coord = data.frame(coord, stringsAsFactors = TRUE)
 # res = fevd_fixed(x_gev, data=coord, method='MLE', verbose=TRUE, use.phi=FALSE)
-res = fevd_fixed(x_gev, data=coord, location.fun= ~T, method='MLE', type="Gumbel", verbose=FALSE, use.phi=FALSE)
+# res = fevd_fixed(x_gev, data=coord, location.fun= ~T, scale.fun= ~T, method='MLE', type="GEV", verbose=FALSE, use.phi=FALSE)
+res = fevd_fixed(x_gev, data=coord, location.fun= ~T, method='MLE', type="GEV", verbose=FALSE, use.phi=FALSE)
 print(res)
 
 # Some display for the results
@@ -40,11 +45,11 @@ print(res$results$par)
 
 # Bug to solve for the non stationary - the returned parameter do not match with the return level
 v = make.qcov(res, vals = list(mu1 = c(0.0)))
-print(res$results$num.pars$location)
-print(res$results$num.pars$scale)
-print(res$results$num.pars$shape)
-res$results$num.pars$shape = 0
-print(res$results$num.pars$shape)
+# print(res$results$num.pars$location)
+# print(res$results$num.pars$scale)
+# print(res$results$num.pars$shape)
+# res$results$num.pars$shape = 0
+# print(res$results$num.pars$shape)
 
 # r = return.level(res, return.period = 50, qcov = v, type='Gumbel')
 # print(r)
diff --git a/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_mle_extremes.py b/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_mle_extremes.py
index 9fa9b18ff8842b27f80bd0cc434eba4c3c08b50c..3cf68a4673ec6eb6081b125dd02223bc1626a81e 100644
--- a/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_mle_extremes.py
+++ b/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_mle_extremes.py
@@ -30,7 +30,7 @@ class ResultFromMleExtremes(AbstractResultFromExtremes):
                 'method': method_name,
             # xrange = NULL, nint = 20
         }
-        res = r.ci(self.result_from_fit, **mle_ci_parameters, **common_kwargs)
+        res = r('ci.fevd.mle_fixed')(self.result_from_fit, **mle_ci_parameters, **common_kwargs)
         if self.is_non_stationary:
             a = np.array(res)[0]
             lower, mean_estimate, upper, _ = a
diff --git a/extreme_fit/model/utils.py b/extreme_fit/model/utils.py
index 90dfd27877ed7d2702b4bcbd45134e61c6a554cd..658c2fcbade7865c3538945c60baac81638f665d 100644
--- a/extreme_fit/model/utils.py
+++ b/extreme_fit/model/utils.py
@@ -29,9 +29,10 @@ warnings.filterwarnings("ignore")
 # Load ismev
 r.library('ismev')
 # Load fevd fixed
-fevd_fixed_filepath = op.join(get_root_path(), 'extreme_fit', 'distribution', 'gev', 'fevd_fixed.R')
-assert op.exists(fevd_fixed_filepath)
-r.source(fevd_fixed_filepath)
+for filename in ['ci_fevd_fixed.R', 'fevd_fixed.R']:
+    fevd_fixed_filepath = op.join(get_root_path(), 'extreme_fit', 'distribution', 'gev', filename)
+    assert op.exists(fevd_fixed_filepath)
+    r.source(fevd_fixed_filepath)
 # Reactivate warning
 warnings.filters = default_filters
 
diff --git a/test/test_extreme_fit/test_model/test_confidence_interval.py b/test/test_extreme_fit/test_model/test_confidence_interval.py
index 4b3e7038353b9e129bf516c49e2e4cc0f2282a91..1d8795c568a91f99262c2aace1b34841b5a7d783 100644
--- a/test/test_extreme_fit/test_model/test_confidence_interval.py
+++ b/test/test_extreme_fit/test_model/test_confidence_interval.py
@@ -76,10 +76,10 @@ class TestConfidenceInterval(unittest.TestCase):
     def test_ci_normal(self):
         self.fit_method = TemporalMarginFitMethod.extremes_fevd_mle
         self.ci_method = ConfidenceIntervalMethodFromExtremes.ci_mle
-        self.model_class_to_triplet= {
+        self.model_class_to_triplet = {
             StationaryTemporalModel: (-4.703945484843988, 30.482318639674023, 65.66858276419204),
-            NonStationaryLocationTemporalModel: (-4.223086740397132, 30.29842988666537, 64.81994651372787),
-            NonStationaryLocationAndScaleTemporalModel: (-15.17041284612494, 43.69511224410276, 102.56063733433047),
+            NonStationaryLocationTemporalModel: (-30.361576509947707, 4.159940117114796, 38.6814567441773),
+            NonStationaryLocationAndScaleTemporalModel: (-52.797816369170455, 6.0677087210572465, 64.93323381128495),
         }
 
     def test_ci_boot(self):