Commit 19a423d2 authored by Le Roux Erwan's avatar Le Roux Erwan
Browse files

[exceeding] update style of the qqplot. clean...

[exceeding] update style of the qqplot. clean study_visualizer_for_non_stationary_trends.py, add a mode with only selected models that does not crash.
parent e4b82417
No related merge requests found
Showing with 22 additions and 62 deletions
+22 -62
...@@ -372,11 +372,11 @@ class AbstractGevTrendTest(object): ...@@ -372,11 +372,11 @@ class AbstractGevTrendTest(object):
size = 15 size = 15
standard_gumbel_quantiles = self.get_standard_gumbel_quantiles() standard_gumbel_quantiles = self.get_standard_gumbel_quantiles()
unconstrained_empirical_quantiles = self.compute_empirical_quantiles(self.unconstrained_estimator) unconstrained_empirical_quantiles = self.compute_empirical_quantiles(self.unconstrained_estimator)
constrained_empirical_quantiles = self.compute_empirical_quantiles(self.constrained_estimator) # constrained_empirical_quantiles = self.compute_empirical_quantiles(self.constrained_estimator)
all_quantiles = standard_gumbel_quantiles + unconstrained_empirical_quantiles + constrained_empirical_quantiles all_quantiles = standard_gumbel_quantiles + unconstrained_empirical_quantiles
epsilon = 0.5 epsilon = 0.1
ax_lim = [min(all_quantiles) - epsilon, max(all_quantiles) + epsilon] ax_lim = [min(all_quantiles) - epsilon, max(all_quantiles) + epsilon]
ax.plot(standard_gumbel_quantiles, standard_gumbel_quantiles, color='k')
# ax.plot(standard_gumbel_quantiles, constrained_empirical_quantiles, 'x', # ax.plot(standard_gumbel_quantiles, constrained_empirical_quantiles, 'x',
# label='Stationary Gumbel model $\mathcal{M}_0$') # label='Stationary Gumbel model $\mathcal{M}_0$')
...@@ -390,10 +390,12 @@ class AbstractGevTrendTest(object): ...@@ -390,10 +390,12 @@ class AbstractGevTrendTest(object):
ax.legend(loc='lower right', prop={'size': 10}) ax.legend(loc='lower right', prop={'size': 10})
ax.set_xlim(ax_lim) ax.set_xlim(ax_lim)
ax.set_ylim(ax_lim) ax.set_ylim(ax_lim)
ax.plot(ax_lim, ax_lim, color='k')
ticks = [i for i in range(ceil(ax_lim[0]), floor(ax_lim[1]) + 1)] ticks = [i for i in range(ceil(ax_lim[0]), floor(ax_lim[1]) + 1)]
ax.set_xticks(ticks) ax.set_xticks(ticks)
ax.set_yticks(ticks) ax.set_yticks(ticks)
ax.grid() # ax.grid()
ax.tick_params(labelsize=size) ax.tick_params(labelsize=size)
def get_standard_gumbel_quantiles(self): def get_standard_gumbel_quantiles(self):
......
...@@ -144,59 +144,21 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer): ...@@ -144,59 +144,21 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
d = {m: v for m, v in d.items() if self.massif_name_to_psnow[m] >= 0.9} d = {m: v for m, v in d.items() if self.massif_name_to_psnow[m] >= 0.9}
return d return d
@property
def massif_name_to_trend_test_that_minimized_aic(self) -> Dict[str, AbstractGevTrendTest]:
return self.massif_name_to_trend_test_tuple[0]
@property
def massif_name_to_stationary_trend_test_that_minimized_aic(self) -> Dict[str, AbstractGevTrendTest]:
return self.massif_name_to_trend_test_tuple[1]
@property
def massif_name_to_gumbel_trend_test_that_minimized_aic(self) -> Dict[str, AbstractGevTrendTest]:
return self.massif_name_to_trend_test_tuple[2]
@cached_property @cached_property
def massif_name_to_trend_test_tuple(self) -> Tuple[ def massif_name_to_trend_test_that_minimized_aic(self) -> Dict[str, AbstractGevTrendTest]:
Dict[str, AbstractGevTrendTest], Dict[str, AbstractGevTrendTest], Dict[str, AbstractGevTrendTest]]:
print('cached', self.altitude, id(self))
massif_name_to_trend_test_that_minimized_aic = {} massif_name_to_trend_test_that_minimized_aic = {}
massif_name_to_stationary_trend_test_that_minimized_aic = {}
massif_name_to_gumbel_trend_test_that_minimized_aic = {}
for massif_name in self.massif_name_to_years_and_maxima_for_model_fitting.keys(): for massif_name in self.massif_name_to_years_and_maxima_for_model_fitting.keys():
# Compute sorted trend test # Compute sorted trend test
sorted_trend_test = self.get_sorted_trend_test(massif_name) sorted_trend_test = self.get_sorted_trend_test(massif_name)
# Extract the stationary or non-stationary model that minimized AIC # Extract the stationary or non-stationary model that minimized AIC
trend_test_that_minimized_aic = sorted_trend_test[0] trend_test_that_minimized_aic = sorted_trend_test[0]
if self.select_only_model_that_pass_anderson_test and \ if (not self.select_only_model_that_pass_anderson_test) or \
(not trend_test_that_minimized_aic.goodness_of_fit_anderson_test): trend_test_that_minimized_aic.goodness_of_fit_anderson_test:
print('not anderson')
else:
print('ok')
massif_name_to_trend_test_that_minimized_aic[massif_name] = trend_test_that_minimized_aic massif_name_to_trend_test_that_minimized_aic[massif_name] = trend_test_that_minimized_aic
# Extract the stationary model that minimized AIC
stationary_trend_tests_that_minimized_aic = [t for t in sorted_trend_test if type(t) in
[GumbelVersusGumbel, GevStationaryVersusGumbel]]
if len(stationary_trend_tests_that_minimized_aic) == 0:
stationary_trend_test_that_minimized_aic = None
else:
stationary_trend_test_that_minimized_aic = stationary_trend_tests_that_minimized_aic[0]
massif_name_to_stationary_trend_test_that_minimized_aic[
massif_name] = stationary_trend_test_that_minimized_aic
# Extract the Gumbel model that minimized AIC
gumbel_trend_tests = [t for t in sorted_trend_test if
type(t) in [GumbelVersusGumbel, GumbelLocationTrendTest, GumbelScaleTrendTest,
GumbelLocationAndScaleTrendTest]]
if len(gumbel_trend_tests) > 0:
gumbel_trend_test_that_minimized_aic = gumbel_trend_tests[0]
else: else:
gumbel_trend_test_that_minimized_aic = None print('anderson test was not passed')
massif_name_to_gumbel_trend_test_that_minimized_aic[massif_name] = gumbel_trend_test_that_minimized_aic
return massif_name_to_trend_test_that_minimized_aic, massif_name_to_stationary_trend_test_that_minimized_aic, massif_name_to_gumbel_trend_test_that_minimized_aic return massif_name_to_trend_test_that_minimized_aic
def get_sorted_trend_test(self, massif_name): def get_sorted_trend_test(self, massif_name):
x, y = self.massif_name_to_years_and_maxima_for_model_fitting[massif_name] x, y = self.massif_name_to_years_and_maxima_for_model_fitting[massif_name]
...@@ -360,12 +322,6 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer): ...@@ -360,12 +322,6 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
def massif_name_and_model_subset_to_model_class(self, massif_name, model_subset_for_uncertainty): def massif_name_and_model_subset_to_model_class(self, massif_name, model_subset_for_uncertainty):
if model_subset_for_uncertainty is ModelSubsetForUncertainty.stationary_gumbel: if model_subset_for_uncertainty is ModelSubsetForUncertainty.stationary_gumbel:
return GumbelTemporalModel return GumbelTemporalModel
if model_subset_for_uncertainty is ModelSubsetForUncertainty.stationary_gev:
return StationaryTemporalModel
elif model_subset_for_uncertainty is ModelSubsetForUncertainty.stationary_gumbel_and_gev:
return self.massif_name_to_stationary_trend_test_that_minimized_aic[massif_name].unconstrained_model_class
elif model_subset_for_uncertainty is ModelSubsetForUncertainty.non_stationary_gumbel:
return self.massif_name_to_gumbel_trend_test_that_minimized_aic[massif_name].unconstrained_model_class
elif model_subset_for_uncertainty is ModelSubsetForUncertainty.non_stationary_gumbel_and_gev: elif model_subset_for_uncertainty is ModelSubsetForUncertainty.non_stationary_gumbel_and_gev:
return self.massif_name_to_trend_test_that_minimized_aic[massif_name].unconstrained_model_class return self.massif_name_to_trend_test_that_minimized_aic[massif_name].unconstrained_model_class
else: else:
...@@ -378,6 +334,10 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer): ...@@ -378,6 +334,10 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
# Compute for the uncertainty massif names # Compute for the uncertainty massif names
massifs_names = set(self.massif_name_to_years_and_maxima_for_model_fitting.keys()). \ massifs_names = set(self.massif_name_to_years_and_maxima_for_model_fitting.keys()). \
intersection(self.uncertainty_massif_names) intersection(self.uncertainty_massif_names)
# Update the name of the massif (because some massifs might have been removed by anderson test)
if model_subset_for_uncertainty is ModelSubsetForUncertainty.non_stationary_gumbel_and_gev\
and self.select_only_model_that_pass_anderson_test:
massifs_names = massifs_names.intersection(set(self.massif_name_to_trend_test_that_minimized_aic.keys()))
arguments = [ arguments = [
[self.massif_name_to_years_and_maxima_for_model_fitting[m], [self.massif_name_to_years_and_maxima_for_model_fitting[m],
self.massif_name_and_model_subset_to_model_class(m, model_subset_for_uncertainty), self.massif_name_and_model_subset_to_model_class(m, model_subset_for_uncertainty),
......
...@@ -55,7 +55,7 @@ def max_graph_annual_maxima_poster(): ...@@ -55,7 +55,7 @@ def max_graph_annual_maxima_poster():
last_plot = color == "magenta" last_plot = color == "magenta"
label = '{} massif at {}m'.format(massif_name, altitude) label = '{} massif at {}m'.format(massif_name, altitude)
tight_pad = {'h_pad': 0.2} tight_pad = {'h_pad': 0.2}
snow_abbreviation = 'annual maxima of ' + snow_abbreviation snow_abbreviation = 'annual maximum of ' + snow_abbreviation
study_visualizer.visualize_max_graphs_poster(massif_name, altitude, snow_abbreviation, color, label, study_visualizer.visualize_max_graphs_poster(massif_name, altitude, snow_abbreviation, color, label,
last_plot, ax, tight_pad=tight_pad, last_plot, ax, tight_pad=tight_pad,
dpi=dpi_paper1_figure, dpi=dpi_paper1_figure,
......
...@@ -74,8 +74,10 @@ def intermediate_result(altitudes, massif_names=None, ...@@ -74,8 +74,10 @@ def intermediate_result(altitudes, massif_names=None,
plot_uncertainty_massifs(altitude_to_visualizer) plot_uncertainty_massifs(altitude_to_visualizer)
plot_uncertainty_histogram(altitude_to_visualizer) plot_uncertainty_histogram(altitude_to_visualizer)
plot_selection_curves(altitude_to_visualizer) plot_selection_curves(altitude_to_visualizer)
# uncertainty_interval_size(altitude_to_visualizer)
plot_intensity_against_gumbel_quantile_for_3_examples(altitude_to_visualizer) plot_intensity_against_gumbel_quantile_for_3_examples(altitude_to_visualizer)
# Additional plots
# uncertainty_interval_size(altitude_to_visualizer)
# plot_full_diagnostic(altitude_to_visualizer) # plot_full_diagnostic(altitude_to_visualizer)
...@@ -85,13 +87,9 @@ def major_result(): ...@@ -85,13 +87,9 @@ def major_result():
# massif_names = ['Beaufortain', 'Vercors'] # massif_names = ['Beaufortain', 'Vercors']
massif_names = None massif_names = None
study_classes = paper_study_classes[:1] study_classes = paper_study_classes[:1]
# model_subsets_for_uncertainty = [ModelSubsetForUncertainty.stationary_gumbel,
# ModelSubsetForUncertainty.stationary_gumbel_and_gev,
# ModelSubsetForUncertainty.non_stationary_gumbel,
# ModelSubsetForUncertainty.non_stationary_gumbel_and_gev]
model_subsets_for_uncertainty = None model_subsets_for_uncertainty = None
# study_classes = [CrocusSnowLoad3Days, CrocusSnowLoad5Days, CrocusSnowLoad7Days][::-1] # study_classes = [CrocusSnowLoad3Days, CrocusSnowLoad5Days, CrocusSnowLoad7Days][::-1]
altitudes = [300, 600] # altitudes = [300, 600, 900, 1800, 2700][:2]
altitudes = paper_altitudes altitudes = paper_altitudes
# altitudes = [900, 1800, 2700][:1] # altitudes = [900, 1800, 2700][:1]
for study_class in study_classes: for study_class in study_classes:
......
...@@ -19,7 +19,7 @@ from extreme_trend.trend_test_two_parameters.gumbel_test_two_parameters import \ ...@@ -19,7 +19,7 @@ from extreme_trend.trend_test_two_parameters.gumbel_test_two_parameters import \
GumbelLocationAndScaleTrendTest, GumbelLocationQuadraticTrendTest, GumbelScaleQuadraticTrendTest GumbelLocationAndScaleTrendTest, GumbelLocationQuadraticTrendTest, GumbelScaleQuadraticTrendTest
paper_altitudes = ALL_ALTITUDES_WITHOUT_NAN paper_altitudes = ALL_ALTITUDES_WITHOUT_NAN
paper_study_classes = [CrocusSnowLoadTotal, CrocusSnowLoadEurocode, CrocusSnowLoad3Days][:2] paper_study_classes = [CrocusSnowLoadTotal, CrocusSnowLoadEurocode][:]
# dpi_paper1_figure = 700 # dpi_paper1_figure = 700
dpi_paper1_figure = None dpi_paper1_figure = None
NON_STATIONARY_TREND_TEST_PAPER_1 = [GumbelVersusGumbel, NON_STATIONARY_TREND_TEST_PAPER_1 = [GumbelVersusGumbel,
......
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