Commit c7c042d4 authored by Le Roux Erwan's avatar Le Roux Erwan
Browse files

[contrasting] display individual model fitted with all data (300-4800)....

[contrasting] display individual model fitted with all data (300-4800). Display only the best non-stationary model if it passes the Anderson fit test.
parent d1822a0b
No related merge requests found
Showing with 211 additions and 26 deletions
+211 -26
......@@ -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
......
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
......@@ -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,
......
......@@ -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,
......
......@@ -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('_', ''))
......
......@@ -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):
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment