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 7f1381c6ec4ed97920ad8e549623b5f15b542c0c..547f39f1723004072ff4205d600e40eff9fd7e82 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 @@ -37,6 +37,9 @@ class AbstractAltitudinalModel(AbstractSpatioTemporalPolynomialModel): name = self.DISTRIBUTION_STR 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) + shape_str_number = self.dim_to_str_number(GevParams.SHAPE, self.coordinates.idx_temporal_coordinates) + if shape_str_number != '0': + name += shape_str_number if isinstance(self, AbstractAddCrossTermForLocation): name += 'x' if isinstance(self, AbstractAddCrossTermForScale): @@ -126,6 +129,14 @@ class NonStationaryAltitudinalLocationQuadraticScaleLinear(AbstractAltitudinalMo class AbstractAddCrossTermForLocation(AbstractAltitudinalModel): + # @property + # def param_name_to_list_dim_and_degree_for_margin_function(self): + # d = self.param_name_to_list_dim_and_degree + # assert 1 <= len(d[GevParams.LOC]) <= 2 + # assert self.coordinates.idx_x_coordinates == d[GevParams.LOC][0][0] + # d[GevParams.LOC].insert(1, ((self.coordinates.idx_x_coordinates, self.coordinates.idx_temporal_coordinates), 1)) + # return d + @property def param_name_to_list_dim_and_degree_for_margin_function(self): d = self.param_name_to_list_dim_and_degree diff --git a/extreme_fit/model/margin_model/polynomial_margin_model/gev_altitudinal_models_with_scale_2.py b/extreme_fit/model/margin_model/polynomial_margin_model/gev_altitudinal_models_with_scale_2.py new file mode 100644 index 0000000000000000000000000000000000000000..8139ff23fd1285f156f38b4a3a3b5473a4898dd9 --- /dev/null +++ b/extreme_fit/model/margin_model/polynomial_margin_model/gev_altitudinal_models_with_scale_2.py @@ -0,0 +1,129 @@ +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 NonStationaryAltitudinalScaleLinear(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, 1)], + } + +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 NonStationaryAltitudinalScaleLinearLocationLinear(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, 1)], + } + +class NonStationaryAltitudinalScaleQuadraticLocationLinear(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 NonStationaryAltitudinalLocationQuadraticScaleLinear(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, 1)], + } + +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)], + } + +class NonStationaryAltitudinalLocationCubicScaleLinear(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, 3)], + GevParams.SCALE: [(self.coordinates.idx_x_coordinates, 1), (self.coordinates.idx_temporal_coordinates, 1)], + } + +class NonStationaryAltitudinalLocationCubicScaleQuadratic(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, 3)], + GevParams.SCALE: [(self.coordinates.idx_x_coordinates, 1), (self.coordinates.idx_temporal_coordinates, 2)], + } + +# Add cross terms + + + +class NonStationaryAltitudinalScaleLinearCrossTermForLocation(AbstractAddCrossTermForLocation, NonStationaryAltitudinalScaleLinear): + pass + + +class NonStationaryAltitudinalScaleLinearLocationLinearCrossTermForLocation(AbstractAddCrossTermForLocation, + NonStationaryAltitudinalScaleLinearLocationLinear): + pass + + +class NonStationaryAltitudinalLocationQuadraticScaleLinearCrossTermForLocation(AbstractAddCrossTermForLocation, + NonStationaryAltitudinalLocationQuadraticScaleLinear): + pass +class NonStationaryAltitudinalLocationCubicScaleLinearCrossTermForLocation(AbstractAddCrossTermForLocation, + NonStationaryAltitudinalLocationCubicScaleLinear): + pass + + + + +class NonStationaryAltitudinalScaleQuadraticCrossTermForLocation(AbstractAddCrossTermForLocation, NonStationaryAltitudinalScaleQuadratic): + pass + + +class NonStationaryAltitudinalScaleQuadraticLocationLinearCrossTermForLocation(AbstractAddCrossTermForLocation, + NonStationaryAltitudinalScaleQuadraticLocationLinear): + pass + + + +class NonStationaryAltitudinalLocationQuadraticScaleQuadraticCrossTermForLocation(AbstractAddCrossTermForLocation, + NonStationaryAltitudinalLocationQuadraticScaleQuadratic): + pass + + +class NonStationaryAltitudinalLocationCubicScaleQuadraticCrossTermForLocation(AbstractAddCrossTermForLocation, + NonStationaryAltitudinalLocationCubicScaleQuadratic): + 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 e0288df748c957bd61e564c2205420247d0790ad..7a61879eac82d1c93f243d4265164d24faca06f8 100644 --- a/extreme_fit/model/margin_model/polynomial_margin_model/utils.py +++ b/extreme_fit/model/margin_model/polynomial_margin_model/utils.py @@ -28,6 +28,13 @@ from extreme_fit.model.margin_model.polynomial_margin_model.gev_altitudinal_mode NonStationaryAltitudinalLocationLinearScaleQuadratic, NonStationaryAltitudinalLocationQuadraticScaleQuadratic, \ NonStationaryAltitudinalLocationLinearScaleQuadraticCrossTermForLocation, \ NonStationaryAltitudinalLocationQuadraticScaleQuadraticCrossTermForLocation +from extreme_fit.model.margin_model.polynomial_margin_model.gev_altitudinal_models_with_scale_2 import \ + NonStationaryAltitudinalScaleLinearCrossTermForLocation, \ + NonStationaryAltitudinalScaleLinearLocationLinearCrossTermForLocation, \ + NonStationaryAltitudinalLocationQuadraticScaleLinearCrossTermForLocation, \ + NonStationaryAltitudinalScaleQuadraticLocationLinearCrossTermForLocation, \ + NonStationaryAltitudinalLocationCubicScaleQuadraticCrossTermForLocation, \ + NonStationaryAltitudinalLocationCubicScaleLinearCrossTermForLocation from extreme_fit.model.margin_model.polynomial_margin_model.gumbel_altitudinal_models import \ StationaryGumbelAltitudinal, NonStationaryGumbelAltitudinalScaleLinear, \ NonStationaryGumbelAltitudinalLocationLinear, NonStationaryGumbelAltitudinalLocationQuadratic, \ @@ -44,37 +51,61 @@ from extreme_fit.model.margin_model.polynomial_margin_model.spatio_temporal_poly NonStationaryLocationSpatioTemporalLinearityModelAssertError2, \ NonStationaryLocationSpatioTemporalLinearityModelAssertError3, NonStationaryLocationSpatioTemporalLinearityModel6 - ALTITUDINAL_GEV_MODELS_LOCATION = [ StationaryAltitudinal, - NonStationaryAltitudinalLocationLinear, + # NonStationaryAltitudinalLocationLinear, + # NonStationaryAltitudinalLocationQuadratic, + # NonStationaryAltitudinalLocationCubic, + # NonStationaryAltitudinalLocationOrder4, + # # First order cross term + # NonStationaryCrossTermForLocation, + + # NonStationaryAltitudinalLocationLinearCrossTermForLocation, + NonStationaryAltitudinalLocationQuadraticCrossTermForLocation, + NonStationaryAltitudinalLocationCubicCrossTermForLocation, + # NonStationaryAltitudinalLocationOrder4CrossTermForLocation, + # NonStationaryAltitudinalScaleLinearCrossTermForLocation, + # NonStationaryAltitudinalScaleLinearLocationLinearCrossTermForLocation, + + # NonStationaryAltitudinalScaleQuadraticCrossTermForLocation, + # NonStationaryAltitudinalScaleQuadraticLocationLinearCrossTermForLocation, + + # NonStationaryAltitudinalLocationQuadraticScaleLinearCrossTermForLocation, + # NonStationaryAltitudinalLocationCubicScaleLinearCrossTermForLocation, + # NonStationaryAltitudinalLocationQuadraticScaleQuadraticCrossTermForLocation, + # NonStationaryAltitudinalLocationCubicScaleQuadraticCrossTermForLocation, + +] +ALTITUDINAL_GEV_MODELS_LOCATION_QUADRATIC_MINIMUM = [ + StationaryAltitudinal, + # NonStationaryAltitudinalLocationLinear, NonStationaryAltitudinalLocationQuadratic, NonStationaryAltitudinalLocationCubic, NonStationaryAltitudinalLocationOrder4, - # First order cross term - NonStationaryCrossTermForLocation, - NonStationaryAltitudinalLocationLinearCrossTermForLocation, + # # First order cross term + # NonStationaryCrossTermForLocation, + # NonStationaryAltitudinalLocationLinearCrossTermForLocation, NonStationaryAltitudinalLocationQuadraticCrossTermForLocation, NonStationaryAltitudinalLocationCubicCrossTermForLocation, NonStationaryAltitudinalLocationOrder4CrossTermForLocation, ] + ALTITUDINAL_GEV_MODELS_LOCATION_ONLY_SCALE_ALTITUDES = [ -StationaryAltitudinalOnlyScale, -NonStationaryAltitudinalOnlyScaleLocationLinear, -NonStationaryAltitudinalOnlyScaleLocationQuadratic, -NonStationaryAltitudinalOnlyScaleLocationCubic, -NonStationaryAltitudinalOnlyScaleLocationOrder4, -# Cross terms -NonStationaryOnlyScaleCrossTermForLocation, -NonStationaryAltitudinalOnlyScaleLocationLinearCrossTermForLocation, -NonStationaryAltitudinalOnlyScaleLocationQuadraticCrossTermForLocation, -NonStationaryAltitudinalOnlyScaleLocationCubicCrossTermForLocation, -NonStationaryAltitudinalOnlyScaleLocationOrder4CrossTermForLocation, + StationaryAltitudinalOnlyScale, + NonStationaryAltitudinalOnlyScaleLocationLinear, + NonStationaryAltitudinalOnlyScaleLocationQuadratic, + NonStationaryAltitudinalOnlyScaleLocationCubic, + NonStationaryAltitudinalOnlyScaleLocationOrder4, + # Cross terms + NonStationaryOnlyScaleCrossTermForLocation, + NonStationaryAltitudinalOnlyScaleLocationLinearCrossTermForLocation, + NonStationaryAltitudinalOnlyScaleLocationQuadraticCrossTermForLocation, + NonStationaryAltitudinalOnlyScaleLocationCubicCrossTermForLocation, + NonStationaryAltitudinalOnlyScaleLocationOrder4CrossTermForLocation, ] - ALTITUDINAL_GEV_MODELS = [ StationaryAltitudinal, NonStationaryAltitudinalScaleLinear, @@ -108,7 +139,6 @@ ALTITUDINAL_GEV_MODELS = [ ][:] - 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 ab7989561979b760521f8cf46a6e5adf59be609a..5508f06c0f8a5fc6a0ef68a20a1ce1e418522835 100644 --- a/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py +++ b/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py @@ -3,7 +3,7 @@ import matplotlib.pyplot as plt import numpy as np from extreme_fit.model.margin_model.polynomial_margin_model.utils import ALTITUDINAL_GEV_MODELS, \ - ALTITUDINAL_GEV_MODELS_LOCATION, ALTITUDINAL_GEV_MODELS_LOCATION_ONLY_SCALE_ALTITUDES + ALTITUDINAL_GEV_MODELS_LOCATION_QUADRATIC_MINIMUM, ALTITUDINAL_GEV_MODELS_LOCATION_ONLY_SCALE_ALTITUDES, ALTITUDINAL_GEV_MODELS_LOCATION from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.plot_total_aic import plot_total_aic, \ plot_individual_aic @@ -35,16 +35,18 @@ def plot_moments(studies, massif_names=None): def plot_altitudinal_fit(studies, massif_names=None): # model_classes = ALTITUDINAL_GEV_MODELS_LOCATION_ONLY_SCALE_ALTITUDES + # model_classes = ALTITUDINAL_GEV_MODELS_LOCATION_QUADRATIC_MINIMUM model_classes = ALTITUDINAL_GEV_MODELS_LOCATION + # model_classes = ALTITUDINAL_GEV_MODELS_LOCATION_CUBIC_MINIMUM # model_classes = ALTITUDINAL_GEV_MODELS_QUADRATIC 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) + # plot_total_aic(model_classes, visualizer) def main(): @@ -54,6 +56,9 @@ def main(): # mais pris en compte pour le calcul de mon aic # altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][:] altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900, 4200, 4500, 4800][:] + # altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900][:] + # altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][:] + # altitudes = [600, 900, 1200, 1500, 1800, 2100][:] # altitudes = [1500, 1800][:] study_classes = [SafranSnowfall1Day, SafranSnowfall3Days, SafranSnowfall5Days, SafranSnowfall7Days][:2] study_classes = [SafranPrecipitation1Day, SafranPrecipitation3Days, SafranPrecipitation5Days, 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 70215b28a55252ffcd6b96a79dbab428c94ca268..db0e3634ea867332ec751352701ceaaa75f2bddc 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=False): + display_only_model_that_pass_anderson_test=True): super().__init__(studies.study, show=show, save_to_file=not show) self.studies = studies self.non_stationary_models = model_classes @@ -100,7 +100,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer): dim_to_coordinate_name = dict(zip([0, 1], coordinate_names)) for dim in [0, 1, (0, 1)]: coordinate_name = LinearCoef.coefficient_name(dim, dim_to_coordinate_name) - for degree in range(3): + for degree in range(4): coef_name = ' '.join([param_name + coordinate_name + str(degree)]) massif_name_to_best_coef = {} for massif_name, one_fold_fit in self.massif_name_to_one_fold_fit.items(): @@ -109,7 +109,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer): massif_name_to_best_coef[massif_name] = best_coef if len(massif_name_to_best_coef) > 0: - for evaluate_coordinate in [False, True]: + for evaluate_coordinate in [False, True][:1]: if evaluate_coordinate: coef_name += 'evaluated at coordinates' for massif_name in massif_name_to_best_coef.keys(): @@ -144,10 +144,11 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer): def plot_year_for_the_peak(self): t_list = 1959 + np.arange(141) + # t_list = 1800 + np.arange(400) for massif_name, one_fold_fit in self.massif_name_to_one_fold_fit.items(): ax = plt.gca() # One plot for each altitude - altitudes = np.arange(500, 3000, 500) + altitudes = np.arange(500, min(3000, max(self.studies.altitudes)), 500) for altitude in altitudes: i = 0 while self.studies.altitudes[i] < altitude: @@ -163,6 +164,10 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer): label = '{} m'.format(altitude) ax.plot(t_list, mean_list, label=label) ax.legend() + # Modify the limits of the y axis + lim_down, lim_up = ax.get_ylim() + ax_lim = (0, lim_up) + ax.set_ylim(ax_lim) ax.set_xlabel('Year') ax.set_ylabel('Mean annual maxima of {}'.format(SCM_STUDY_CLASS_TO_ABBREVIATION[type(self.study)])) plot_name = '{}/{}/Peak year for {}'.format(OneFoldFit.folder_for_plots, 'Peak year', massif_name.replace('_', '')) 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 a5d3c9f97f09f7731d3856cd78ecbf67b12f2e4d..2687b7ce0bf27679bf6fb9d6fca252b2ef35bbe2 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 @@ -97,6 +97,10 @@ class OneFoldFit(object): sorted_estimators = sorted([estimator for estimator in estimators], key=lambda e: e.aic()) return sorted_estimators + @cached_property + def sorted_estimators_without_stationary(self): + return [e for e in self.sorted_estimators if not isinstance(e.margin_model, StationaryAltitudinal)] + @property def model_class_to_estimator_with_finite_aic(self): return {type(estimator.margin_model): estimator for estimator in self.sorted_estimators} @@ -106,7 +110,8 @@ class OneFoldFit(object): if self.best_estimator_minimizes_total_aic and self.best_estimator_class_for_total_aic is not None: return self.model_class_to_estimator[self.best_estimator_class_for_total_aic] else: - return self.sorted_estimators[0] + best_estimator = self.sorted_estimators_without_stationary[0] + return best_estimator @property def best_margin_model(self):