diff --git a/extreme_data/meteo_france_data/scm_models_data/abstract_study.py b/extreme_data/meteo_france_data/scm_models_data/abstract_study.py
index e6de599ebf7334a5468d1e1faccd135d89e13b66..44228b3f7289fcd6a4aaf88cd83c886e510dca77 100644
--- a/extreme_data/meteo_france_data/scm_models_data/abstract_study.py
+++ b/extreme_data/meteo_france_data/scm_models_data/abstract_study.py
@@ -27,6 +27,7 @@ from extreme_data.meteo_france_data.scm_models_data.utils import ALTITUDES, ZS_I
first_day_and_last_day, FrenchRegion, ZS_INT_MASK_PYRENNES, alps_massif_order, ZS_INT_MASK_PYRENNES_LIST, \
season_to_str
from extreme_data.meteo_france_data.scm_models_data.visualization.utils import get_km_formatter
+from extreme_fit.estimator.margin_estimator.utils import fitted_stationary_gev
from extreme_fit.function.margin_function.abstract_margin_function import \
AbstractMarginFunction
from extreme_data.meteo_france_data.scm_models_data.visualization.create_shifted_cmap import create_colorbase_axis, \
@@ -47,6 +48,7 @@ filled_marker_legend_list2 = ['Filled marker = Selected', 'model is significant'
YEAR_MIN = 1959
YEAR_MAX = 2019
+
class AbstractStudy(object):
"""
A Study is defined by:
@@ -86,7 +88,7 @@ class AbstractStudy(object):
self.multiprocessing = multiprocessing
self.season = season
if split_years is None:
- split_years = list(range(year_min, year_max+1))
+ split_years = list(range(year_min, year_max + 1))
self.split_years = set(split_years)
# Add some attributes, for the "allslopes" reanalysis
assert orientation is None or orientation in ORIENTATIONS
@@ -793,3 +795,27 @@ class AbstractStudy(object):
mask_massif = np.array(img)
mask_french_alps += mask_massif
return ~np.array(mask_french_alps, dtype=bool)
+
+ @cached_property
+ def massif_name_to_stationary_gev_params(self):
+ massif_name_to_stationary_gev_params = {}
+ for massif_name, annual_maxima in self.massif_name_to_annual_maxima.items():
+ gev_params = fitted_stationary_gev(annual_maxima)
+ massif_name_to_stationary_gev_params[massif_name] = gev_params
+ return massif_name_to_stationary_gev_params
+
+ def massif_name_to_gev_param_list(self, year_min_and_max_list):
+ year_to_index = {y: i for i, y in enumerate(self.ordered_years)}
+ self.index_min_and_max_list = [(year_to_index[mi], year_to_index[ma]) for mi, ma in year_min_and_max_list]
+ return self._massif_name_to_gev_param_list
+
+ @cached_property
+ def _massif_name_to_gev_param_list(self):
+ massif_name_to_stationary_gev_params = {}
+ for massif_name, annual_maxima in self.massif_name_to_annual_maxima.items():
+ gev_params_list = []
+ for index_min, index_max in self.index_min_and_max_list:
+ annual_maxima_sliced = annual_maxima[index_min: index_max+1]
+ gev_params_list.append(fitted_stationary_gev(annual_maxima_sliced))
+ massif_name_to_stationary_gev_params[massif_name] = gev_params_list
+ return massif_name_to_stationary_gev_params
diff --git a/extreme_data/meteo_france_data/scm_models_data/visualization/plot_utils.py b/extreme_data/meteo_france_data/scm_models_data/visualization/plot_utils.py
index 525f0cfc4af99bd1cfdac4af1c83058276cd62ac..0c22e10d84c1ade073a3782be3f29c86903d98da 100644
--- a/extreme_data/meteo_france_data/scm_models_data/visualization/plot_utils.py
+++ b/extreme_data/meteo_france_data/scm_models_data/visualization/plot_utils.py
@@ -7,7 +7,7 @@ from extreme_data.meteo_france_data.scm_models_data.visualization.create_shifted
get_half_colormap
-def plot_against_altitude(altitudes, ax, massif_id, massif_name, values):
+def plot_against_altitude(altitudes, ax, massif_id, massif_name, values, fill=False):
di = massif_id // 8
if di == 0:
linestyle = '-'
@@ -19,7 +19,13 @@ def plot_against_altitude(altitudes, ax, massif_id, massif_name, values):
colors[-3:-1] = [] # remove gray and olive
color = colors[massif_id % 8]
massif_name_str = ' '.join(massif_name.split('_'))
- ax.plot(altitudes, values, color=color, linewidth=2, label=massif_name_str, linestyle=linestyle)
+ label = massif_name_str
+ if not fill:
+ ax.plot(altitudes, values, color=color, linewidth=2, label=label, linestyle=linestyle)
+ else:
+ lower_bound, upper_bound = zip(*values)
+ # ax.fill_between(altitudes, lower_bound, upper_bound, color=color, alpha=0.2, label=label + '95\% confidence interval')
+ ax.fill_between(altitudes, lower_bound, upper_bound, color=color, alpha=0.2)
def load_plot(cmap, graduation, label, massif_name_to_value, altitude, fit_method, add_x_label=True,
diff --git a/extreme_fit/distribution/gev/gev_params.py b/extreme_fit/distribution/gev/gev_params.py
index d73a76d64fec5711e3b6afeaf5ff9b4a3653997b..384509c31a5b7e8437c7dfab6e6f030060875e26 100644
--- a/extreme_fit/distribution/gev/gev_params.py
+++ b/extreme_fit/distribution/gev/gev_params.py
@@ -24,6 +24,7 @@ class GevParams(AbstractExtremeParams):
self.block_size = block_size
if accept_zero_scale_parameter and scale == 0.0:
self.has_undefined_parameters = False
+ self.param_name_to_confidence_interval = None
def sample(self, n) -> np.ndarray:
if self.has_undefined_parameters:
diff --git a/extreme_fit/estimator/margin_estimator/utils.py b/extreme_fit/estimator/margin_estimator/utils.py
index 75f4315e59a262029a91ea76ccfcd2d2d9fb1232..7d5962e70319df8c886ffa2540ade8104ebccaed 100644
--- a/extreme_fit/estimator/margin_estimator/utils.py
+++ b/extreme_fit/estimator/margin_estimator/utils.py
@@ -36,6 +36,18 @@ def fitted_stationary_gev(x_gev, fit_method=MarginFitMethod.is_mev_gev_fit, mode
estimator = fitted_linear_margin_estimator(model_class, coordinates, dataset, starting_year, fit_method)
first_coordinate = coordinates.coordinates_values()[0]
gev_param = estimator.function_from_fit.get_params(first_coordinate)
+
+ # Build confidence interval
+ fit_method_with_confidence_interval_implemented = [MarginFitMethod.is_mev_gev_fit]
+ if fit_method in fit_method_with_confidence_interval_implemented:
+ param_name_to_confidence_interval = {}
+ for i, (param_name, param_value) in enumerate(gev_param.to_dict().items()):
+ confidence_interval_half_size = estimator.result_from_model_fit.confidence_interval_half_sizes[i]
+ confidence_interval = (param_value - confidence_interval_half_size, param_value + confidence_interval_half_size)
+ param_name_to_confidence_interval[param_name] = confidence_interval
+ gev_param.param_name_to_confidence_interval = param_name_to_confidence_interval
+
+ # Warning
if not -0.5 < gev_param.shape < 0.5:
warnings.warn('fitted shape parameter is outside physical bounds {}'.format(gev_param.shape))
return gev_param
diff --git a/extreme_fit/model/result_from_model_fit/abstract_result_from_model_fit.py b/extreme_fit/model/result_from_model_fit/abstract_result_from_model_fit.py
index dcc8bed991fa972337639af63d49b57ed5b2c840..ced1aa90415b2f45d50436262dfc98e924199066 100644
--- a/extreme_fit/model/result_from_model_fit/abstract_result_from_model_fit.py
+++ b/extreme_fit/model/result_from_model_fit/abstract_result_from_model_fit.py
@@ -1,5 +1,7 @@
+import numpy as np
from rpy2 import robjects
+
class AbstractResultFromModelFit(object):
def __init__(self, result_from_fit: robjects.ListVector) -> None:
@@ -9,6 +11,19 @@ class AbstractResultFromModelFit(object):
else:
self.name_to_value = {}
+ @property
+ def variance_covariance_matrix(self):
+ raise NotImplementedError
+
+ @property
+ def standard_errors_for_mle(self):
+ """ See Coles 2001 page 41, for an example"""
+ return np.sqrt(np.diagonal(self.variance_covariance_matrix))
+
+ @property
+ def confidence_interval_half_sizes(self):
+ return [1.96 * s for s in self.standard_errors_for_mle]
+
@staticmethod
def get_python_dictionary(r_dictionary):
return {name: r_dictionary.rx2(name) for name in r_dictionary.names}
@@ -56,8 +71,3 @@ class AbstractResultFromModelFit(object):
@property
def covariance(self):
raise NotImplementedError
-
-
-
-
-
diff --git a/extreme_fit/model/result_from_model_fit/result_from_ismev.py b/extreme_fit/model/result_from_model_fit/result_from_ismev.py
index 8db506758817ac6ee46ea489333b6dc5d05e44f3..0e424454526505cf2108569636fe1c53412e4198 100644
--- a/extreme_fit/model/result_from_model_fit/result_from_ismev.py
+++ b/extreme_fit/model/result_from_model_fit/result_from_ismev.py
@@ -14,6 +14,10 @@ class ResultFromIsmev(AbstractResultFromModelFit):
super().__init__(result_from_fit)
self.param_name_to_dim = param_name_to_dim
+ @property
+ def variance_covariance_matrix(self):
+ return np.array(self.name_to_value['cov'])
+
@property
def margin_coef_ordered_dict(self):
return get_margin_coef_ordered_dict(self.param_name_to_dim, self.name_to_value['mle'])
diff --git a/projects/contrasting_trends_in_snow_loads/post thesis committee/pointwise_gev_study_visualizer.py b/projects/contrasting_trends_in_snow_loads/post thesis committee/pointwise_gev_study_visualizer.py
new file mode 100644
index 0000000000000000000000000000000000000000..4549b0ded27dd233a9f30eee9d9dd6805cb08583
--- /dev/null
+++ b/projects/contrasting_trends_in_snow_loads/post thesis committee/pointwise_gev_study_visualizer.py
@@ -0,0 +1,95 @@
+import matplotlib.pyplot as plt
+from cached_property import cached_property
+
+from extreme_data.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
+from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day
+from extreme_data.meteo_france_data.scm_models_data.visualization.plot_utils import plot_against_altitude
+from extreme_data.meteo_france_data.scm_models_data.visualization.study_visualizer import StudyVisualizer
+from extreme_fit.distribution.gev.gev_params import GevParams
+from projects.altitude_spatial_model.altitudes_fit.altitudes_studies import AltitudesStudies
+from projects.exceeding_snow_loads.utils import paper_altitudes
+
+
+class PointwiseGevStudyVisualizer(AltitudesStudies):
+
+ def plot_gev_params_against_altitude(self):
+ for param_name in GevParams.PARAM_NAMES[:]:
+ ax = plt.gca()
+ for massif_name in self.study.all_massif_names()[:]:
+ self._plot_gev_params_against_altitude_one_massif(ax, massif_name, param_name)
+ ax.legend(prop={'size': 7}, ncol=3)
+ ax.set_xlabel('Altitude')
+ ax.set_ylabel(param_name)
+ plot_name = '{} change /with altitude'.format(param_name)
+ self.show_or_save_to_file(plot_name, no_title=True, tight_layout=True, show=False)
+ ax.clear()
+ plt.close()
+
+ def _plot_gev_params_against_altitude_one_massif(self, ax, massif_name, param_name):
+ altitudes = []
+ params = []
+ confidence_intervals = []
+ for altitude, study in self.altitude_to_study.items():
+ if massif_name in study.study_massif_names:
+ altitudes.append(altitude)
+ gev_params = study.massif_name_to_stationary_gev_params[massif_name]
+ params.append(gev_params.to_dict()[param_name])
+ confidence_intervals.append(gev_params.param_name_to_confidence_interval[param_name])
+ massif_id = self.study.all_massif_names().index(massif_name)
+ plot_against_altitude(altitudes, ax, massif_id, massif_name, params, fill=False)
+ # plot_against_altitude(altitudes, ax, massif_id, massif_name, confidence_intervals, fill=True)
+
+ # Plot against the time
+
+ @property
+ def year_min_and_max_list(self):
+ l = []
+ year_min, year_max = 1959, 1989
+ for shift in range(0, 7):
+ l.append((year_min + 5 * shift, year_max + 5 * shift))
+ return l
+
+ @property
+ def min_years_for_plot_x_axis(self):
+ return [c[0] for c in self.year_min_and_max_list]
+
+ def plot_gev_params_against_time(self):
+ for altitude in self.altitudes[:]:
+ self._plot_gev_params_against_time(altitude)
+
+ def _plot_gev_params_against_time(self, altitude):
+ for param_name in GevParams.PARAM_NAMES[:]:
+ ax = plt.gca()
+ for massif_name in self.study.all_massif_names()[:]:
+ self._plot_gev_params_against_time_one_massif(ax, massif_name, param_name, altitude)
+ ax.legend(prop={'size': 7}, ncol=3)
+ ax.set_xlabel('Year')
+ ax.set_ylabel(param_name + ' for altitude={}'.format(altitude))
+ xlabels = ['-'.join([str(e) for e in t]) for t in self.year_min_and_max_list]
+ ax.set_xticks(self.min_years_for_plot_x_axis)
+ ax.set_xticklabels(xlabels)
+ # ax.tick_params(labelsize=5)
+ plot_name = '{} change /with years for altitude={}'.format(param_name, altitude)
+ self.show_or_save_to_file(plot_name, no_title=True, tight_layout=True, show=False)
+ ax.clear()
+ plt.close()
+
+ def _plot_gev_params_against_time_one_massif(self, ax, massif_name, param_name, altitude):
+ study = self.altitude_to_study[altitude]
+ if massif_name in study.study_massif_names:
+ gev_params_list = study.massif_name_to_gev_param_list(self.year_min_and_max_list)[massif_name]
+ params = [gev_params.to_dict()[param_name] for gev_params in gev_params_list]
+ confidence_intervals = [gev_params.param_name_to_confidence_interval[param_name] for gev_params in
+ gev_params_list]
+ massif_id = self.study.all_massif_names().index(massif_name)
+ plot_against_altitude(self.min_years_for_plot_x_axis, ax, massif_id, massif_name, params, fill=False)
+ # plot_against_altitude(self.years, ax, massif_id, massif_name, confidence_intervals, fill=True)
+
+
+if __name__ == '__main__':
+ altitudes = [600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600]
+ # altitudes = paper_altitudes
+ # altitudes = [1800, 2100]
+ visualizer = PointwiseGevStudyVisualizer(SafranSnowfall1Day, altitudes=altitudes)
+ visualizer.plot_gev_params_against_altitude()
+ visualizer.plot_gev_params_against_time()