[Thesis committee drawing] improve GEV non stationary plots. rename ci_delta...

[Thesis committee drawing] improve GEV non stationary plots.  rename ci_delta to ci_mle. catch some ci runtime exceptions. remove hatch/stripes drawing temporarily. fix unit for snow load. add new crocus classes & variable classes for snow load. modify conversion factor only for the new plots.

experiment/eurocode_data/eurocode_visualizer.py

 from typing import Dict, List, Tuple
 
 import matplotlib.pyplot as plt
+import numpy as np
 
 from extreme_fit.model.result_from_model_fit.result_from_extremes.eurocode_return_level_uncertainties import \
     EurocodeConfidenceIntervalFromExtremes
@@ -17,7 +18,7 @@ def get_label_name(model_name, ci_method_name: str):
     model_name += model_symbol
     model_name += '}}'
     model_name += parameter
-    model_name += ')_{ \\textrm{' + ci_method_name.upper() + '}} $ '
+    model_name += ')_{ \\textrm{' + ci_method_name.upper().split(' ')[1] + '}} $ '
     return model_name
 
 
@@ -25,7 +26,7 @@ def get_model_name(model_class):
     return get_display_name_from_object_type(model_class).split('Stationary')[0] + 'Stationary'
 
 
-def plot_massif_name_to_model_name_to_uncertainty_method_to_ordered_dict(d, nb_massif_names, nb_model_names, show=True):
+def plot_massif_name_to_model_name_to_uncertainty_method_to_ordered_dict(d, nb_massif_names, nb_model_names):
     """
     Rows correspond to massif names
     Columns correspond to stationary/non stationary model name for a given date
@@ -39,11 +40,7 @@ def plot_massif_name_to_model_name_to_uncertainty_method_to_ordered_dict(d, nb_m
         plot_model_name_to_uncertainty_method_to_ordered_dict(model_name_to_uncertainty_level,
                                                               massif_name, ax)
 
-    plt.suptitle('50-year return levels of extreme snow loads in France for several confiance interval methods.')
-
-    if show:
-        plt.show()
-
+    # plt.suptitle('50-year return levels of extreme snow loads in France for several confiance interval methods.')
 
 def plot_model_name_to_uncertainty_method_to_ordered_dict(d, massif_name, axes):
     if len(d) == 1:
@@ -71,17 +68,30 @@ def plot_label_to_ordered_return_level_uncertainties(ax, massif_name, model_name
         # Plot eurocode standards only for the first loop
         if j == 0:
             eurocode_region.plot_max_loading(ax, altitudes=altitudes)
-        conversion_factor = 100  # We divide by 100 to transform the snow water equivalent into snow load
-        mean = [r.mean_estimate / conversion_factor for r in ordered_return_level_uncertaines]
+        conversion_factor = 9.8 / 1000
+        mean = [r.mean_estimate * conversion_factor for r in ordered_return_level_uncertaines]
+        # Filter and keep only non nan values
+        not_nan_index = [not np.isnan(m) for m in mean]
+        mean = list(np.array(mean)[not_nan_index])
+        altitudes = list(np.array(altitudes)[not_nan_index])
+        ordered_return_level_uncertaines = list(np.array(ordered_return_level_uncertaines)[not_nan_index])
+
         ci_method_name = str(label).split('.')[1].replace('_', ' ')
-        ax.plot(altitudes, mean, '-', color=color, label=get_label_name(model_name, ci_method_name))
-        lower_bound = [r.confidence_interval[0] / conversion_factor for r in ordered_return_level_uncertaines]
-        upper_bound = [r.confidence_interval[1] / conversion_factor for r in ordered_return_level_uncertaines]
+        ax.plot(altitudes, mean, linestyle='--', marker='o', color=color, label=get_label_name(model_name, ci_method_name))
+        lower_bound = [r.confidence_interval[0] * conversion_factor for r in ordered_return_level_uncertaines]
+        upper_bound = [r.confidence_interval[1] * conversion_factor for r in ordered_return_level_uncertaines]
         ax.fill_between(altitudes, lower_bound, upper_bound, color=color, alpha=alpha)
     ax.legend(loc=2)
-    ax.set_ylim([0.0, 4.0])
-    ax.set_title(massif_name + ' ' + model_name)
-    ax.set_ylabel('50-year return level (N $m^-2$)')
+    ax.set_ylim([0.0, 16])
+    massif_name_str = massif_name.replace('_', ' ')
+    eurocode_region_str = get_display_name_from_object_type(type(eurocode_region))
+    if 'Non' in model_name:
+        model_name = 'non-stationary'
+    else:
+        model_name = 'stationary'
+    title = '{} ({} Eurocodes area) with a {} model'.format(massif_name_str, eurocode_region_str, model_name)
+    ax.set_title(title)
+    ax.set_ylabel('50-year return level (kN $m^-2$)')
     ax.set_xlabel('Altitude (m)')
     ax.grid()
 

experiment/eurocode_data/main_eurocode_drawing.py

 import time
+import os.path as op
+import matplotlib.pyplot as plt
 from collections import OrderedDict
 
+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
 from extreme_fit.model.result_from_model_fit.result_from_extremes.eurocode_return_level_uncertainties import \
     ConfidenceIntervalMethodFromExtremes
 from experiment.eurocode_data.eurocode_visualizer import \
     plot_massif_name_to_model_name_to_uncertainty_method_to_ordered_dict, get_model_name
 from experiment.eurocode_data.massif_name_to_departement import MASSIF_NAMES_ALPS
 from experiment.eurocode_data.utils import EUROCODE_ALTITUDES, LAST_YEAR_FOR_EUROCODE
-from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSwe3Days
+from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSwe3Days, CrocusSweTotal
 from experiment.meteo_france_data.scm_models_data.visualization.hypercube_visualization.altitude_hypercube_visualizer import \
     AltitudeHypercubeVisualizer
 from experiment.meteo_france_data.scm_models_data.visualization.hypercube_visualization.utils_hypercube import \
@@ -18,6 +23,8 @@ from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_m
 # Model class
 import matplotlib as mpl
 
+from root_utils import VERSION_TIME
+
 mpl.rcParams['text.usetex'] = True
 mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
 
@@ -32,7 +39,7 @@ def massif_name_to_ordered_return_level_uncertainties(model_class, last_year_for
                                                    only_first_one=False, save_to_file=False,
                                                    exact_starting_year=1958,
                                                    first_starting_year=None,
-                                                   study_classes=[CrocusSwe3Days],
+                                                   study_classes=[CrocusSwe3Days, CrocusSweTotal][1:],
                                                    trend_test_class=None)  # type: AltitudeHypercubeVisualizer
     # Loop on the data
     assert isinstance(altitude_visualizer.tuple_to_study_visualizer, OrderedDict)
@@ -51,29 +58,7 @@ def massif_name_to_ordered_return_level_uncertainties(model_class, last_year_for
     return {model_name: massif_name_to_ordered_eurocode_level_uncertainty}
 
 
-def main_drawing():
-    fast_plot = [True, False][0]
-    temporal_covariate = 2017
-    # Select parameters
-    massif_names = MASSIF_NAMES_ALPS[:]
-    model_class_and_last_year = [
-                                    (StationaryTemporalModel, 2017),
-                                    (NonStationaryLocationAndScaleTemporalModel, 2017),
-                                    # Add the temperature here
-                                ][1:]
-    altitudes = EUROCODE_ALTITUDES[:]
-    uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.my_bayes,
-                           ConfidenceIntervalMethodFromExtremes.ci_normal]
-    show = True
-
-    if fast_plot:
-        show = True
-        model_class_and_last_year = model_class_and_last_year[:1]
-        altitudes = altitudes[2:]
-        # altitudes = altitudes[:]
-        massif_names = massif_names[:1]
-        uncertainty_methods = uncertainty_methods[1:]
-
+def plot_ci_graphs(altitudes, massif_names, model_class_and_last_year, show, temporal_covariate, uncertainty_methods):
     model_name_to_massif_name_to_ordered_return_level = {}
     for model_class, last_year_for_the_data in model_class_and_last_year:
         start = time.time()
@@ -91,8 +76,80 @@ def main_drawing():
     # Plot graph
     plot_massif_name_to_model_name_to_uncertainty_method_to_ordered_dict(
         massif_name_to_model_name_to_ordered_return_level_uncertainties, nb_massif_names=len(massif_names),
-        nb_model_names=len(model_class_and_last_year), show=show)
+        nb_model_names=len(model_class_and_last_year))
+    if show:
+        plt.show()
+    else:
+        massif_names_str = '_'.join(massif_names)
+        model_names_str = '_'.join(
+            [model_name for model_name in model_name_to_massif_name_to_ordered_return_level.keys()])
+        filename = op.join(VERSION_TIME, model_names_str + '_' + massif_names_str)
+        StudyVisualizer.savefig_in_results(filename)
+
+
+def main_drawing():
+    fast_plot = [True, False][0]
+    temporal_covariate = 2017
+    # Select parameters
+    massif_names = MASSIF_NAMES_ALPS[:]
+    model_class_and_last_year = [
+                                    (StationaryTemporalModel, 2017),
+                                    (NonStationaryLocationAndScaleTemporalModel, 2017),
+                                    # Add the temperature here
+                                ][:]
+    altitudes = EUROCODE_ALTITUDES[:]
+    uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.my_bayes,
+                           ConfidenceIntervalMethodFromExtremes.ci_mle]
+    show = False
+
+    if fast_plot:
+        show = True
+        model_class_and_last_year = model_class_and_last_year[:]
+        altitudes = altitudes[-2:]
+        # altitudes = altitudes[:]
+        massif_names = massif_names[:1]
+        uncertainty_methods = uncertainty_methods[:]
+
+    plot_ci_graphs(altitudes, massif_names, model_class_and_last_year, show, temporal_covariate, uncertainty_methods)
+
+
+# Create 5 main plots
+def main_5_drawings():
+    model_class_and_last_year = [
+                                    (StationaryTemporalModel, 2017),
+                                    (NonStationaryLocationAndScaleTemporalModel, 2017),
+                                    # Add the temperature here
+                                ][:1]
+    altitudes = EUROCODE_ALTITUDES[:]
+    uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.my_bayes,
+                           ConfidenceIntervalMethodFromExtremes.ci_mle]
+    show = False
+    massif_names = MASSIF_NAMES_ALPS[:]
+    temporal_covariate = 2017
+    m = 4
+    n = (23 // m) + 1
+    for i in list(range(n))[:]:
+        massif_names = MASSIF_NAMES_ALPS[m * i: m * (i+1)]
+        print(massif_names)
+        plot_ci_graphs(altitudes, massif_names, model_class_and_last_year, show, temporal_covariate, uncertainty_methods)
+
+
+def main_3_massif_of_interest():
+    massif_names = ['Parpaillon', 'Chartreuse', 'Maurienne'][1:]
+    model_class_and_last_year = [
+                                    (StationaryTemporalModel, 2017),
+                                    (NonStationaryLocationAndScaleTemporalModel, 2017),
+                                    # Add the temperature here
+                                ][:]
+    altitudes = EUROCODE_ALTITUDES[1:]
+    uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.my_bayes,
+                           ConfidenceIntervalMethodFromExtremes.ci_mle][:]
+    temporal_covariate = 2017
+    show = False
+    plot_ci_graphs(altitudes, massif_names, model_class_and_last_year, show, temporal_covariate, uncertainty_methods)
 
 
 if __name__ == '__main__':
-    main_drawing()
+    # main_drawing()
+    # main_5_drawings()
+    main_3_massif_of_interest()
\ No newline at end of file

experiment/eurocode_data/slide_plot.py

+import matplotlib.pyplot as plt
+import numpy as np
+
+from experiment.eurocode_data.eurocode_region import C2, E
+from root_utils import get_display_name_from_object_type
+
+if __name__ == '__main__':
+    ax = plt.gca()
+    altitudes = np.linspace(200, 2000)
+    for region_class in [C2, E][1:]:
+        region_object = region_class()
+        region_object.plot_max_loading(ax, altitudes)
+        # ax.set_title(get_display_name_from_object_type(region_object) + ' Eurocodes region')
+        ax.set_ylabel('50-year return level (kN $m^-2$)')
+        ax.set_xlabel('Altitude (m)')
+        ax.set_ylim([0.0, 11.0])
+        ax.grid()
+        plt.show()
\ No newline at end of file

experiment/eurocode_data/utils.py

 
 # Eurocode quantile correspond to a 50 year return period
 EUROCODE_QUANTILE = 0.98
-# Altitudes (between low and mid altitudes) < 2000m
-EUROCODE_ALTITUDES = [0, 300, 600, 900, 1200, 1500, 1800]
+# Altitudes (between low and mid altitudes) < 2000m and should be > 200m
+EUROCODE_ALTITUDES = [300, 600, 900, 1200, 1500, 1800]
 #  Last year taken into account for the Eurocode
 # Date of publication was 2014, therefore the winter 2013/2014 could not have been measured
 # Therefore, the winter 2012/2013 was the last one. Thus, 2012 is the last year for the Eurocode

experiment/meteo_france_data/scm_models_data/abstract_study.py

@@ -335,15 +335,16 @@ class AbstractStudy(object):
                 else:
                     ax.fill(*coords_list, **{'color': default_color_for_missing_massif})
 
-                # Add a hatch to visualize the mean & variance variation sign
-                hatch_list = ['//', '\\\\']
-                if massif_name_to_hatch_boolean_list is not None:
-                    if massif_name in massif_name_to_hatch_boolean_list:
-                        a = np.array(coords_list).transpose()
-                        hatch_boolean_list = massif_name_to_hatch_boolean_list[massif_name]
-                        for hatch, is_hatch in zip(hatch_list, hatch_boolean_list):
-                            if is_hatch:
-                                ax.add_patch(Polygon(xy=a, fill=False, hatch=hatch))
+                # For the moment we comment all the part of this code
+                # # Add a hatch to visualize the mean & variance variation sign
+                # hatch_list = ['//', '\\\\']
+                # if massif_name_to_hatch_boolean_list is not None:
+                #     if massif_name in massif_name_to_hatch_boolean_list:
+                #         a = np.array(coords_list).transpose()
+                #         hatch_boolean_list = massif_name_to_hatch_boolean_list[massif_name]
+                #         for hatch, is_hatch in zip(hatch_list, hatch_boolean_list):
+                #             if is_hatch:
+                #                 ax.add_patch(Polygon(xy=a, fill=False, hatch=hatch))
 
         # Display the center of the massif
         masssif_coordinate_for_display = cls.massifs_coordinates_for_display(massif_names)
@@ -404,9 +405,9 @@ class AbstractStudy(object):
     def map_full_path(self) -> str:
         return op.join(self.full_path, 'map')
 
-    @property
-    def result_full_path(self) -> str:
-        return op.join(self.full_path, 'results')
+    @classproperty
+    def result_full_path(cls) -> str:
+        return op.join(cls.full_path, 'results')
 
     @property
     def study_full_path(self) -> str:

experiment/meteo_france_data/scm_models_data/crocus/crocus.py

@@ -3,7 +3,7 @@ import numpy as np
 from experiment.meteo_france_data.scm_models_data.abstract_extended_study import AbstractExtendedStudy
 from experiment.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
 from experiment.meteo_france_data.scm_models_data.crocus.crocus_variables import CrocusTotalSweVariable, \
-    CrocusDepthVariable, CrocusRecentSweVariable
+    CrocusDepthVariable, CrocusRecentSweVariable, TotalSnowLoadVariable, RecentSnowLoadVariable
 
 
 class Crocus(AbstractStudy):
@@ -12,7 +12,8 @@ class Crocus(AbstractStudy):
     """
 
     def __init__(self, variable_class, *args, **kwargs):
-        assert variable_class in [CrocusTotalSweVariable, CrocusDepthVariable, CrocusRecentSweVariable]
+        assert variable_class in [CrocusTotalSweVariable, CrocusDepthVariable, CrocusRecentSweVariable,
+                                  RecentSnowLoadVariable, TotalSnowLoadVariable]
         super().__init__(variable_class, *args, **kwargs)
         self.model_name = 'Crocus'
 
@@ -44,6 +45,19 @@ class CrocusSweTotal(Crocus):
         return self.winter_annual_aggregation(time_serie)
 
 
+# Create some class that enables to deal directly with the snow load
+
+
+class CrocusSnowLoadTotal(Crocus):
+    def __init__(self, *args, **kwargs):
+        Crocus.__init__(self, TotalSnowLoadVariable, *args, **kwargs)
+
+
+class CrocusSnowLoad3Days(CrocusSweTotal):
+    def __init__(self, *args, **kwargs):
+        Crocus.__init__(self, RecentSnowLoadVariable, *args, **kwargs)
+
+
 class ExtendedCrocusSweTotal(AbstractExtendedStudy, CrocusSweTotal):
     pass
 

experiment/meteo_france_data/scm_models_data/crocus/crocus_variables.py

@@ -27,6 +27,22 @@ class CrocusRecentSweVariable(CrocusTotalSweVariable):
         return 'SWE_3DY_ISBA'
 
 
+class AbstractSnowLoadVariable(CrocusVariable):
+    UNIT = 'kN $m^{-2}$'
+    snow_load_multiplication_factor = 9.81 / 1000
+
+    @property
+    def daily_time_serie_array(self) -> np.ndarray:
+        return self.snow_load_multiplication_factor * super().daily_time_serie_array
+
+
+class RecentSnowLoadVariable(AbstractSnowLoadVariable, CrocusRecentSweVariable):
+    NAME = 'Snow load last 3 days'
+
+class TotalSnowLoadVariable(AbstractSnowLoadVariable, CrocusTotalSweVariable):
+    NAME = 'Snow load total'
+
+
 class CrocusDepthVariable(CrocusVariable):
     NAME = 'Snow Depth'
     UNIT = 'm'

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

@@ -6,7 +6,7 @@ from experiment.meteo_france_data.scm_models_data.visualization.study_visualizat
 from experiment.trend_analysis.abstract_score import MannKendall
 from experiment.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
 from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusDepth, CrocusSweTotal, ExtendedCrocusDepth, \
-    ExtendedCrocusSweTotal, CrocusDaysWithSnowOnGround, CrocusSwe3Days
+    ExtendedCrocusSweTotal, CrocusDaysWithSnowOnGround, CrocusSwe3Days, CrocusSnowLoad3Days, CrocusSnowLoadTotal
 from experiment.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall, ExtendedSafranSnowfall, \
     SafranRainfall, \
     SafranTemperature, SafranTotalPrecip
@@ -23,9 +23,11 @@ SCM_STUDIES_NAMES = [get_display_name_from_object_type(k) for k in SCM_STUDIES]
 SCM_STUDY_NAME_TO_SCM_STUDY = dict(zip(SCM_STUDIES_NAMES, SCM_STUDIES))
 SCM_STUDY_CLASS_TO_ABBREVIATION = {
     SafranSnowfall: 'SF3',
-    CrocusSweTotal: 'SWET',
+    CrocusSweTotal: 'SWE',
     CrocusSwe3Days: 'SWE3',
     CrocusDepth: 'SD',
+    CrocusSnowLoadTotal: 'SL',
+    CrocusSnowLoad3Days: 'SL3',
 }
 
 altitude_massif_name_and_study_class_for_poster = [
@@ -40,6 +42,12 @@ altitude_massif_name_and_study_class_for_poster_evan = [
     (2700, 'Parpaillon', CrocusSwe3Days),
 ]
 
+altitude_massif_name_and_study_class_for_committee= [
+    (900, 'Chartreuse', CrocusSnowLoad3Days),
+    (1800, 'Vanoise', CrocusSnowLoad3Days),
+    (2700, 'Parpaillon', CrocusSnowLoad3Days),
+]
+
 
 SCM_STUDY_NAME_TO_ABBREVIATION = {get_display_name_from_object_type(k): v for k, v in
                                   SCM_STUDY_CLASS_TO_ABBREVIATION.items()}
@@ -250,7 +258,8 @@ def max_graph_annual_maxima_poster():
     choice_tuple = [
         altitude_massif_name_and_study_class_for_poster,
         altitude_massif_name_and_study_class_for_poster_evan,
-    ][1]
+        altitude_massif_name_and_study_class_for_committee,
+    ][2]
     for altitude, massif_name, study_class in choice_tuple:
         for study in study_iterator_global([study_class], altitudes=[altitude]):
             study_visualizer = StudyVisualizer(study, save_to_file=save_to_file,

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

@@ -625,7 +625,10 @@ class StudyVisualizer(VisualizationParameters):
         ax = plt.gca()
         x, y = self.smooth_maxima_x_y(massif_names.index(massif_name))
         ax.plot(x, y, color=color, linewidth=5)
-        ax.set_ylabel('{} (in {})'.format(snow_abbreviation, self.study.variable_unit), color=color, fontsize=15)
+        # ax.set_ylabel('{} (in {})'.format(snow_abbreviation, self.study.variable_unit), color=color, fontsize=15)
+        ax.set_ylabel('{} (in {})'.format(snow_abbreviation, self.study.variable_unit), fontsize=15)
+        ax.set_xlabel('years', fontsize=15)
+        ax.set_title('{} at {} m'.format(massif_name, altitude))
         ax.xaxis.set_ticks(x[2::10])
         ax.tick_params(axis='both', which='major', labelsize=13)
 
@@ -824,11 +827,15 @@ class StudyVisualizer(VisualizationParameters):
             if not self.only_one_graph:
                 filename += "{}".format('_'.join(self.plot_name.split())) + '_'
             filename += specific_title
-            filepath = op.join(self.study.result_full_path, filename + '.png')
-            dirname = op.dirname(filepath)
-            if not op.exists(dirname):
-                os.makedirs(dirname, exist_ok=True)
-            plt.savefig(filepath)
+            self.savefig_in_results(filename)
+
+    @classmethod
+    def savefig_in_results(cls, filename):
+        filepath = op.join(AbstractStudy.result_full_path, filename + '.png')
+        dirname = op.dirname(filepath)
+        if not op.exists(dirname):
+            os.makedirs(dirname, exist_ok=True)
+        plt.savefig(filepath)
 
     def visualize_independent_margin_fits(self, threshold=None, axes=None, show=True):
         # Fit either a GEV or a GPD

experiment/paper1/poster_EVAN2019/main_poster_EVAN2019.py

@@ -17,7 +17,7 @@ mpl.rcParams['hatch.linewidth'] = 0.3
 
 def main_poster_A_non_stationary_model_choice():
     nb = 1
-    for altitude in POSTER_ALTITUDES[:nb]:
+    for altitude in POSTER_ALTITUDES[:]:
         for trend_test_class in [GevLocationTrendTest, GevScaleTrendTest, GevLocationAndScaleTrendTest][-nb:]:
             vizualiser = get_full_altitude_visualizer(Altitude_Hypercube_Year_Visualizer, altitude=altitude,
                                                       exact_starting_year=1958, reduce_strength_array=False,
@@ -133,7 +133,7 @@ def main_poster_D_other_quantities_analysis():
 
 
 if __name__ == '__main__':
-    # main_poster_A_non_stationary_model_choice()
-    main_poster_B_starting_years_analysis()
+    main_poster_A_non_stationary_model_choice()
+    # main_poster_B_starting_years_analysis()
     # main_poster_C_orientation_analysis()
     # main_poster_D_other_quantities_analysis()

extreme_fit/model/result_from_model_fit/result_from_extremes/confidence_interval_method.py

@@ -4,7 +4,7 @@ from enum import Enum
 class ConfidenceIntervalMethodFromExtremes(Enum):
     # Confidence interval from the ci function
     ci_bayes = 0
-    ci_normal = 1
+    ci_mle = 1
     ci_boot = 2
     ci_proflik = 3
     # Confidence interval from my functions
@@ -12,7 +12,7 @@ class ConfidenceIntervalMethodFromExtremes(Enum):
 
 
 ci_method_to_method_name = {
-    ConfidenceIntervalMethodFromExtremes.ci_normal: 'normal',
+    ConfidenceIntervalMethodFromExtremes.ci_mle: 'normal',
     ConfidenceIntervalMethodFromExtremes.ci_boot: 'boot',
     ConfidenceIntervalMethodFromExtremes.ci_proflik: 'proflik',
 }

extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_mle_extremes.py

 import numpy as np
+import rpy2
 
 from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_result_from_extremes import \
     AbstractResultFromExtremes
@@ -23,7 +24,10 @@ class ResultFromMleExtremes(AbstractResultFromExtremes):
                 'method': method_name,
             # xrange = NULL, nint = 20
         }
-        res = r.ci(self.result_from_fit, **mle_ci_parameters, **common_kwargs)
+        try:
+            res = r.ci(self.result_from_fit, **mle_ci_parameters, **common_kwargs)
+        except rpy2.rinterface.RRuntimeError:
+            return np.nan, (np.nan, np.nan)
         if self.is_non_stationary:
             a = np.array(res)[0]
             lower, mean_estimate, upper, _ = a

test/test_extreme_fit/test_model/test_confidence_interval.py

@@ -74,7 +74,7 @@ class TestConfidenceInterval(unittest.TestCase):
 
     def test_ci_normal(self):
         self.fit_method = TemporalMarginFitMethod.extremes_fevd_mle
-        self.ci_method = ConfidenceIntervalMethodFromExtremes.ci_normal
+        self.ci_method = ConfidenceIntervalMethodFromExtremes.ci_mle
         self.model_class_to_triplet= {
             StationaryTemporalModel: (-4.703945484843988, 30.482318639674023, 65.66858276419204),
             NonStationaryLocationTemporalModel: (-4.223086740397132, 30.29842988666537, 64.81994651372787),

thesis_report/gev_plot.py

@@ -38,14 +38,23 @@ def gev_plot_big():
 def gev_plot_big_non_stationary_location():
     lim = 5
     x = np.linspace(-lim, lim, 100)
-    scale, shape = 1, 1
-    locs = [-1, 1]
-    colors = ['red', 'green']
-    y = r.dgev(x, 0.0, scale, shape)
-    plt.plot(x, y, label='distribution of $Y(1968)$', linewidth=5)
+    scale, shape = 1, 0.0
+    locs = [1, 2, 3]
+    inverse_loc_with_scale = True
+    colors = ['red','k', 'green']
+    greek_leeter = ' $\{}_1'.format('mu' if not inverse_loc_with_scale else 'sigma')
+    plt.title('Density for the distribution of Y(t) with different{}$'.format(greek_leeter))
+    template = greek_leeter + '{} 0$'
     for loc, color in zip(locs, colors):
-        sign_str = '>' if loc > 0 else '<'
-        label = 'distribution of $Y(1969)$ with $\mu_1 {} 0$'.format(sign_str)
+        if loc == locs[1]:
+            sign_str = '='
+        elif loc == locs[2]:
+            sign_str = '>'
+        else:
+            sign_str = '<'
+        label = template.format(sign_str)
+        if inverse_loc_with_scale:
+            loc, scale = scale, loc
         y = r.dgev(x, loc, scale, shape)
         plt.plot(x, y, label=label, linewidth=5, color=color)
     plt.legend(prop={'size': 20})