From 198e83295ee4ff1c5b8746ce498d9c468330f2e9 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Tue, 19 Mar 2019 15:45:02 +0100
Subject: [PATCH] [EXTREME MODEL] add max_stable_with_spline.R & prepare test
---
.../main_study_visualizer.py | 8 ++-
.../study_visualization/study_visualizer.py | 10 +++-
.../max_stable_model/max_stable_with_spline.R | 54 +++++++++++++++++++
.../test_fitmaxstab_with_margin.py | 36 ++++++++++++-
4 files changed, 105 insertions(+), 3 deletions(-)
create mode 100644 extreme_estimator/extreme_models/max_stable_model/max_stable_with_spline.R
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 5b9dbe4f..0ec421ae 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
@@ -62,6 +62,11 @@ def annual_mean_vizu_compare_durand_study(safran=True, take_mean_value=True, alt
study_visualizer.visualize_annual_mean_values(take_mean_value=take_mean_value)
+def max_stable_process_vizu_compare_gaume_study(altitude=1800, nb_days=1):
+ study = SafranSnowfall(altitude=altitude, nb_consecutive_days=nb_days)
+ study_visualizer = StudyVisualizer(study)
+ study_visualizer.visualize_brown_resnick_fit()
+
def normal_visualization():
save_to_file = False
only_first_one = True
@@ -91,6 +96,7 @@ def complete_analysis(only_first_one=False):
if __name__ == '__main__':
# annual_mean_vizu_compare_durand_study(safran=True, take_mean_value=True, altitude=2400)
- normal_visualization()
+ # normal_visualization()
+ 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/study_visualizer.py b/experiment/meteo_france_SCM_study/visualization/study_visualization/study_visualizer.py
index a8741c2b..fe6ab50a 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
@@ -16,6 +16,7 @@ from extreme_estimator.estimator.margin_estimator.abstract_margin_estimator impo
from extreme_estimator.extreme_models.margin_model.param_function.param_function import ParamFunction
from extreme_estimator.extreme_models.margin_model.smooth_margin_model import LinearAllParametersAllDimsMarginModel
from extreme_estimator.extreme_models.max_stable_model.abstract_max_stable_model import CovarianceFunction
+from extreme_estimator.extreme_models.max_stable_model.max_stable_models import BrownResnick
from extreme_estimator.margin_fits.abstract_params import AbstractParams
from extreme_estimator.margin_fits.gev.gev_params import GevParams
from extreme_estimator.margin_fits.gev.gevmle_fit import GevMleFit
@@ -209,6 +210,10 @@ class StudyVisualizer(object):
ax.set_xlabel('year')
ax.set_title(self.study.safran_massif_names[massif_id])
+ def visualize_brown_resnick_fit(self):
+ pass
+
+
def visualize_linear_margin_fit(self, only_first_max_stable=False):
default_covariance_function = CovarianceFunction.cauchy
plot_name = 'Full Likelihood with Linear marginals and max stable dependency structure'
@@ -217,9 +222,12 @@ class StudyVisualizer(object):
self.plot_name = plot_name
max_stable_models = load_test_max_stable_models(default_covariance_function=default_covariance_function)
if only_first_max_stable:
- max_stable_models = max_stable_models[:1]
+ # Keep only the BrownResnick model
+ max_stable_models = max_stable_models[1:2]
if only_first_max_stable is None:
max_stable_models = []
+
+
fig, axes = plt.subplots(len(max_stable_models) + 2, len(GevParams.SUMMARY_NAMES), figsize=self.figsize)
fig.subplots_adjust(hspace=self.subplot_space, wspace=self.subplot_space)
margin_class = LinearAllParametersAllDimsMarginModel
diff --git a/extreme_estimator/extreme_models/max_stable_model/max_stable_with_spline.R b/extreme_estimator/extreme_models/max_stable_model/max_stable_with_spline.R
new file mode 100644
index 00000000..d074d5c4
--- /dev/null
+++ b/extreme_estimator/extreme_models/max_stable_model/max_stable_with_spline.R
@@ -0,0 +1,54 @@
+library(SpatialExtremes)
+
+## Not run:
+##Define the coordinate of each location
+n.site <- 30
+locations <- matrix(runif(2*n.site, 0, 10), ncol = 2)
+colnames(locations) <- c("lon", "lat")
+
+## Fitting a penalized smoothing splines for the margins with the
+## Smith's model
+data <- rmaxstab(100, locations, cov.mod = "gauss", cov11 = 100, cov12 =
+ 25, cov22 = 220)
+
+## And transform it to ordinary GEV margins with a non-linear
+## function
+fun <- function(x)
+ 2 * sin(pi * x / 4) + 10
+fun2 <- function(x)
+ (fun(x) - 7 ) / 15
+
+param.loc <- fun(locations[,2])
+param.scale <- fun(locations[,2])
+param.shape <- fun2(locations[,1])
+
+##Transformation from unit Frechet to common GEV margins
+for (i in 1:n.site)
+ data[,i] <- frech2gev(data[,i], param.loc[i], param.scale[i],
+param.shape[i])
+
+##Defining the knots, penalty, degree for the splines
+n.knots_x = 2
+n.knots_y = 2
+knots = quantile(locations[,1], prob=1:n.knots_x/(n.knots_x+1))
+knots2 = quantile(locations[,2], prob=1:n.knots_y/(n.knots_y+1))
+knots_tot = cbind(knots,knots2)
+print(knots)
+print(knots2)
+print(knots_tot)
+
+##Be careful the choice of the penalty (i.e. the smoothing parameter)
+##may strongly affect the result Here we use p-splines for each GEV
+##parameter - so it's really CPU demanding but one can use 1 p-spline
+##and 2 linear models.
+##A simple linear model will be clearly faster...
+loc.form <- y ~ rb(locations[,1], knots = knots, degree = 3, penalty = .5)
+scale.form <- y ~ rb(locations[,2], knots = knots2, degree = 3, penalty = .5)
+shape.form <- y ~ rb(locations, knots = knots_tot, degree = 3, penalty = .5)
+
+fitted <- fitmaxstab(data, locations, "gauss", loc.form, scale.form, shape.form,
+ method = "BFGS")
+fitted
+
+
+## End(Not run)
\ No newline at end of file
diff --git a/test/test_unitary/test_fitmaxstab/test_fitmaxstab_with_margin.py b/test/test_unitary/test_fitmaxstab/test_fitmaxstab_with_margin.py
index a4e29995..0598d4e2 100644
--- a/test/test_unitary/test_fitmaxstab/test_fitmaxstab_with_margin.py
+++ b/test/test_unitary/test_fitmaxstab/test_fitmaxstab_with_margin.py
@@ -1,6 +1,7 @@
import unittest
-from extreme_estimator.estimator.full_estimator.abstract_full_estimator import FullEstimatorInASingleStepWithSmoothMargin
+from extreme_estimator.estimator.full_estimator.abstract_full_estimator import \
+ FullEstimatorInASingleStepWithSmoothMargin
from extreme_estimator.extreme_models.margin_model.smooth_margin_model import ConstantMarginModel, \
LinearMarginModelExample
from extreme_estimator.extreme_models.max_stable_model.abstract_max_stable_model import CovarianceFunction
@@ -65,5 +66,38 @@ class TestMaxStableFitWithLinearMargin(TestUnitaryAbstract):
self.compare()
+# class TestMaxStableFitWithSpline(TestUnitaryAbstract):
+#
+# @property
+# def r_output(self):
+# TestRMaxStabWithMarginConstant.r_code()
+# r("""
+# # Code inspired from Johan code
+# n.knots_x = 1
+# n.knots_y = 2
+# knots = quantile(locations[,1], prob=1:n.knots_x/(n.knots_x+1))
+# knots2 = quantile(locations[,2], prob=1:n.knots_y/(n.knots_y+1))
+# knots_tot = cbind(knots,knots2)
+#
+# loc.form <- y ~ rb(locations[,1], knots = knots, degree = 3, penalty = .5)
+# scale.form <- y ~ rb(locations[,2], knots = knots2, degree = 3, penalty = .5)
+# shape.form <- y ~ rb(locations, knots = knots_tot, degree = 3, penalty = .5)
+# """)
+# return self.r_fitted_values_from_res_variable
+#
+# @property
+# def python_output(self):
+# dataset = TestRMaxStabWithMarginConstant.python_code()
+# max_stable_model = Schlather(covariance_function=CovarianceFunction.whitmat, use_start_value=False)
+# margin_model = LinearMarginModelExample(dataset.coordinates)
+# full_estimator = FullEstimatorInASingleStepWithSmoothMargin(dataset, margin_model,
+# max_stable_model)
+# full_estimator.fit()
+# return full_estimator.fitted_values
+#
+# def test_max_stable_fit_with_spline_margin(self):
+# self.compare()
+
+
if __name__ == '__main__':
unittest.main()
--
GitLab