[CONFIDENCE INTERVAL] add confidence interval test for the normal (i.e. delta)...

[CONFIDENCE INTERVAL] add confidence interval test for the normal (i.e. delta) confidence interval method

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

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:

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

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',
+}

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",

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)
+
+

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()