From 743308abb34467c6e029e4f8439ab6e47f618bd0 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Thu, 31 Oct 2019 16:14:45 +0100
Subject: [PATCH] [CONFIDENCE INTERVAL] add confidence interval test for the
normal (i.e. delta) confidence interval method
---
extreme_fit/distribution/gev/main_fevd_mle.R | 17 +++++++--------
.../abstract_extract_eurocode_return_level.py | 3 ++-
.../abstract_result_from_extremes.py | 8 +++----
.../confidence_interval_method.py | 7 +++++++
.../result_from_bayesian_extremes.py | 2 +-
.../result_from_mle_extremes.py | 21 +++++++++++++++++++
.../test_model/test_confidence_interval.py | 13 ++++++++++++
7 files changed, 55 insertions(+), 16 deletions(-)
diff --git a/extreme_fit/distribution/gev/main_fevd_mle.R b/extreme_fit/distribution/gev/main_fevd_mle.R
index 5a83f9b4..96901909 100644
--- a/extreme_fit/distribution/gev/main_fevd_mle.R
+++ b/extreme_fit/distribution/gev/main_fevd_mle.R
@@ -9,18 +9,15 @@ library(SpatialExtremes)
source('fevd_fixed.R')
# Sample from a GEV
set.seed(42)
-N <- 100
-loc = 0; scale = 1; shape <- 0.1
+N <- 50
+loc = 0; scale = 1; shape <- 1
x_gev <- rgev(N, loc = loc, scale = scale, shape = shape)
-
-
-# Add covariate
coord <- matrix(ncol=1, nrow = N)
-coord[,1]=seq(1,N,1)
+coord[,1]=seq(0,N-1,1)
colnames(coord) = c("T")
coord = data.frame(coord, stringsAsFactors = TRUE)
-# res = fevd_fixed(x_gev, data=coord, method='MLE', verbose=TRUE, use.phi=FALSE)
-res = fevd_fixed(x_gev, data=coord, location.fun= ~T, method='MLE', verbose=TRUE, use.phi=FALSE, time.units = "years", units = "years")
+res = fevd_fixed(x_gev, data=coord, method='MLE', verbose=TRUE, use.phi=FALSE)
+# res = fevd_fixed(x_gev, data=coord, location.fun= ~T, method='MLE', verbose=TRUE, use.phi=FALSE, time.units = "years", units = "years")
# print(res)
# Some display for the results
@@ -34,10 +31,10 @@ print(res$results$par)
# Confidence interval staionary
-# method = "proflik"
+method = "normal"
res_ci = ci(res, alpha = 0.05, type = c("return.level", "parameter"),
return.period = 50, method = method, xrange = NULL, nint = 20, verbose = FALSE,
- tscale = FALSE, return.samples = FALSE)
+ tscale = FALSE)
print(res_ci)
# Bug to solve for the non stationary - the returned parameter do not match with the return level
diff --git a/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_extract_eurocode_return_level.py b/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_extract_eurocode_return_level.py
index f47b10dd..c0cd76d2 100644
--- a/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_extract_eurocode_return_level.py
+++ b/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_extract_eurocode_return_level.py
@@ -45,7 +45,8 @@ class ExtractEurocodeReturnLevelFromCiMethod(AbstractExtractEurocodeReturnLevel)
def confidence_interval_method(self):
return self.result_from_fit.confidence_interval_method(self.eurocode_quantile,
self.alpha_for_confidence_interval,
- self.transformed_temporal_covariate)
+ self.transformed_temporal_covariate,
+ self.ci_method)
@property
def mean_estimate(self) -> np.ndarray:
diff --git a/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_result_from_extremes.py b/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_result_from_extremes.py
index 0a93a88f..cbd16d9a 100644
--- a/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_result_from_extremes.py
+++ b/extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_result_from_extremes.py
@@ -16,13 +16,13 @@ class AbstractResultFromExtremes(AbstractResultFromModelFit):
@property
def is_non_stationary(self):
- return len(self.gev_param_name_to_dim) == 0
+ return len(self.gev_param_name_to_dim) > 0
def load_dataframe_from_r_matrix(self, name):
r_matrix = self.name_to_value[name]
return pd.DataFrame(np.array(r_matrix), columns=r.colnames(r_matrix))
- def confidence_interval_method(self, quantile_level, alpha_interval, transformed_temporal_covariate):
+ def confidence_interval_method(self, quantile_level, alpha_interval, transformed_temporal_covariate, ci_method):
return_period = round(1 / (1 - quantile_level))
common_kwargs = {
'return.period': return_period,
@@ -40,8 +40,8 @@ class AbstractResultFromExtremes(AbstractResultFromModelFit):
qcov = r("make.qcov")(self.result_from_fit,
**kwargs)
common_kwargs['qcov'] = qcov
- mean_estimate, confidence_interval = self._confidence_interval_method(common_kwargs)
+ mean_estimate, confidence_interval = self._confidence_interval_method(common_kwargs, ci_method)
return mean_estimate, confidence_interval
- def _confidence_interval_method(self, common_kwargs):
+ def _confidence_interval_method(self, common_kwargs, ci_method):
raise NotImplementedError
diff --git a/extreme_fit/model/result_from_model_fit/result_from_extremes/confidence_interval_method.py b/extreme_fit/model/result_from_model_fit/result_from_extremes/confidence_interval_method.py
index 3fb5da38..e151009e 100644
--- a/extreme_fit/model/result_from_model_fit/result_from_extremes/confidence_interval_method.py
+++ b/extreme_fit/model/result_from_model_fit/result_from_extremes/confidence_interval_method.py
@@ -9,3 +9,10 @@ class ConfidenceIntervalMethodFromExtremes(Enum):
ci_proflik = 3
# Confidence interval from my functions
my_bayes = 4
+
+
+ci_method_to_method_name = {
+ ConfidenceIntervalMethodFromExtremes.ci_normal: 'normal',
+ ConfidenceIntervalMethodFromExtremes.ci_boot: 'boot',
+ ConfidenceIntervalMethodFromExtremes.ci_proflik: 'proflik',
+}
diff --git a/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_bayesian_extremes.py b/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_bayesian_extremes.py
index 851d0837..b83915d8 100644
--- a/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_bayesian_extremes.py
+++ b/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_bayesian_extremes.py
@@ -47,7 +47,7 @@ class ResultFromBayesianExtremes(AbstractResultFromExtremes):
mean_posterior_parameters = self.df_posterior_samples.iloc[:, :-2].mean(axis=0)
return self.get_coef_dict_from_posterior_sample(mean_posterior_parameters)
- def _confidence_interval_method(self, common_kwargs):
+ def _confidence_interval_method(self, common_kwargs, ci_method):
bayesian_ci_parameters = {
'burn.in': self.burn_in_nb,
'FUN': "mean",
diff --git a/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_mle_extremes.py b/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_mle_extremes.py
index 079c4d63..08bbd976 100644
--- a/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_mle_extremes.py
+++ b/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_mle_extremes.py
@@ -2,7 +2,10 @@ import numpy as np
from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_result_from_extremes import \
AbstractResultFromExtremes
+from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import \
+ ci_method_to_method_name
from extreme_fit.model.result_from_model_fit.utils import get_margin_coef_ordered_dict
+from extreme_fit.model.utils import r
class ResultFromMleExtremes(AbstractResultFromExtremes):
@@ -13,3 +16,21 @@ class ResultFromMleExtremes(AbstractResultFromExtremes):
d = self.get_python_dictionary(values)
values = {i: param for i, param in enumerate(np.array(d['par']))}
return get_margin_coef_ordered_dict(self.gev_param_name_to_dim, values)
+
+ def _confidence_interval_method(self, common_kwargs, ci_method):
+ method_name = ci_method_to_method_name[ci_method]
+ mle_ci_parameters = {
+ 'method': method_name,
+ # xrange = NULL, nint = 20
+ }
+ res = r.ci(self.result_from_fit, **mle_ci_parameters, **common_kwargs)
+ if self.is_non_stationary:
+ a = np.array(res)[0]
+ lower, mean_estimate, upper, _ = a
+ else:
+ d = self.get_python_dictionary(res)
+ keys = ['50-year return level', '95% lower CI', '95% upper CI']
+ mean_estimate, lower, upper = [np.array(d[k])[0] for k in keys]
+ return mean_estimate, (lower, upper)
+
+
diff --git a/test/test_extreme_fit/test_model/test_confidence_interval.py b/test/test_extreme_fit/test_model/test_confidence_interval.py
index b7b4c528..55cc6290 100644
--- a/test/test_extreme_fit/test_model/test_confidence_interval.py
+++ b/test/test_extreme_fit/test_model/test_confidence_interval.py
@@ -72,6 +72,15 @@ class TestConfidenceInterval(unittest.TestCase):
self.ci_method = ConfidenceIntervalMethodFromExtremes.ci_bayes
self.model_class_to_triplet = self.bayesian_ci
+ def test_ci_normal(self):
+ self.fit_method = TemporalMarginFitMethod.extremes_fevd_mle
+ self.ci_method = ConfidenceIntervalMethodFromExtremes.ci_normal
+ self.model_class_to_triplet= {
+ StationaryTemporalModel: (-4.703945484843988, 30.482318639674023, 65.66858276419204),
+ NonStationaryLocationTemporalModel: (-4.223086740397132, 30.29842988666537, 64.81994651372787),
+ NonStationaryLocationAndScaleTemporalModel: (-15.17041284612494, 43.69511224410276, 102.56063733433047),
+ }
+
def tearDown(self) -> None:
for model_class, expected_triplet in self.model_class_to_triplet.items():
eurocode_ci = self.compute_eurocode_ci(model_class)
@@ -100,6 +109,10 @@ class TestConfidenceIntervalModifiedCoordinates(TestConfidenceInterval):
def test_ci_bayes(self):
super().test_ci_bayes()
+ def test_ci_normal(self):
+ self.model_class_to_triplet = {}
+ self.assertTrue(True)
+
if __name__ == '__main__':
unittest.main()
--
GitLab