From eea97371b166e423739745d9dacb79ed7d072594 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Tue, 7 May 2019 18:08:58 +0200
Subject: [PATCH] [SCM][NON STATIONARITY] add non stationary trend object. to
create specific plot for this topic
---
.../main_study_visualizer.py | 11 ++-
.../non_stationary_trends.py | 79 +++++++++++++++++++
.../study_visualization/study_visualizer.py | 32 ++++++--
.../margin_model/linear_margin_model.py | 2 +-
.../parametric_margin_function.py | 1 -
.../extreme_models/result_from_fit.py | 19 ++++-
.../test_estimator/test_margin_estimators.py | 1 +
7 files changed, 135 insertions(+), 10 deletions(-)
create mode 100644 experiment/meteo_france_SCM_study/visualization/study_visualization/non_stationary_trends.py
diff --git a/experiment/meteo_france_SCM_study/visualization/study_visualization/main_study_visualizer.py b/experiment/meteo_france_SCM_study/visualization/study_visualization/main_study_visualizer.py
index 19a01c23..17ab7a99 100644
--- a/experiment/meteo_france_SCM_study/visualization/study_visualization/main_study_visualizer.py
+++ b/experiment/meteo_france_SCM_study/visualization/study_visualization/main_study_visualizer.py
@@ -95,9 +95,18 @@ def complete_analysis(only_first_one=False):
study_visualizer.visualize_linear_margin_fit()
+def trend_analysis():
+ save_to_file = False
+ only_first_one = True
+ for study_class in [CrocusDepth, SafranSnowfall, SafranRainfall, SafranTemperature][1:2]:
+ for study in study_iterator(study_class, only_first_one=only_first_one):
+ study_visualizer = StudyVisualizer(study, save_to_file=save_to_file)
+ study_visualizer.visualize_temporal_trend_relevance(complete_analysis=False)
+
if __name__ == '__main__':
# annual_mean_vizu_compare_durand_study(safran=True, take_mean_value=True, altitude=2100)
- normal_visualization(temporal_non_stationarity=True)
+ # normal_visualization(temporal_non_stationarity=True)
+ trend_analysis()
# max_stable_process_vizu_compare_gaume_study(altitude=1800, nb_days=1)
# extended_visualization()
# complete_analysis()
diff --git a/experiment/meteo_france_SCM_study/visualization/study_visualization/non_stationary_trends.py b/experiment/meteo_france_SCM_study/visualization/study_visualization/non_stationary_trends.py
new file mode 100644
index 00000000..89fe4351
--- /dev/null
+++ b/experiment/meteo_france_SCM_study/visualization/study_visualization/non_stationary_trends.py
@@ -0,0 +1,79 @@
+from extreme_estimator.estimator.abstract_estimator import AbstractEstimator
+from scipy.stats import chi2
+from extreme_estimator.estimator.full_estimator.abstract_full_estimator import \
+ FullEstimatorInASingleStepWithSmoothMargin
+from extreme_estimator.estimator.margin_estimator.abstract_margin_estimator import LinearMarginEstimator
+from extreme_estimator.extreme_models.margin_model.linear_margin_model import \
+ LinearAllParametersTwoFirstCoordinatesMarginModel, LinearAllTwoStatialCoordinatesLocationLinearMarginModel, \
+ LinearStationaryMarginModel, LinearNonStationaryLocationMarginModel
+from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
+
+
+class AbstractNonStationaryTrendTest(object):
+ RESULT_ATTRIBUTE_METRIC = 'deviance'
+
+ def __init__(self, dataset: AbstractDataset, estimator_class,
+ stationary_margin_model_class, non_stationary_margin_model_class):
+ self.dataset = dataset
+ self.estimator_class = estimator_class
+ self.stationary_margin_model_class = stationary_margin_model_class
+ self.non_stationary_margin_model_class = non_stationary_margin_model_class
+
+ def load_estimator(self, margin_model) -> AbstractEstimator:
+ return self.estimator_class(self.dataset, margin_model)
+
+ def get_metric(self, margin_model_class, starting_point):
+ margin_model = margin_model_class(coordinates=self.dataset.coordinates, starting_point=starting_point)
+ estimator = self.load_estimator(margin_model) # type: AbstractEstimator
+ estimator.fit()
+ metric = estimator.result_from_fit.__getattribute__(self.RESULT_ATTRIBUTE_METRIC)
+ assert isinstance(metric, float)
+ return metric
+
+ def visualize(self, ax, complete_analysis=True):
+ # Define the year_min and year_max for the starting point
+ if complete_analysis:
+ year_min, year_max, step = 1960, 1990, 1
+ else:
+ year_min, year_max, step = 1960, 1990, 10
+ # Fit the stationary model
+ stationary_metric = self.get_metric(self.stationary_margin_model_class, starting_point=None)
+ # Fit the non stationary model
+ years = list(range(year_min, year_max + 1, step))
+ non_stationary_metrics = [self.get_metric(self.non_stationary_margin_model_class, year) for year in years]
+ difference = [m - stationary_metric for m in non_stationary_metrics]
+ # Plot some lines
+ ax.axhline(y=0, xmin=year_min, xmax=year_max)
+ # Significative line
+ significative_deviance = chi2.ppf(q=0.95, df=1)
+ ax.axhline(y=significative_deviance, xmin=year_min, xmax=year_max)
+ # todo: plot the line corresponding to the significance 0.05
+ ax.plot(years, difference, 'ro-')
+
+
+class IndependenceLocationTrendTest(AbstractNonStationaryTrendTest):
+
+ def __init__(self, dataset, coordinate_idx):
+ pass
+
+
+class ConditionalIndedendenceLocationTrendTest(AbstractNonStationaryTrendTest):
+
+ def __init__(self, dataset):
+ super().__init__(dataset=dataset,
+ estimator_class=LinearMarginEstimator,
+ stationary_margin_model_class=LinearStationaryMarginModel,
+ non_stationary_margin_model_class=LinearNonStationaryLocationMarginModel)
+
+
+class MaxStableLocationTrendTest(AbstractNonStationaryTrendTest):
+
+ def __init__(self, dataset, max_stable_model):
+ super().__init__(dataset=dataset,
+ estimator_class=FullEstimatorInASingleStepWithSmoothMargin,
+ stationary_margin_model_class=LinearStationaryMarginModel,
+ non_stationary_margin_model_class=LinearNonStationaryLocationMarginModel)
+ self.max_stable_model = max_stable_model
+
+ def load_estimator(self, margin_model) -> AbstractEstimator:
+ return self.estimator_class(self.dataset, margin_model, self.max_stable_model)
diff --git a/experiment/meteo_france_SCM_study/visualization/study_visualization/study_visualizer.py b/experiment/meteo_france_SCM_study/visualization/study_visualization/study_visualizer.py
index 94e46bcd..ea9d5b5c 100644
--- a/experiment/meteo_france_SCM_study/visualization/study_visualization/study_visualizer.py
+++ b/experiment/meteo_france_SCM_study/visualization/study_visualization/study_visualizer.py
@@ -9,13 +9,15 @@ import pandas as pd
import seaborn as sns
from experiment.meteo_france_SCM_study.abstract_study import AbstractStudy
+from experiment.meteo_france_SCM_study.visualization.study_visualization.non_stationary_trends import \
+ ConditionalIndedendenceLocationTrendTest, MaxStableLocationTrendTest
from experiment.meteo_france_SCM_study.visualization.utils import create_adjusted_axes
from experiment.utils import average_smoothing_with_sliding_window
from extreme_estimator.estimator.full_estimator.abstract_full_estimator import \
FullEstimatorInASingleStepWithSmoothMargin
from extreme_estimator.estimator.margin_estimator.abstract_margin_estimator import LinearMarginEstimator
from extreme_estimator.extreme_models.margin_model.linear_margin_model import LinearAllParametersAllDimsMarginModel, \
- LinearNonStationaryMarginModel, LinearStationaryMarginModel
+ LinearNonStationaryLocationMarginModel, LinearStationaryMarginModel
from extreme_estimator.extreme_models.margin_model.margin_function.abstract_margin_function import \
AbstractMarginFunction
from extreme_estimator.extreme_models.margin_model.param_function.param_function import AbstractParamFunction
@@ -51,6 +53,8 @@ class StudyVisualizer(object):
self._coordinates = None
self._observations = None
+ self.default_covariance_function = CovarianceFunction.powexp
+
# KDE PLOT ARGUMENTS
self.vertical_kde_plot = vertical_kde_plot
self.year_for_kde_plot = year_for_kde_plot
@@ -132,6 +136,25 @@ class StudyVisualizer(object):
'for the year {}'.format(self.year_for_kde_plot)
self.show_or_save_to_file()
+ def visualize_temporal_trend_relevance(self, complete_analysis):
+ self.temporal_non_stationarity = True
+ trend_tests = [ConditionalIndedendenceLocationTrendTest(self.dataset)]
+
+ max_stable_models = load_test_max_stable_models(default_covariance_function=self.default_covariance_function)
+ for max_stable_model in max_stable_models[:1]:
+ trend_tests.append(MaxStableLocationTrendTest(self.dataset, max_stable_model))
+
+ nb_trend_tests = len(trend_tests)
+ fig, axes = plt.subplots(1, nb_trend_tests, figsize=self.figsize)
+ if nb_trend_tests == 1:
+ axes = [axes]
+ fig.subplots_adjust(hspace=self.subplot_space, wspace=self.subplot_space)
+ for ax, trend_test in zip(axes, trend_tests):
+ trend_test.visualize(ax, complete_analysis=complete_analysis)
+
+ self.plot_name = 'trend tests'
+ self.show_or_save_to_file()
+
def visualize_experimental_law(self, ax, massif_id):
# Display the experimental law for a given massif
all_massif_data = self.get_all_massif_data(massif_id)
@@ -248,15 +271,14 @@ class StudyVisualizer(object):
pass
def visualize_linear_margin_fit(self, only_first_max_stable=False):
- default_covariance_function = CovarianceFunction.powexp
- margin_class = LinearNonStationaryMarginModel if self.temporal_non_stationarity else LinearStationaryMarginModel
+ margin_class = LinearNonStationaryLocationMarginModel if self.temporal_non_stationarity else LinearStationaryMarginModel
plot_name = 'Full Likelihood with Linear marginals and max stable dependency structure'
plot_name += '\n(with {} covariance structure when a covariance is needed)'.format(
- str(default_covariance_function).split('.')[-1])
+ str(self.default_covariance_function).split('.')[-1])
self.plot_name = plot_name
# Load max stable models
- max_stable_models = load_test_max_stable_models(default_covariance_function=default_covariance_function)
+ max_stable_models = load_test_max_stable_models(default_covariance_function=self.default_covariance_function)
if only_first_max_stable:
# Keep only the BrownResnick model
max_stable_models = max_stable_models[1:2]
diff --git a/extreme_estimator/extreme_models/margin_model/linear_margin_model.py b/extreme_estimator/extreme_models/margin_model/linear_margin_model.py
index adbb16d9..14b1b7d9 100644
--- a/extreme_estimator/extreme_models/margin_model/linear_margin_model.py
+++ b/extreme_estimator/extreme_models/margin_model/linear_margin_model.py
@@ -127,5 +127,5 @@ class LinearStationaryMarginModel(LinearAllParametersTwoFirstCoordinatesMarginMo
pass
-class LinearNonStationaryMarginModel(LinearAllTwoStatialCoordinatesLocationLinearMarginModel):
+class LinearNonStationaryLocationMarginModel(LinearAllTwoStatialCoordinatesLocationLinearMarginModel):
pass
diff --git a/extreme_estimator/extreme_models/margin_model/margin_function/parametric_margin_function.py b/extreme_estimator/extreme_models/margin_model/margin_function/parametric_margin_function.py
index 7476aa13..09337e6d 100644
--- a/extreme_estimator/extreme_models/margin_model/margin_function/parametric_margin_function.py
+++ b/extreme_estimator/extreme_models/margin_model/margin_function/parametric_margin_function.py
@@ -62,7 +62,6 @@ class ParametricMarginFunction(IndependentMarginFunction):
raise NotImplementedError
def get_gev_params(self, coordinate: np.ndarray) -> GevParams:
- print('here get gev', self.starting_point)
if self.starting_point is not None:
# Shift temporal coordinate to enable to model temporal trend with starting point
assert self.coordinates.has_temporal_coordinates
diff --git a/extreme_estimator/extreme_models/result_from_fit.py b/extreme_estimator/extreme_models/result_from_fit.py
index bc01dd8b..669c2f87 100644
--- a/extreme_estimator/extreme_models/result_from_fit.py
+++ b/extreme_estimator/extreme_models/result_from_fit.py
@@ -1,5 +1,6 @@
from typing import Dict
+import numpy as np
from rpy2 import robjects
from extreme_estimator.extreme_models.margin_model.param_function.linear_coef import LinearCoef
@@ -9,6 +10,7 @@ from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoo
class ResultFromFit(object):
+
def __init__(self, result_from_fit: robjects.ListVector) -> None:
if hasattr(result_from_fit, 'names'):
self.name_to_value = {name: result_from_fit.rx2(name) for name in result_from_fit.names}
@@ -27,6 +29,14 @@ class ResultFromFit(object):
def margin_coef_dict(self):
raise NotImplementedError
+ @property
+ def nllh(self):
+ raise NotImplementedError
+
+ @property
+ def deviance(self):
+ raise NotImplementedError
+
class ResultFromIsmev(ResultFromFit):
@@ -47,7 +57,8 @@ class ResultFromIsmev(ResultFromFit):
i += 1
# Add a potential linear temporal trend
if gev_param_name in self.gev_param_name_to_dim:
- temporal_coef_name = LinearCoef.coef_template_str(gev_param_name, AbstractCoordinates.COORDINATE_T).format(1)
+ temporal_coef_name = LinearCoef.coef_template_str(gev_param_name,
+ AbstractCoordinates.COORDINATE_T).format(1)
coef_dict[temporal_coef_name] = mle_values[i]
i += 1
return coef_dict
@@ -58,7 +69,7 @@ class ResultFromIsmev(ResultFromFit):
@property
def nllh(self):
- return self.res['nllh']
+ return self.name_to_value['nllh']
class ResultFromSpatialExtreme(ResultFromFit):
@@ -68,6 +79,10 @@ class ResultFromSpatialExtreme(ResultFromFit):
FITTED_VALUES_NAME = 'fitted.values'
CONVERGENCE_NAME = 'convergence'
+ @property
+ def deviance(self):
+ return np.array(self.name_to_value['deviance'])[0]
+
@property
def convergence(self) -> str:
convergence_value = self.name_to_value[self.CONVERGENCE_NAME]
diff --git a/test/test_extreme_estimator/test_estimator/test_margin_estimators.py b/test/test_extreme_estimator/test_estimator/test_margin_estimators.py
index 309d8701..11f2813f 100644
--- a/test/test_extreme_estimator/test_estimator/test_margin_estimators.py
+++ b/test/test_extreme_estimator/test_estimator/test_margin_estimators.py
@@ -30,6 +30,7 @@ class TestSmoothMarginEstimator(unittest.TestCase):
# Fit estimator
estimator = LinearMarginEstimator(dataset=dataset, margin_model=margin_model)
estimator.fit()
+ print(estimator.result_from_fit.name_to_value.keys())
# Plot
if self.DISPLAY:
margin_model.margin_function_sample.visualize_function(show=True)
--
GitLab