From 864f2e1b65525a77a8ba4419775f208c14c59dde Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Wed, 17 Jul 2019 13:40:04 +0200
Subject: [PATCH] [STRENGTH EVOLUTION] fix strength computation for the scale
test.
---
.../main1_good_stationary_gev_fit.py | 2 +-
...main3_non_stationary_strength_evolution.py | 23 ++++++++++---------
experiment/paper1_steps/utils.py | 6 ++---
.../abstract_gev_change_point_test.py | 5 ++++
.../margin_fits/gev/gev_params.py | 3 ++-
5 files changed, 23 insertions(+), 16 deletions(-)
diff --git a/experiment/paper1_steps/hard extreme evolution - annual maxima/main1_good_stationary_gev_fit.py b/experiment/paper1_steps/hard extreme evolution - annual maxima/main1_good_stationary_gev_fit.py
index 4bcf26cb..0cfb9bcb 100644
--- a/experiment/paper1_steps/hard extreme evolution - annual maxima/main1_good_stationary_gev_fit.py
+++ b/experiment/paper1_steps/hard extreme evolution - annual maxima/main1_good_stationary_gev_fit.py
@@ -6,7 +6,7 @@ from experiment.meteo_france_data.scm_models_data.visualization.study_visualizat
def maxima_analysis():
- save_to_file = True
+ save_to_file = False
only_first_one = False
durand_altitude = [900, 1500, 2100, 2700]
altitudes = durand_altitude
diff --git a/experiment/paper1_steps/hard extreme evolution - annual maxima/main3_non_stationary_strength_evolution.py b/experiment/paper1_steps/hard extreme evolution - annual maxima/main3_non_stationary_strength_evolution.py
index eccbb692..4cbedf1e 100644
--- a/experiment/paper1_steps/hard extreme evolution - annual maxima/main3_non_stationary_strength_evolution.py
+++ b/experiment/paper1_steps/hard extreme evolution - annual maxima/main3_non_stationary_strength_evolution.py
@@ -2,6 +2,8 @@ import time
from experiment.meteo_france_data.scm_models_data.visualization.hypercube_visualization.altitude_year_hypercube_visualizer import \
Altitude_Hypercube_Year_Visualizer
+from experiment.trend_analysis.univariate_test.abstract_gev_change_point_test import GevLocationChangePointTest, \
+ GevScaleChangePointTest
"""
Visualize the 0.99 quantile initial value and its evolution
@@ -12,25 +14,24 @@ from experiment.paper1_steps.utils import get_full_altitude_visualizer, FULL_ALT
def main_fast_spatial_risk_evolution():
for altitude in FULL_ALTITUDES[-1:]:
vizualiser = get_full_altitude_visualizer(Altitude_Hypercube_Year_Visualizer, altitude=altitude,
- exact_starting_year=1958, reduce_strength_array=True)
+ exact_starting_year=1958, reduce_strength_array=True,
+ trend_test_class=GevScaleChangePointTest)
vizualiser.save_to_file = False
vizualiser.visualize_massif_trend_test_one_altitude()
-def main_full_spatial_repartition():
+def main_full_spatial_risk_evolution():
for altitude in FULL_ALTITUDES[:]:
- # Compute for the most likely starting year
- # vizualiser = get_full_altitude_visualizer(Altitude_Hypercube_Year_Visualizer, altitude=altitude)
- # vizualiser.visualize_massif_trend_test_one_altitude()
- # Compute the trend for a linear trend
- vizualiser = get_full_altitude_visualizer(Altitude_Hypercube_Year_Visualizer, altitude=altitude,
- exact_starting_year=1958)
- vizualiser.visualize_massif_trend_test_one_altitude()
+ for trend_test_class in [GevLocationChangePointTest, GevScaleChangePointTest][1:]:
+ vizualiser = get_full_altitude_visualizer(Altitude_Hypercube_Year_Visualizer, altitude=altitude,
+ exact_starting_year=1958, reduce_strength_array=True,
+ trend_test_class=trend_test_class)
+ vizualiser.visualize_massif_trend_test_one_altitude()
def main_run():
- # main_full_spatial_repartition()
- main_fast_spatial_risk_evolution()
+ main_full_spatial_risk_evolution()
+ # main_fast_spatial_risk_evolution()
if __name__ == '__main__':
diff --git a/experiment/paper1_steps/utils.py b/experiment/paper1_steps/utils.py
index 699d710f..3aa74c75 100644
--- a/experiment/paper1_steps/utils.py
+++ b/experiment/paper1_steps/utils.py
@@ -14,13 +14,13 @@ FULL_ALTITUDES = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000]
def get_full_altitude_visualizer(altitude_hypercube_class, exact_starting_year=None, altitude=900,
- reduce_strength_array=False):
- altitudes, first_starting_year, last_starting_year, nb_data_reduced_for_speed, only_first_one, save_to_file, trend_test_class = get_full_parameters(
+ reduce_strength_array=False,
+ trend_test_class = GevLocationChangePointTest):
+ altitudes, first_starting_year, last_starting_year, nb_data_reduced_for_speed, only_first_one, save_to_file, _ = get_full_parameters(
altitude=altitude)
if exact_starting_year is not None:
first_starting_year, last_starting_year = None, None
study_classes = [CrocusRecentSwe]
- trend_test_class = GevLocationChangePointTest
visualizer = load_altitude_visualizer(altitude_hypercube_class, altitudes, last_starting_year,
nb_data_reduced_for_speed, only_first_one, save_to_file, study_classes,
trend_test_class, first_starting_year=first_starting_year,
diff --git a/experiment/trend_analysis/univariate_test/abstract_gev_change_point_test.py b/experiment/trend_analysis/univariate_test/abstract_gev_change_point_test.py
index f7acb27e..6417491e 100644
--- a/experiment/trend_analysis/univariate_test/abstract_gev_change_point_test.py
+++ b/experiment/trend_analysis/univariate_test/abstract_gev_change_point_test.py
@@ -169,6 +169,11 @@ class GevScaleChangePointTest(AbstractGevChangePointTest):
super().__init__(years, maxima, starting_year,
NonStationaryScaleStationModel, GevParams.SCALE)
+ def _slope_strength(self):
+ return self.non_stationary_constant_gev_params.quantile_strength_evolution_ratio(
+ p=self.quantile_for_strength,
+ sigma1=self.non_stationary_linear_coef)
+
class GevShapeChangePointTest(AbstractGevChangePointTest):
diff --git a/extreme_estimator/margin_fits/gev/gev_params.py b/extreme_estimator/margin_fits/gev/gev_params.py
index 90f6f946..650d4dca 100644
--- a/extreme_estimator/margin_fits/gev/gev_params.py
+++ b/extreme_estimator/margin_fits/gev/gev_params.py
@@ -58,7 +58,8 @@ class GevParams(ExtremeParams):
initial_quantile = self.quantile(p)
quantity_increased = mu1
if sigma1 != 0:
- quantity_increased += (sigma1 / self.shape) * (1 - (- np.float_power(np.log(p), -self.shape)))
+ power = np.float_power(- np.log(p), -self.shape)
+ quantity_increased -= (sigma1 / self.shape) * (1 - power)
return quantity_increased / initial_quantile
# Compute some indicators (such as the mean and the variance)
--
GitLab