[PAPER 1] add gumbel type for mle fit. add...

[PAPER 1] add gumbel type for mle fit. add study_visualizer_for_shape_repartition both for gumbel VS stationary and stationary VS non-stationary

experiment/meteo_france_data/scm_models_data/visualization/study_visualization/study_visualizer.py

@@ -632,7 +632,10 @@ class StudyVisualizer(VisualizationParameters):
         # Display the graph of the max on top
         if ax is None:
             ax = plt.gca()
-        x, y = self.smooth_maxima_x_y(massif_names.index(massif_name))
+        if massif_name is None:
+            x, y = self.study.ordered_years, self.study.observations_annual_maxima.df_maxima_gev.mean(axis=0)
+        else:
+            x, y = self.smooth_maxima_x_y(massif_names.index(massif_name))
         ax.plot(x, y, color=color, linewidth=5, label=label, linestyle=linestyle)
         # ax.set_ylabel('{} (in {})'.format(snow_abbreviation, self.study.variable_unit), color=color, fontsize=15)
 
@@ -642,7 +645,7 @@ class StudyVisualizer(VisualizationParameters):
         self.plot_name = plot_name
         ax.set_ylabel('{} ({})'.format(snow_abbreviation, self.study.variable_unit), fontsize=15)
         ax.set_xlabel('years', fontsize=15)
-        if label is None:
+        if label is None and massif_name is not None:
             ax.set_title('{} at {} m'.format(massif_name, altitude))
         if last_plot:
             ax.legend()

experiment/paper_past_snow_loads/check_mle_convergence_for_trends/main_shape_repartition.py

+from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal
+from experiment.paper_past_snow_loads.check_mle_convergence_for_trends.study_visualizer_for_shape_repartition import \
+    StudyVisualizerForShape, StudyVisualizerGumbel
+from experiment.paper_past_snow_loads.paper_main_utils import load_altitude_to_visualizer
+
+
+def main_shape_repartition(altitudes, massif_names=None,
+                           non_stationary_uncertainty=None, uncertainty_methods=None,
+                           study_class=CrocusSnowLoadTotal,
+                           study_visualizer_class=StudyVisualizerForShape, save_to_file=True):
+    # Load altitude to visualizer
+    altitude_to_visualizer = load_altitude_to_visualizer(altitudes, massif_names, non_stationary_uncertainty,
+                                                         study_class, uncertainty_methods,
+                                                         study_visualizer_class=study_visualizer_class,
+                                                         save_to_file=save_to_file)
+    altitudes_for_plot_trend = altitudes  #[900, 1800, 2700]
+    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=True)
+
+
+if __name__ == '__main__':
+    # main_shape_repartition([900], save_to_file=False)
+    # main_shape_repartition([900, 1800, 2700])
+    # main_shape_repartition([300, 600, 900, 1200, 1500, 1800, 2700])
+    main_shape_repartition([900], study_visualizer_class=StudyVisualizerGumbel)

experiment/paper_past_snow_loads/check_mle_convergence_for_trends/study_visualizer_for_shape_repartition.py

+from cached_property import cached_property
+
+from experiment.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
+from experiment.paper_past_snow_loads.study_visualizer_for_non_stationary_trends import \
+    StudyVisualizerForNonStationaryTrends
+from experiment.trend_analysis.abstract_score import MeanScore
+from experiment.trend_analysis.univariate_test.gev_trend_test_one_parameter import GevLocationTrendTest, \
+    GevStationaryVersusGumbel
+
+
+class StudyVisualizerForShape(StudyVisualizerForNonStationaryTrends):
+
+    def __init__(self, study: AbstractStudy, show=True, save_to_file=False, only_one_graph=False, only_first_row=False,
+                 vertical_kde_plot=False, year_for_kde_plot=None, plot_block_maxima_quantiles=False,
+                 temporal_non_stationarity=False, transformation_class=None, verbose=False, multiprocessing=False,
+                 complete_non_stationary_trend_analysis=False, normalization_under_one_observations=True,
+                 score_class=MeanScore, uncertainty_methods=None, non_stationary_contexts=None,
+                 uncertainty_massif_names=None, effective_temporal_covariate=2017, relative_change_trend_plot=True):
+        super().__init__(study, show, save_to_file, only_one_graph, only_first_row, vertical_kde_plot,
+                         year_for_kde_plot, plot_block_maxima_quantiles, temporal_non_stationarity,
+                         transformation_class, verbose, multiprocessing, complete_non_stationary_trend_analysis,
+                         normalization_under_one_observations, score_class, uncertainty_methods,
+                         non_stationary_contexts, uncertainty_massif_names, effective_temporal_covariate,
+                         relative_change_trend_plot)
+
+    @cached_property
+    def massif_name_to_unconstrained_shape_parameter(self):
+        return {m: t.unconstrained_estimator_gev_params.shape
+                for m, t in self.massif_name_to_minimized_aic_non_stationary_trend_test.items()}
+
+    @cached_property
+    def massif_name_to_change_value(self):
+        print(self.massif_name_to_unconstrained_shape_parameter)
+        return self.massif_name_to_unconstrained_shape_parameter
+
+    @property
+    def label(self):
+        return 'Shape parameter value'
+
+    @property
+    def graduation(self):
+        return 0.1
+
+
+class StudyVisualizerGumbel(StudyVisualizerForNonStationaryTrends):
+
+    def __init__(self, study: AbstractStudy, show=True, save_to_file=False, only_one_graph=False, only_first_row=False,
+                 vertical_kde_plot=False, year_for_kde_plot=None, plot_block_maxima_quantiles=False,
+                 temporal_non_stationarity=False, transformation_class=None, verbose=False, multiprocessing=False,
+                 complete_non_stationary_trend_analysis=False, normalization_under_one_observations=True,
+                 score_class=MeanScore, uncertainty_methods=None, non_stationary_contexts=None,
+                 uncertainty_massif_names=None, effective_temporal_covariate=2017, relative_change_trend_plot=True):
+        super().__init__(study, show, save_to_file, only_one_graph, only_first_row, vertical_kde_plot,
+                         year_for_kde_plot, plot_block_maxima_quantiles, temporal_non_stationarity,
+                         transformation_class, verbose, multiprocessing, complete_non_stationary_trend_analysis,
+                         normalization_under_one_observations, score_class, uncertainty_methods,
+                         non_stationary_contexts, uncertainty_massif_names, effective_temporal_covariate,
+                         relative_change_trend_plot)
+        # Assign default argument for the non stationary trends
+        self.non_stationary_trend_test = [GevStationaryVersusGumbel]
+        self.non_stationary_trend_test_to_marker = dict(zip(self.non_stationary_trend_test, ["o"]))

experiment/paper_past_snow_loads/discussion_data_comparison_with_eurocode/main_comparison_with_eurocode_global.py

+from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusDepth
+from experiment.meteo_france_data.scm_models_data.crocus.crocus_variables import CrocusDepthVariable
+from experiment.meteo_france_data.scm_models_data.visualization.study_visualization.main_study_visualizer import \
+    study_iterator_global, SCM_STUDY_CLASS_TO_ABBREVIATION, snow_density_str
+from experiment.meteo_france_data.scm_models_data.visualization.study_visualization.study_visualizer import \
+    StudyVisualizer
+import matplotlib.pyplot as plt
+
+from experiment.paper_past_snow_loads.discussion_data_comparison_with_eurocode.crocus_study_comparison_with_eurocode import \
+    CrocusDifferenceSnowLoad, \
+    CrocusSnowDensityAtMaxofSwe, CrocusDifferenceSnowLoadRescaledAndEurocodeToSeeSynchronization, \
+    CrocusSnowDepthAtMaxofSwe, CrocusSnowDepthDifference
+from experiment.paper_past_snow_loads.paper_utils import dpi_paper1_figure
+
+
+def max_graph_annual_maxima_comparison():
+    """
+    We choose these massif because each represents a different eurocode region
+    we also choose them because they belong to a different climatic area
+    :return:
+    """
+    save_to_file = False
+    study_classes = [
+                        CrocusSnowDensityAtMaxofSwe,
+                        CrocusDifferenceSnowLoad,
+                        CrocusSnowDepthDifference,
+                    ][:]
+    study_class_to_ylim_and_yticks = {
+        CrocusSnowDensityAtMaxofSwe: ([150, 500], [50 * i for i in range(3, 11)]),
+        CrocusDifferenceSnowLoad: ([0, 10], [2 * i for i in range(0, 6)]),
+        CrocusSnowDepthDifference: ([0, 0.6], [0.2 * i for i in range(0, 4)]),
+    }
+    for study_class in study_classes:
+        ylim, yticks = study_class_to_ylim_and_yticks[study_class]
+
+        marker_altitude = [
+                              ('deepskyblue', 900),
+                              ('dodgerblue', 1800),
+                              ('blue', 2700),
+                          ][:]
+        ax = plt.gca()
+        for color, altitude in marker_altitude:
+            for study in study_iterator_global([study_class], altitudes=[altitude]):
+                study_visualizer = StudyVisualizer(study, save_to_file=save_to_file,
+                                                   verbose=True,
+                                                   multiprocessing=True)
+                snow_abbreviation = SCM_STUDY_CLASS_TO_ABBREVIATION[study_class]
+                nb_massifs = len(study_visualizer.study.study_massif_names)
+                label = 'Mean value at {}m (out of {} massifs)'.format(altitude, nb_massifs)
+                study_visualizer.visualize_max_graphs_poster(None, altitude, snow_abbreviation, color, label,
+                                                             False, ax)
+
+                last_plot = altitude == 2700
+                if last_plot:
+                    # if study_class == CrocusSnowDensityAtMaxofSwe:
+                    #     label = '{} for French standards'.format(snow_density_str)
+                    #     snow_density_eurocode = [150 for _ in study.ordered_years]
+                    #     ax.plot(study.ordered_years, snow_density_eurocode, color='k', label=label)
+                    ax.legend()
+                    tight_pad = {'h_pad': 0.2}
+                    ax.set_ylim(ylim)
+                    ax.set_xlim([1957, 2018])
+                    ax.yaxis.set_ticks(yticks)
+                    study_visualizer.show_or_save_to_file(no_title=True, tight_layout=True,
+                                                          tight_pad=tight_pad, dpi=dpi_paper1_figure)
+                    ax.clear()
+
+
+if __name__ == '__main__':
+    max_graph_annual_maxima_comparison()

experiment/paper_past_snow_loads/paper_main_utils.py

@@ -4,11 +4,13 @@ from experiment.paper_past_snow_loads.study_visualizer_for_non_stationary_trends
     StudyVisualizerForNonStationaryTrends
 
 
-def load_altitude_to_visualizer(altitudes, massif_names, non_stationary_uncertainty, study_class, uncertainty_methods):
+def load_altitude_to_visualizer(altitudes, massif_names, non_stationary_uncertainty, study_class, uncertainty_methods,
+                                study_visualizer_class=StudyVisualizerForNonStationaryTrends,
+                                save_to_file=True):
     altitude_to_visualizer = OrderedDict()
     for altitude in altitudes:
-        altitude_to_visualizer[altitude] = StudyVisualizerForNonStationaryTrends(
-            study=study_class(altitude=altitude), multiprocessing=True, save_to_file=True,
+        altitude_to_visualizer[altitude] = study_visualizer_class(
+            study=study_class(altitude=altitude), multiprocessing=True, save_to_file=save_to_file,
             uncertainty_massif_names=massif_names, uncertainty_methods=uncertainty_methods,
             non_stationary_contexts=non_stationary_uncertainty)
     return altitude_to_visualizer

experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py

@@ -133,7 +133,6 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
         ax.get_xaxis().set_visible(True)
         ax.set_xticks([])
         ax.set_xlabel('Altitude = {}m'.format(self.study.altitude), fontsize=15)
-
         self.plot_name = 'tdlr_trends'
         self.show_or_save_to_file(add_classic_title=False, tight_layout=True, no_title=True,
                                   dpi=500)
@@ -157,16 +156,20 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
 
     @property
     def ticks_values_and_labels(self):
-        graduation = 10 if self.relative_change_trend_plot else 0.2
         positive_ticks = []
-        tick = graduation
+        tick = self.graduation
         while tick < self._max_abs_change:
             positive_ticks.append(round(tick, 1))
-            tick += graduation
+            tick += self.graduation
         all_ticks_labels = [-t for t in positive_ticks] + [0] + positive_ticks
         ticks_values = [((t / self._max_abs_change) + 1) / 2 for t in all_ticks_labels]
         return ticks_values, all_ticks_labels
 
+    @property
+    def graduation(self):
+        graduation = 10 if self.relative_change_trend_plot else 0.2
+        return graduation
+
     @cached_property
     def massif_name_to_tdrl_value(self):
         return {m: t.time_derivative_of_return_level for m, t in
@@ -291,3 +294,8 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
                        uncertainty.confidence_interval[1] > eurocode)
                       for eurocode, uncertainty in eurocode_and_uncertainties])
         return 100 * np.mean(a, axis=0)
+
+
+
+
+

experiment/trend_analysis/univariate_test/abstract_gev_trend_test.py

@@ -132,13 +132,13 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
                                                                               AbstractCoordinates.COORDINATE_T)
 
     @cached_property
-    def non_stationary_constant_gev_params(self) -> GevParams:
+    def unconstrained_estimator_gev_params(self) -> GevParams:
         # Constant parameters correspond to the gev params in 1958
         return self.unconstrained_estimator.margin_function_from_fit.get_gev_params(coordinate=np.array([1958]),
                                                                                     is_transformed=False)
 
     @cached_property
-    def stationary_constant_gev_params(self) -> GevParams:
+    def constrained_estimator_gev_params(self) -> GevParams:
         # Constant parameters correspond to any gev params
         return self.constrained_estimator.margin_function_from_fit.get_gev_params(coordinate=np.array([1958]),
                                                                                   is_transformed=False)
@@ -192,4 +192,4 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
         if self.crashed:
             return 0.0
         else:
-            return self.non_stationary_constant_gev_params.quantile(p=self.quantile_level)
+            return self.unconstrained_estimator_gev_params.quantile(p=self.quantile_level)

experiment/trend_analysis/univariate_test/gev_trend_test_one_parameter.py

 from experiment.eurocode_data.utils import EUROCODE_QUANTILE
 from experiment.trend_analysis.univariate_test.abstract_gev_trend_test import AbstractGevTrendTest
 from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import \
-    NonStationaryLocationTemporalModel, NonStationaryScaleTemporalModel, NonStationaryShapeTemporalModel
+    NonStationaryLocationTemporalModel, NonStationaryScaleTemporalModel, NonStationaryShapeTemporalModel, \
+    StationaryTemporalModel, GumbelTemporalModel
 from extreme_fit.distribution.gev.gev_params import GevParams
 
 
 class GevTrendTestOneParameter(AbstractGevTrendTest):
 
-    def __init__(self, years, maxima, starting_year, unconstrained_model_class, gev_param_name, quantile_level=EUROCODE_QUANTILE):
-        super().__init__(years, maxima, starting_year,
-                         unconstrained_model_class=unconstrained_model_class,
-                         quantile_level=quantile_level)
-        self.gev_param_name = gev_param_name
-
     @property
     def degree_freedom_chi2(self) -> int:
         return 1
 
+
+class GevTrendTestOneParameterAgainstStationary(GevTrendTestOneParameter):
+
+    def __init__(self, years, maxima, starting_year, unconstrained_model_class, gev_param_name,
+                 quantile_level=EUROCODE_QUANTILE,
+                 constrained_model_class=StationaryTemporalModel):
+        super().__init__(years, maxima, starting_year,
+                         unconstrained_model_class=unconstrained_model_class,
+                         quantile_level=quantile_level,
+                         constrained_model_class=constrained_model_class)
+        self.gev_param_name = gev_param_name
+
     @property
     def non_stationary_linear_coef(self):
         return self.get_non_stationary_linear_coef(gev_param_name=self.gev_param_name)
 
 
-class GevLocationTrendTest(GevTrendTestOneParameter):
+class GevLocationTrendTest(GevTrendTestOneParameterAgainstStationary):
 
     def __init__(self, years, maxima, starting_year, quantile_level=EUROCODE_QUANTILE):
         super().__init__(years, maxima, starting_year,
@@ -31,12 +38,12 @@ class GevLocationTrendTest(GevTrendTestOneParameter):
                          quantile_level=quantile_level)
 
     def _slope_strength(self):
-        return self.non_stationary_constant_gev_params.time_derivative_of_return_level(p=self.quantile_level,
+        return self.unconstrained_estimator_gev_params.time_derivative_of_return_level(p=self.quantile_level,
                                                                                        mu1=self.non_stationary_linear_coef)
 
     @property
     def mean_difference_same_sign_as_slope_strenght(self) -> bool:
-        zeta0 = self.non_stationary_constant_gev_params.shape
+        zeta0 = self.unconstrained_estimator_gev_params.shape
         mean_difference = self.mean_difference(zeta0=zeta0, mu1=self.non_stationary_linear_coef)
         return self.same_sign(mean_difference, self._slope_strength())
 
@@ -45,7 +52,7 @@ class GevLocationTrendTest(GevTrendTestOneParameter):
         return False
 
 
-class GevScaleTrendTest(GevTrendTestOneParameter):
+class GevScaleTrendTest(GevTrendTestOneParameterAgainstStationary):
 
     def __init__(self, years, maxima, starting_year, quantile_level=EUROCODE_QUANTILE):
         super().__init__(years, maxima, starting_year,
@@ -54,13 +61,13 @@ class GevScaleTrendTest(GevTrendTestOneParameter):
                          quantile_level=quantile_level)
 
     def _slope_strength(self):
-        return self.non_stationary_constant_gev_params.time_derivative_of_return_level(
+        return self.unconstrained_estimator_gev_params.time_derivative_of_return_level(
             p=self.quantile_level,
             sigma1=self.non_stationary_linear_coef)
 
     @property
     def mean_difference_same_sign_as_slope_strenght(self) -> bool:
-        zeta0 = self.non_stationary_constant_gev_params.shape
+        zeta0 = self.unconstrained_estimator_gev_params.shape
         mean_difference = self.mean_difference(zeta0=zeta0, sigma1=self.non_stationary_linear_coef)
         return self.same_sign(mean_difference, self._slope_strength())
 
@@ -70,10 +77,19 @@ class GevScaleTrendTest(GevTrendTestOneParameter):
         return self.same_sign(sigma1, self._slope_strength())
 
 
-class GevShapeTrendTest(GevTrendTestOneParameter):
+class GevShapeTrendTest(GevTrendTestOneParameterAgainstStationary):
 
     def __init__(self, years, maxima, starting_year, quantile_level=EUROCODE_QUANTILE):
         super().__init__(years, maxima, starting_year,
                          unconstrained_model_class=NonStationaryShapeTemporalModel,
                          gev_param_name=GevParams.SHAPE,
                          quantile_level=quantile_level)
+
+
+class GevStationaryVersusGumbel(GevTrendTestOneParameter):
+
+    def __init__(self, years, maxima, starting_year, quantile_level=EUROCODE_QUANTILE):
+        super().__init__(years, maxima, starting_year,
+                         unconstrained_model_class=StationaryTemporalModel,
+                         constrained_model_class=GumbelTemporalModel,
+                         quantile_level=quantile_level)

experiment/trend_analysis/univariate_test/gev_trend_test_two_parameters.py

@@ -29,13 +29,13 @@ class GevLocationAndScaleTrendTest(GevTrendTestTwoParameters):
         return self.get_non_stationary_linear_coef(gev_param_name=GevParams.SCALE)
 
     def _slope_strength(self):
-        return self.non_stationary_constant_gev_params.time_derivative_of_return_level(p=self.quantile_level,
+        return self.unconstrained_estimator_gev_params.time_derivative_of_return_level(p=self.quantile_level,
                                                                                        mu1=self.mu1,
                                                                                        sigma1=self.sigma1)
 
     @property
     def mean_difference_same_sign_as_slope_strenght(self) -> bool:
-        zeta0 = self.non_stationary_constant_gev_params.shape
+        zeta0 = self.unconstrained_estimator_gev_params.shape
         mean_difference = self.mean_difference(zeta0=zeta0, mu1=self.mu1, sigma1=self.sigma1)
         return self.same_sign(mean_difference, self._slope_strength())
 

extreme_fit/model/margin_model/linear_margin_model/abstract_temporal_linear_margin_model.py

@@ -25,8 +25,10 @@ class AbstractTemporalLinearMarginModel(LinearMarginModel):
 
     def __init__(self, coordinates: AbstractCoordinates, use_start_value=False, params_start_fit=None,
                  params_sample=None, starting_point=None, fit_method=TemporalMarginFitMethod.is_mev_gev_fit,
-                 nb_iterations_for_bayesian_fit=5000, params_start_fit_bayesian=None):
+                 nb_iterations_for_bayesian_fit=5000, params_start_fit_bayesian=None,
+                 type_for_MLE="GEV"):
         super().__init__(coordinates, use_start_value, params_start_fit, params_sample, starting_point)
+        self.type_for_mle = type_for_MLE
         self.params_start_fit_bayesian = params_start_fit_bayesian
         self.nb_iterations_for_bayesian_fit = nb_iterations_for_bayesian_fit
         assert isinstance(fit_method, TemporalMarginFitMethod)
@@ -59,10 +61,12 @@ class AbstractTemporalLinearMarginModel(LinearMarginModel):
         res = safe_run_r_estimator(function=r('fevd_fixed'),
                                    x=x,
                                    data=y,
-                                   method='MLE',
+                                   type=self.type_for_mle,
+                                   method="MLE",
                                    **r_type_argument_kwargs
                                    )
-        return ResultFromMleExtremes(res, self.margin_function_start_fit.gev_param_name_to_dims)
+        return ResultFromMleExtremes(res, self.margin_function_start_fit.gev_param_name_to_dims,
+                                     type_for_mle=self.type_for_mle)
 
     def extremes_fevd_bayesian_fit(self, x, df_coordinates_temp) -> AbstractResultFromExtremes:
         r_type_argument_kwargs, y = self.extreme_arguments(df_coordinates_temp)

extreme_fit/model/margin_model/linear_margin_model/temporal_linear_margin_models.py

 from extreme_fit.model.margin_model.linear_margin_model.abstract_temporal_linear_margin_model import \
-    AbstractTemporalLinearMarginModel
+    AbstractTemporalLinearMarginModel, TemporalMarginFitMethod
 from extreme_fit.model.utils import r
 from extreme_fit.distribution.gev.gev_params import GevParams
+from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
 
 
 class StationaryTemporalModel(AbstractTemporalLinearMarginModel):
@@ -10,6 +11,15 @@ class StationaryTemporalModel(AbstractTemporalLinearMarginModel):
         super().load_margin_functions({})
 
 
+class GumbelTemporalModel(StationaryTemporalModel):
+
+    def __init__(self, coordinates: AbstractCoordinates, use_start_value=False, params_start_fit=None,
+                 params_sample=None, starting_point=None, fit_method=TemporalMarginFitMethod.is_mev_gev_fit,
+                 nb_iterations_for_bayesian_fit=5000, params_start_fit_bayesian=None):
+        super().__init__(coordinates, use_start_value, params_start_fit, params_sample, starting_point, fit_method,
+                         nb_iterations_for_bayesian_fit, params_start_fit_bayesian, type_for_MLE="Gumbel")
+
+
 class NonStationaryLocationTemporalModel(AbstractTemporalLinearMarginModel):
 
     def load_margin_functions(self, gev_param_name_to_dims=None):
@@ -62,4 +72,4 @@ class NonStationaryLocationAndScaleTemporalModel(AbstractTemporalLinearMarginMod
 
     @property
     def sigl(self):
-        return 1
+        return 1
\ No newline at end of file

extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_mle_extremes.py

 import numpy as np
 import rpy2
+from rpy2 import robjects
 
 from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_result_from_extremes import \
     AbstractResultFromExtremes
@@ -11,12 +12,17 @@ from extreme_fit.model.utils import r
 
 class ResultFromMleExtremes(AbstractResultFromExtremes):
 
+    def __init__(self, result_from_fit: robjects.ListVector, gev_param_name_to_dim=None,
+                 type_for_mle="GEV") -> None:
+        super().__init__(result_from_fit, gev_param_name_to_dim)
+        self.type_for_mle = type_for_mle
+
     @property
     def margin_coef_ordered_dict(self):
         values = self.name_to_value['results']
         d = self.get_python_dictionary(values)
         values = {i: param for i, param in enumerate(np.array(d['par']))}
-        return get_margin_coef_ordered_dict(self.gev_param_name_to_dim, values)
+        return get_margin_coef_ordered_dict(self.gev_param_name_to_dim, values, self.type_for_mle)
 
     def _confidence_interval_method(self, common_kwargs, ci_method, return_period):
         method_name = ci_method_to_method_name[ci_method]

extreme_fit/model/result_from_model_fit/utils.py

@@ -11,15 +11,19 @@ def convertFloatVector_to_float(f):
     return np.array(f)[0]
 
 
-def get_margin_coef_ordered_dict(gev_param_name_to_dim, mle_values):
+def get_margin_coef_ordered_dict(gev_param_name_to_dim, mle_values, type_for_mle="GEV"):
     assert gev_param_name_to_dim is not None
     # Build the Coeff dict from gev_param_name_to_dim
     coef_dict = OrderedDict()
     i = 0
     for gev_param_name in GevParams.PARAM_NAMES:
-        # Add intercept
+        # Add intercept (i.e. stationary parameter)
         intercept_coef_name = LinearCoef.coef_template_str(gev_param_name, LinearCoef.INTERCEPT_NAME).format(1)
-        coef_dict[intercept_coef_name] = mle_values[i]
+        if type_for_mle == "Gumbel" and i == 2:
+            coef_value = 0
+        else:
+            coef_value = mle_values[i]
+        coef_dict[intercept_coef_name] = coef_value
         i += 1
         # Add a potential linear temporal trend
         if gev_param_name in gev_param_name_to_dim:

test/test_extreme_fit/test_distribution/test_gev/test_gev_temporal_extremes_gumbel.py

+import unittest
+
+import numpy as np
+import pandas as pd
+
+from experiment.trend_analysis.univariate_test.abstract_gev_trend_test import fitted_linear_margin_estimator
+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, \
+    NonStationaryLocationTemporalModel, NonStationaryLocationAndScaleTemporalModel, GumbelTemporalModel
+from extreme_fit.model.utils import r, set_seed_r
+from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
+from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_coordinates import \
+    AbstractTemporalCoordinates
+from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
+from spatio_temporal_dataset.spatio_temporal_observations.abstract_spatio_temporal_observations import \
+    AbstractSpatioTemporalObservations
+
+
+class TestGevTemporalExtremesGumbel(unittest.TestCase):
+
+    def setUp(self) -> None:
+        set_seed_r()
+        r("""
+        N <- 50
+        loc = 0; scale = 1; shape <- 0
+        x_gev <- rgev(N, loc = loc, scale = scale, shape = shape)
+        start_loc = 0; start_scale = 1; start_shape = 1
+        """)
+        # Compute the stationary temporal margin with isMev
+        self.start_year = 0
+        df = pd.DataFrame({AbstractCoordinates.COORDINATE_T: range(self.start_year, self.start_year + 50)})
+        self.coordinates = AbstractTemporalCoordinates.from_df(df)
+        df2 = pd.DataFrame(data=np.array(r['x_gev']), index=df.index)
+        observations = AbstractSpatioTemporalObservations(df_maxima_gev=df2)
+        self.dataset = AbstractDataset(observations=observations, coordinates=self.coordinates)
+
+    def test_gev_temporal_margin_fit(self):
+        # Create estimator
+        estimator = fitted_linear_margin_estimator(GumbelTemporalModel,
+                                                   self.coordinates, self.dataset,
+                                                   starting_year=0,
+                                                   fit_method=TemporalMarginFitMethod.extremes_fevd_mle)
+        ref = {'loc': -0.0862185692806497, 'scale': 1.0818465357627252, 'shape': 0}
+        for year in range(1, 3):
+            mle_params_estimated = estimator.margin_function_from_fit.get_gev_params(np.array([year])).to_dict()
+            for key in ref.keys():
+                self.assertAlmostEqual(ref[key], mle_params_estimated[key], places=3)
+
+if __name__ == '__main__':
+    unittest.main()