From 68560237d0f963899e58fbdef7f7e29b21f6cf8c Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Mon, 27 Jul 2020 20:42:12 +0200
Subject: [PATCH] [exceeding] add plot year for the peak. add many altitudinal
models.
---
.../gev_altitudinal_models.py | 17 +++++-
..._altitudinal_models_cross_term_in_scale.py | 45 +++++++++++++++
.../gev_altitudinal_models_with_scale.py | 57 +++++++++++++++++++
.../polynomial_margin_model/utils.py | 31 ++++++++++
.../altitudes_fit/main_altitudes_studies.py | 13 ++---
...es_visualizer_for_non_stationary_models.py | 34 +++++++++--
.../one_fold_analysis/one_fold_fit.py | 4 +-
.../one_fold_analysis/plot_total_aic.py | 15 +++--
8 files changed, 197 insertions(+), 19 deletions(-)
create mode 100644 extreme_fit/model/margin_model/polynomial_margin_model/gev_altitudinal_models_cross_term_in_scale.py
create mode 100644 extreme_fit/model/margin_model/polynomial_margin_model/gev_altitudinal_models_with_scale.py
diff --git a/extreme_fit/model/margin_model/polynomial_margin_model/gev_altitudinal_models.py b/extreme_fit/model/margin_model/polynomial_margin_model/gev_altitudinal_models.py
index 5b67925c..218b2780 100644
--- a/extreme_fit/model/margin_model/polynomial_margin_model/gev_altitudinal_models.py
+++ b/extreme_fit/model/margin_model/polynomial_margin_model/gev_altitudinal_models.py
@@ -36,7 +36,9 @@ class AbstractAltitudinalModel(AbstractSpatioTemporalPolynomialModel):
name += self.dim_to_str_number(GevParams.LOC, self.coordinates.idx_temporal_coordinates)
name += self.dim_to_str_number(GevParams.SCALE, self.coordinates.idx_temporal_coordinates)
if isinstance(self, AbstractAddCrossTermForLocation):
- name += 'x'
+ name += 'l'
+ if isinstance(self, AbstractAddCrossTermForScale):
+ name += 's'
return name
@@ -112,6 +114,19 @@ class AbstractAddCrossTermForLocation(AbstractAltitudinalModel):
d[GevParams.LOC].insert(1, ((self.coordinates.idx_x_coordinates, self.coordinates.idx_temporal_coordinates), 1))
return d
+class AbstractAddCrossTermForScale(AbstractAltitudinalModel):
+
+ @property
+ def param_name_to_list_dim_and_degree_for_margin_function(self):
+ d = self.param_name_to_list_dim_and_degree
+ # The two insert below enable to check that the insert_index should indeed be 1
+ assert 1 <= len(d[GevParams.SCALE]) <= 2
+ assert self.coordinates.idx_x_coordinates == d[GevParams.SCALE][0][0]
+ insert_index = 1
+ d[GevParams.SCALE].insert(insert_index, ((self.coordinates.idx_x_coordinates, self.coordinates.idx_temporal_coordinates), 1))
+ return d
+
+
class NonStationaryCrossTermForLocation(AbstractAddCrossTermForLocation, StationaryAltitudinal):
pass
diff --git a/extreme_fit/model/margin_model/polynomial_margin_model/gev_altitudinal_models_cross_term_in_scale.py b/extreme_fit/model/margin_model/polynomial_margin_model/gev_altitudinal_models_cross_term_in_scale.py
new file mode 100644
index 00000000..77ce775b
--- /dev/null
+++ b/extreme_fit/model/margin_model/polynomial_margin_model/gev_altitudinal_models_cross_term_in_scale.py
@@ -0,0 +1,45 @@
+from extreme_fit.distribution.gev.gev_params import GevParams
+from extreme_fit.model.margin_model.polynomial_margin_model.gev_altitudinal_models import AbstractAltitudinalModel, \
+ NonStationaryAltitudinalScaleLinear, NonStationaryAltitudinalLocationLinearScaleLinear, \
+ NonStationaryAltitudinalLocationQuadraticScaleLinear, AbstractAddCrossTermForScale
+from extreme_fit.model.margin_model.polynomial_margin_model.gev_altitudinal_models_with_scale import \
+ NonStationaryAltitudinalLocationQuadraticScaleQuadratic, NonStationaryAltitudinalScaleQuadratic, \
+ NonStationaryAltitudinalLocationLinearScaleQuadratic
+from extreme_fit.model.margin_model.polynomial_margin_model.polynomial_margin_model import PolynomialMarginModel
+from extreme_fit.model.margin_model.polynomial_margin_model.spatio_temporal_polynomial_model import \
+ AbstractSpatioTemporalPolynomialModel
+
+
+
+
+class NonStationaryAltitudinalScaleLinearCrossTermForScale(AbstractAddCrossTermForScale,
+ NonStationaryAltitudinalScaleLinear):
+ pass
+
+
+class NonStationaryAltitudinalLocationLinearScaleLinearCrossTermForScale(AbstractAddCrossTermForScale,
+ NonStationaryAltitudinalLocationLinearScaleLinear):
+ pass
+
+
+class NonStationaryAltitudinalLocationQuadraticScaleLinearCrossTermForScale(AbstractAddCrossTermForScale,
+ NonStationaryAltitudinalLocationQuadraticScaleLinear):
+ pass
+
+class NonStationaryAltitudinalScaleQuadraticCrossTermForScale(AbstractAddCrossTermForScale, NonStationaryAltitudinalScaleQuadratic):
+ pass
+
+
+class NonStationaryAltitudinalLocationLinearScaleQuadraticCrossTermForScale(AbstractAddCrossTermForScale,
+ NonStationaryAltitudinalLocationLinearScaleQuadratic):
+ pass
+
+
+class NonStationaryAltitudinalLocationQuadraticScaleQuadraticCrossTermForScale(AbstractAddCrossTermForScale,
+ NonStationaryAltitudinalLocationQuadraticScaleQuadratic):
+ pass
+
+
+
+
+
diff --git a/extreme_fit/model/margin_model/polynomial_margin_model/gev_altitudinal_models_with_scale.py b/extreme_fit/model/margin_model/polynomial_margin_model/gev_altitudinal_models_with_scale.py
new file mode 100644
index 00000000..1bf7701f
--- /dev/null
+++ b/extreme_fit/model/margin_model/polynomial_margin_model/gev_altitudinal_models_with_scale.py
@@ -0,0 +1,57 @@
+from extreme_fit.distribution.gev.gev_params import GevParams
+from extreme_fit.model.margin_model.polynomial_margin_model.gev_altitudinal_models import AbstractAltitudinalModel, \
+ AbstractAddCrossTermForLocation
+from extreme_fit.model.margin_model.polynomial_margin_model.polynomial_margin_model import PolynomialMarginModel
+from extreme_fit.model.margin_model.polynomial_margin_model.spatio_temporal_polynomial_model import \
+ AbstractSpatioTemporalPolynomialModel
+
+
+
+
+
+
+class NonStationaryAltitudinalScaleQuadratic(AbstractAltitudinalModel):
+
+ @property
+ def param_name_to_list_dim_and_degree(self):
+ return {
+ GevParams.LOC: [(self.coordinates.idx_x_coordinates, 1)],
+ GevParams.SCALE: [(self.coordinates.idx_x_coordinates, 1), (self.coordinates.idx_temporal_coordinates, 2)],
+ }
+
+class NonStationaryAltitudinalLocationLinearScaleQuadratic(AbstractAltitudinalModel):
+
+ @property
+ def param_name_to_list_dim_and_degree(self):
+ return {
+ GevParams.LOC: [(self.coordinates.idx_x_coordinates, 1), (self.coordinates.idx_temporal_coordinates, 1)],
+ GevParams.SCALE: [(self.coordinates.idx_x_coordinates, 1), (self.coordinates.idx_temporal_coordinates, 2)],
+ }
+
+class NonStationaryAltitudinalLocationQuadraticScaleQuadratic(AbstractAltitudinalModel):
+
+ @property
+ def param_name_to_list_dim_and_degree(self):
+ return {
+ GevParams.LOC: [(self.coordinates.idx_x_coordinates, 1), (self.coordinates.idx_temporal_coordinates, 2)],
+ GevParams.SCALE: [(self.coordinates.idx_x_coordinates, 1), (self.coordinates.idx_temporal_coordinates, 2)],
+ }
+
+
+# Add cross terms
+
+
+
+class NonStationaryAltitudinalScaleQuadraticCrossTermForLocation(AbstractAddCrossTermForLocation, NonStationaryAltitudinalScaleQuadratic):
+ pass
+
+
+class NonStationaryAltitudinalLocationLinearScaleQuadraticCrossTermForLocation(AbstractAddCrossTermForLocation,
+ NonStationaryAltitudinalLocationLinearScaleQuadratic):
+ pass
+
+
+class NonStationaryAltitudinalLocationQuadraticScaleQuadraticCrossTermForLocation(AbstractAddCrossTermForLocation,
+ NonStationaryAltitudinalLocationQuadraticScaleQuadratic):
+ pass
+
diff --git a/extreme_fit/model/margin_model/polynomial_margin_model/utils.py b/extreme_fit/model/margin_model/polynomial_margin_model/utils.py
index 0ff5aaad..1a002b60 100644
--- a/extreme_fit/model/margin_model/polynomial_margin_model/utils.py
+++ b/extreme_fit/model/margin_model/polynomial_margin_model/utils.py
@@ -6,6 +6,18 @@ from extreme_fit.model.margin_model.polynomial_margin_model.gev_altitudinal_mode
NonStationaryAltitudinalLocationLinearScaleLinearCrossTermForLocation, \
NonStationaryAltitudinalLocationQuadraticScaleLinearCrossTermForLocation, NonStationaryAltitudinalScaleLinear, \
NonStationaryAltitudinalScaleLinearCrossTermForLocation
+from extreme_fit.model.margin_model.polynomial_margin_model.gev_altitudinal_models_cross_term_in_scale import \
+ NonStationaryAltitudinalLocationLinearScaleQuadraticCrossTermForScale, \
+ NonStationaryAltitudinalLocationQuadraticScaleQuadraticCrossTermForScale, \
+ NonStationaryAltitudinalScaleLinearCrossTermForScale, \
+ NonStationaryAltitudinalLocationLinearScaleLinearCrossTermForScale, \
+ NonStationaryAltitudinalLocationQuadraticScaleLinearCrossTermForScale, \
+ NonStationaryAltitudinalScaleQuadraticCrossTermForScale
+from extreme_fit.model.margin_model.polynomial_margin_model.gev_altitudinal_models_with_scale import \
+ NonStationaryAltitudinalScaleQuadraticCrossTermForLocation, NonStationaryAltitudinalScaleQuadratic, \
+ NonStationaryAltitudinalLocationLinearScaleQuadratic, NonStationaryAltitudinalLocationQuadraticScaleQuadratic, \
+ NonStationaryAltitudinalLocationLinearScaleQuadraticCrossTermForLocation, \
+ NonStationaryAltitudinalLocationQuadraticScaleQuadraticCrossTermForLocation
from extreme_fit.model.margin_model.polynomial_margin_model.gumbel_altitudinal_models import \
StationaryGumbelAltitudinal, NonStationaryGumbelAltitudinalScaleLinear, \
NonStationaryGumbelAltitudinalLocationLinear, NonStationaryGumbelAltitudinalLocationQuadratic, \
@@ -22,6 +34,7 @@ from extreme_fit.model.margin_model.polynomial_margin_model.spatio_temporal_poly
NonStationaryLocationSpatioTemporalLinearityModelAssertError2, \
NonStationaryLocationSpatioTemporalLinearityModelAssertError3, NonStationaryLocationSpatioTemporalLinearityModel6
+
ALTITUDINAL_GEV_MODELS = [
StationaryAltitudinal,
NonStationaryAltitudinalScaleLinear,
@@ -36,8 +49,26 @@ ALTITUDINAL_GEV_MODELS = [
NonStationaryAltitudinalLocationQuadraticCrossTermForLocation,
NonStationaryAltitudinalLocationLinearScaleLinearCrossTermForLocation,
NonStationaryAltitudinalLocationQuadraticScaleLinearCrossTermForLocation,
+
+ # Quadratic in the scale
+ NonStationaryAltitudinalScaleQuadratic,
+ NonStationaryAltitudinalLocationLinearScaleQuadratic,
+ NonStationaryAltitudinalLocationQuadraticScaleQuadratic,
+ NonStationaryAltitudinalScaleQuadraticCrossTermForLocation,
+ NonStationaryAltitudinalLocationLinearScaleQuadraticCrossTermForLocation,
+ NonStationaryAltitudinalLocationQuadraticScaleQuadraticCrossTermForLocation,
+
+ # Cross term for the scale
+ NonStationaryAltitudinalScaleLinearCrossTermForScale,
+ NonStationaryAltitudinalLocationLinearScaleLinearCrossTermForScale,
+ NonStationaryAltitudinalLocationQuadraticScaleLinearCrossTermForScale,
+ NonStationaryAltitudinalScaleQuadraticCrossTermForScale,
+ NonStationaryAltitudinalLocationLinearScaleQuadraticCrossTermForScale,
+ NonStationaryAltitudinalLocationQuadraticScaleQuadraticCrossTermForScale,
+
][:]
+
ALTITUDINAL_GUMBEL_MODELS = [
StationaryGumbelAltitudinal,
NonStationaryGumbelAltitudinalScaleLinear,
diff --git a/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py b/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
index 5b013a5f..56b470b2 100644
--- a/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
+++ b/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
@@ -2,6 +2,7 @@ import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
+from extreme_fit.model.margin_model.polynomial_margin_model.utils import ALTITUDINAL_GEV_MODELS
from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.plot_total_aic import plot_total_aic, \
plot_individual_aic
@@ -16,8 +17,6 @@ from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranS
SafranSnowfall5Days, SafranSnowfall7Days, SafranPrecipitation1Day, SafranPrecipitation3Days, \
SafranPrecipitation5Days, SafranPrecipitation7Days
from extreme_data.meteo_france_data.scm_models_data.utils import Season
-from extreme_fit.model.margin_model.polynomial_margin_model.utils import ALTITUDINAL_GEV_MODELS, \
- ALTITUDINAL_GUMBEL_MODELS
from projects.altitude_spatial_model.altitudes_fit.altitudes_studies import AltitudesStudies
from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitudes_studies_visualizer_for_non_stationary_models import \
AltitudesStudiesVisualizerForNonStationaryModels
@@ -34,24 +33,23 @@ def plot_moments(studies, massif_names=None):
def plot_altitudinal_fit(studies, massif_names=None):
- model_classes = ALTITUDINAL_GEV_MODELS + ALTITUDINAL_GUMBEL_MODELS
+ model_classes = ALTITUDINAL_GEV_MODELS
visualizer = AltitudesStudiesVisualizerForNonStationaryModels(studies=studies,
model_classes=model_classes,
massif_names=massif_names,
show=False)
# Plot the results for the model that minimizes the individual aic
- plot_individual_aic(visualizer)
+ # plot_individual_aic(visualizer)
# Plot the results for the model that minimizes the total aic
plot_total_aic(model_classes, visualizer)
def main():
- altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][4:7]
- # altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][:]
+ altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][4:6]
# todo: l ecart pour les saisons de l automne ou de sprint
# vient probablement de certains zéros pas pris en compte pour le fit,
# mais pris en compte pour le calcul de mon aic
- # altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900][:]
+ altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900][:]
study_classes = [SafranSnowfall1Day, SafranSnowfall3Days, SafranSnowfall5Days, SafranSnowfall7Days][:2]
study_classes = [SafranPrecipitation1Day, SafranPrecipitation3Days, SafranPrecipitation5Days,
SafranPrecipitation7Days][:]
@@ -59,6 +57,7 @@ def main():
SafranSnowfall3Days, SafranPrecipitation3Days][:1]
# seasons = [Season.automn, Season.winter, Season.spring][::-1]
seasons = [Season.winter]
+ # seasons = [Season.winter_extended]
massif_names = None
# massif_names = ['Aravis']
diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
index 859e661f..351f9c7c 100644
--- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
+++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
@@ -26,7 +26,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
show=False,
massif_names=None,
fit_method=MarginFitMethod.extremes_fevd_mle,
- display_only_model_that_pass_anderson_test=True):
+ display_only_model_that_pass_anderson_test=False):
super().__init__(studies.study, show=show, save_to_file=not show)
self.studies = studies
self.non_stationary_models = model_classes
@@ -127,8 +127,9 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
self.plot_abstract_fast(massif_name_to_best_coef,
label='{}/Coef/{} plot for {} {}'.format(OneFoldFit.folder_for_plots,
coef_name,
- SCM_STUDY_CLASS_TO_ABBREVIATION[type(self.study)],
- self.study.variable_unit),
+ SCM_STUDY_CLASS_TO_ABBREVIATION[
+ type(self.study)],
+ self.study.variable_unit),
add_x_label=False, graduation=graduation, massif_name_to_text=self.massif_name_to_name,
negative_and_positive_values=negative_and_positive_values)
@@ -138,5 +139,30 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
self.study.variable_unit),
add_x_label=False, graduation=0.1, massif_name_to_text=self.massif_name_to_name)
- def plot_altitude_change_of_sign(self):
+ def plot_altitude_for_the_peak(self):
pass
+
+ def plot_year_for_the_peak(self):
+ t_list = 1900 + np.arange(200)
+ for massif_name, one_fold_fit in self.massif_name_to_one_fold_fit.items():
+ ax = plt.gca()
+ # One plot for each altitude
+ print(one_fold_fit.dataset.df_dataset.head())
+ altitudes = self.studies.altitudes
+ print(altitudes)
+ for altitude in altitudes:
+
+ mean_list = []
+ for t in t_list:
+ coordinate = np.array([altitude, t])
+ gev_params = one_fold_fit.best_function_from_fit.get_params(coordinate, is_transformed=False)
+ mean_list.append(gev_params.mean)
+ label = '{} m'.format(altitude)
+ ax.plot(t_list, mean_list, label=label)
+ ax.legend()
+ ax.set_xlabel('Year')
+ ax.set_xlabel('Mean annual maxima of ')
+ self.plot_name = 'Peak year for {}'.format(massif_name)
+ self.show_or_save_to_file(no_title=True)
+ plt.close()
+
diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py
index 1a64a7a3..62d514de 100644
--- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py
+++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py
@@ -163,7 +163,9 @@ class OneFoldFit(object):
quantiles = self.compute_empirical_quantiles(estimator=self.best_estimator)
goodness_of_fit_anderson_test = goodness_of_fit_anderson(quantiles, self.SIGNIFICANCE_LEVEL)
if not goodness_of_fit_anderson_test:
- print('{} with {} does not pass the anderson test'.format(self.massif_name, self.folder_for_plots))
+ print('{} with {} does not pass the anderson test for model {}'.format(self.massif_name,
+ self.folder_for_plots,
+ type(self.best_margin_model)))
return goodness_of_fit_anderson_test
def compute_empirical_quantiles(self, estimator):
diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/plot_total_aic.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/plot_total_aic.py
index 0c9bba51..9c50cc3e 100644
--- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/plot_total_aic.py
+++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/plot_total_aic.py
@@ -10,11 +10,16 @@ from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.one_fold_fi
from projects.exceeding_snow_loads.utils import dpi_paper1_figure
-def plot_individual_aic(visualizer):
- OneFoldFit.best_estimator_minimizes_total_aic = False
+def plots(visualizer):
visualizer.plot_moments()
- # visualizer.plot_best_coef_maps()
+ visualizer.plot_best_coef_maps()
visualizer.plot_shape_map()
+ visualizer.plot_year_for_the_peak()
+
+
+def plot_individual_aic(visualizer):
+ OneFoldFit.best_estimator_minimizes_total_aic = False
+ plots(visualizer)
def plot_total_aic(model_classes, visualizer):
@@ -45,9 +50,7 @@ def plot_total_aic(model_classes, visualizer):
# Plot the ranking of the model based on their total aic
plot_total_aic_repartition(visualizer, sorted_labels, sorted_scores)
# Plot the results for the model that minimizes the mean aic
- visualizer.plot_moments()
- visualizer.plot_shape_map()
- visualizer.plot_best_coef_maps()
+ plots(visualizer)
def plot_total_aic_repartition(visualizer, labels, scores):
--
GitLab