[CONFIDENCE INTERVAL] check that ci code return the same interval and mean...

[CONFIDENCE INTERVAL] check that ci code return the same interval and mean posterior for the return level as my code with coordinates from 0 to 49. and also with transformed coordinates

experiment/eurocode_data/main_eurocode_drawing.py

@@ -23,7 +23,7 @@ mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
 
 
 def massif_name_to_ordered_return_level_uncertainties(model_class, last_year_for_the_data, altitudes, massif_names,
-                                                      uncertainty_methods):
+                                                      uncertainty_methods, temporal_covariate):
     # Load model name
     model_name = get_model_name(model_class)
     # Load altitude visualizer
@@ -42,7 +42,8 @@ def massif_name_to_ordered_return_level_uncertainties(model_class, last_year_for
         print('{} processing altitude = {} '.format(model_name, altitude))
         for ci_method in uncertainty_methods:
             d = visualizer.massif_name_to_altitude_and_eurocode_level_uncertainty(model_class, last_year_for_the_data,
-                                                                                  massif_names, ci_method)
+                                                                                  massif_names, ci_method,
+                                                                                  temporal_covariate)
             # Append the altitude one by one
             for massif_name, return_level_uncertainty in d.items():
                 massif_name_to_ordered_eurocode_level_uncertainty[massif_name][ci_method].append(
@@ -52,6 +53,7 @@ def massif_name_to_ordered_return_level_uncertainties(model_class, last_year_for
 
 def main_drawing():
     fast_plot = [True, False][0]
+    temporal_covariate = 2017
     # Select parameters
     massif_names = MASSIF_NAMES_ALPS[:]
     model_class_and_last_year = [
@@ -76,9 +78,9 @@ def main_drawing():
         start = time.time()
         model_name_to_massif_name_to_ordered_return_level.update(
             massif_name_to_ordered_return_level_uncertainties(model_class, last_year_for_the_data, altitudes,
-                                                              massif_names, uncertainty_methods))
+                                                              massif_names, uncertainty_methods, temporal_covariate))
         duration = time.time() - start
-        print(model_class, duration)
+        print('Duration:', model_class, duration)
     # Transform the dictionary into the desired format
     massif_name_to_model_name_to_ordered_return_level_uncertainties = {}
     for massif_name in massif_names:

experiment/meteo_france_data/scm_models_data/visualization/study_visualization/study_visualizer.py

@@ -354,9 +354,9 @@ class StudyVisualizer(VisualizationParameters):
         start_year, stop_year = self.study.start_year_and_stop_year
         return list(range(start_year, stop_year))
 
-    def massif_name_to_altitude_and_eurocode_level_uncertainty(self, model_class, last_year_for_the_data, massif_names, ci_method) -> Dict[str, Tuple[int, EurocodeConfidenceIntervalFromExtremes]]:
+    def massif_name_to_altitude_and_eurocode_level_uncertainty(self, model_class, last_year_for_the_data, massif_names, ci_method, temporal_covariate) -> Dict[str, Tuple[int, EurocodeConfidenceIntervalFromExtremes]]:
         massif_ids_and_names = [(massif_id, massif_name) for massif_id, massif_name in enumerate(self.study.study_massif_names) if massif_name in  massif_names]
-        arguments = [[last_year_for_the_data, self.smooth_maxima_x_y(massif_id), model_class, ci_method] for massif_id, _ in massif_ids_and_names]
+        arguments = [[last_year_for_the_data, self.smooth_maxima_x_y(massif_id), model_class, ci_method, temporal_covariate] for massif_id, _ in massif_ids_and_names]
         if self.multiprocessing:
             with Pool(NB_CORES) as p:
                 res = p.starmap(compute_eurocode_confidence_interval, arguments)

extreme_fit/distribution/gev/main_fevd_bayesian.R

@@ -9,9 +9,10 @@ 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)
+start_loc = 0; start_scale = 1; start_shape = 1
 
 # fevdPriorMy <- function (theta, q, p, log = FALSE){
 #     x = theta["shape"] + 0.5
@@ -43,11 +44,25 @@ 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)
+print(coord)
 colnames(coord) = c("T")
 coord = data.frame(coord, stringsAsFactors = TRUE)
-# res = fevd_fixed(x_gev, data=coord, location.fun= ~T, scale.fun= ~T, method='Bayesian', priorFun="fevdPriorCustom", priorParams=list(q=c(6), p=c(9)), iter=5000, verbose=TRUE, use.phi=FALSE)
-res = fevd_fixed(x_gev, data=coord, method='Bayesian', priorFun="fevdPriorCustom", priorParams=list(q=c(6), p=c(9)), iter=5000, verbose=TRUE, use.phi=FALSE)
+res = fevd_fixed(x_gev, data=coord, location.fun= ~T, method='Bayesian', priorFun="fevdPriorCustom", priorParams=list(q=c(6), p=c(9)), iter=5000, verbose=TRUE, use.phi=FALSE)
+# res = fevd_fixed(x_gev, data=coord, method='Bayesian', priorFun="fevdPriorCustom", priorParams=list(q=c(6), p=c(9)), iter=5000, verbose=TRUE, use.phi=FALSE)
+print(res)
+
+v = make.qcov(res, vals = list(mu1 = c(0.0)))
+# # res_find = findpars(res, FUN = "mean" , tscale = FALSE, qcov = v,
+# #     burn.in=4998)
+# # print(res_find)
+# #
+res_ci = ci(res, alpha = 0.05, type = c("return.level"),
+    burn.in=499,
+    return.period = 50, FUN = "mean", tscale = FALSE,  qcov=v)
+print(res_ci)
+print(res_ci[1, ])
+
 
 # Some display for the results
 # print(res)
@@ -67,9 +82,9 @@ res = fevd_fixed(x_gev, data=coord, method='Bayesian', priorFun="fevdPriorCustom
 
 
 # Confidence interval
-res_ci = ci(res, alpha = 0.05, type = c("return.level"),
-    return.period = 50, FUN = "mean", tscale = FALSE)
-print(res_ci)
+# res_ci = ci(res, alpha = 0.05, type = c("return.level"),
+#     return.period = 50, FUN = "mean", tscale = FALSE)
+# print(res_ci)
 # print(attributes(res_ci))
 # print(dimnames(res_ci))
 # print(res_ci[1])
@@ -77,15 +92,16 @@ print(res_ci)
 
 # covariates = c(1.0, 100.0, 1000.0)
 # v = make.qcov(res, vals = list(mu1 = covariates, sigma1 = covariates))
-# v = make.qcov(res, vals = list(mu1 = c(0.0), sigma1 = c(0.0)))
-# res_find = findpars(res, FUN = "mean" , tscale = FALSE, qcov = v,
-#     burn.in=4998)
-# print(res_find)
-#
+# v = make.qcov(res, vals = list(mu1 = c(0.0)))
+# # res_find = findpars(res, FUN = "mean" , tscale = FALSE, qcov = v,
+# #     burn.in=4998)
+# # print(res_find)
+# #
 # res_ci = ci(res, alpha = 0.05, type = c("return.level"),
 #     burn.in=4998,
 #     return.period = 50, FUN = "mean", tscale = FALSE,  qcov=v)
 # print(res_ci)
+# print(res_ci[1, ])
 #
 # r = qgev(0.98, 0.1427229 , 1.120554, -0.1008511)
 # print(r)

extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_extract_eurocode_return_level.py

@@ -16,16 +16,12 @@ from extreme_fit.model.result_from_model_fit.result_from_extremes.result_from_ml
 class AbstractExtractEurocodeReturnLevel(object):
     ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY = 0.05
 
-    def __init__(self, estimator: LinearMarginEstimator,
-                 ci_method,
-                 year_of_interest: int = YEAR_OF_INTEREST_FOR_RETURN_LEVEL,
-                 ):
+    def __init__(self, estimator: LinearMarginEstimator, ci_method, temporal_covariate):
         self.ci_method = ci_method
         self.estimator = estimator
         self.result_from_fit = self.estimator.result_from_model_fit
-
-        self.year_of_interest = year_of_interest
-
+        self.temporal_covariate = temporal_covariate
+        # Fixed Parameters
         self.eurocode_quantile = EUROCODE_QUANTILE
         self.alpha_for_confidence_interval = self.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY
 
@@ -41,15 +37,15 @@ class AbstractExtractEurocodeReturnLevel(object):
 class ExtractEurocodeReturnLevelFromCiMethod(AbstractExtractEurocodeReturnLevel):
     result_from_fit: ResultFromMleExtremes
 
-    def __init__(self, estimator: LinearMarginEstimator, ci_method,
-                 year_of_interest: int = YEAR_OF_INTEREST_FOR_RETURN_LEVEL):
-        super().__init__(estimator, ci_method, year_of_interest)
+    @property
+    def transformed_temporal_covariate(self):
+        return self.estimator.dataset.coordinates.transformation.transform_float(self.temporal_covariate)
 
     @cached_property
     def confidence_interval_method(self):
         return self.result_from_fit.confidence_interval_method(self.eurocode_quantile,
                                                                self.alpha_for_confidence_interval,
-                                                               self.year_of_interest)
+                                                               self.transformed_temporal_covariate)
 
     @property
     def mean_estimate(self) -> np.ndarray:
@@ -63,9 +59,8 @@ class ExtractEurocodeReturnLevelFromCiMethod(AbstractExtractEurocodeReturnLevel)
 class ExtractEurocodeReturnLevelFromMyBayesianExtremes(AbstractExtractEurocodeReturnLevel):
     result_from_fit: ResultFromBayesianExtremes
 
-    def __init__(self, estimator: LinearMarginEstimator, ci_method,
-                 year_of_interest: int = YEAR_OF_INTEREST_FOR_RETURN_LEVEL):
-        super().__init__(estimator, ci_method, year_of_interest)
+    def __init__(self, estimator: LinearMarginEstimator, ci_method, temporal_covariate):
+        super().__init__(estimator, ci_method, temporal_covariate)
         assert isinstance(self.result_from_fit, ResultFromBayesianExtremes)
 
     @property
@@ -78,24 +73,24 @@ class ExtractEurocodeReturnLevelFromMyBayesianExtremes(AbstractExtractEurocodeRe
         return margin_functions
 
     @property
-    def gev_params_from_fit_for_year_of_interest(self) -> List[GevParams]:
-        return [margin_function.get_gev_params(coordinate=np.array([self.year_of_interest]), is_transformed=False)
+    def gev_params_from_fit_for_temporal_covariate(self) -> List[GevParams]:
+        return [margin_function.get_gev_params(coordinate=np.array([self.temporal_covariate]), is_transformed=False)
                 for margin_function in self.margin_functions_from_fit]
 
     @cached_property
-    def posterior_eurocode_return_level_samples_for_year_of_interest(self) -> np.ndarray:
+    def posterior_eurocode_return_level_samples_for_temporal_covariate(self) -> np.ndarray:
         return np.array(
-            [p.quantile(self.eurocode_quantile) for p in self.gev_params_from_fit_for_year_of_interest])
+            [p.quantile(self.eurocode_quantile) for p in self.gev_params_from_fit_for_temporal_covariate])
 
     @property
     def mean_estimate(self) -> np.ndarray:
         # Mean posterior value here
-        return np.mean(self.posterior_eurocode_return_level_samples_for_year_of_interest)
+        return np.mean(self.posterior_eurocode_return_level_samples_for_temporal_covariate)
 
     @property
     def confidence_interval(self):
         # Bottom and upper quantile correspond to the quantile
         bottom_quantile = self.alpha_for_confidence_interval / 2
         bottom_and_upper_quantile = (bottom_quantile, 1 - bottom_quantile)
-        return [np.quantile(self.posterior_eurocode_return_level_samples_for_year_of_interest, q=q)
+        return [np.quantile(self.posterior_eurocode_return_level_samples_for_temporal_covariate, q=q)
                 for q in bottom_and_upper_quantile]

extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_result_from_extremes.py

@@ -14,11 +14,15 @@ class AbstractResultFromExtremes(AbstractResultFromModelFit):
         super().__init__(result_from_fit)
         self.gev_param_name_to_dim = gev_param_name_to_dim
 
+    @property
+    def is_non_stationary(self):
+        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, temporal_covariate):
+    def confidence_interval_method(self, quantile_level, alpha_interval, transformed_temporal_covariate):
         return_period = round(1 / (1 - quantile_level))
         common_kwargs = {
             'return.period': return_period,
@@ -27,7 +31,7 @@ class AbstractResultFromExtremes(AbstractResultFromModelFit):
             'type': r.c("return.level")
         }
         if self.gev_param_name_to_dim:
-            d = {GevParams.greek_letter_from_gev_param_name(gev_param_name) + '1': r.c(temporal_covariate) for
+            d = {GevParams.greek_letter_from_gev_param_name(gev_param_name) + '1': r.c(transformed_temporal_covariate) for
                  gev_param_name in self.gev_param_name_to_dim.keys()}
             kwargs = {
                 "vals": r.list(**d

extreme_fit/model/result_from_model_fit/result_from_extremes/eurocode_return_level_uncertainties.py

@@ -11,15 +11,14 @@ from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_int
     ConfidenceIntervalMethodFromExtremes
 
 
-def compute_eurocode_confidence_interval(last_year_for_the_data, smooth_maxima_x_y, model_class, ci_method):
+def compute_eurocode_confidence_interval(last_year_for_the_data, smooth_maxima_x_y, model_class, ci_method, temporal_covariate):
     years, smooth_maxima = smooth_maxima_x_y
     idx = years.index(last_year_for_the_data) + 1
     years, smooth_maxima = years[:idx], smooth_maxima[:idx]
-    return EurocodeConfidenceIntervalFromExtremes.from_maxima_years_model_class(smooth_maxima, years, model_class, ci_method)
+    return EurocodeConfidenceIntervalFromExtremes.from_maxima_years_model_class(smooth_maxima, years, model_class, temporal_covariate, ci_method)
 
 
 class EurocodeConfidenceIntervalFromExtremes(object):
-    YEAR_OF_INTEREST = 2017
 
     def __init__(self, mean_estimate, confidence_interval):
         self.mean_estimate = mean_estimate
@@ -31,15 +30,17 @@ class EurocodeConfidenceIntervalFromExtremes(object):
 
     @classmethod
     def from_estimator_extremes(cls, estimator_extremes: LinearMarginEstimator,
-                                ci_method: ConfidenceIntervalMethodFromExtremes):
+                                ci_method: ConfidenceIntervalMethodFromExtremes,
+                                temporal_covariate: int):
         if ci_method == ConfidenceIntervalMethodFromExtremes.my_bayes:
-            extractor = ExtractEurocodeReturnLevelFromMyBayesianExtremes(estimator_extremes, ci_method, cls.YEAR_OF_INTEREST)
+            extractor = ExtractEurocodeReturnLevelFromMyBayesianExtremes(estimator_extremes, ci_method, temporal_covariate)
         else:
-            extractor = ExtractEurocodeReturnLevelFromCiMethod(estimator_extremes, ci_method, cls.YEAR_OF_INTEREST)
+            extractor = ExtractEurocodeReturnLevelFromCiMethod(estimator_extremes, ci_method, temporal_covariate)
         return cls(extractor.mean_estimate,  extractor.confidence_interval)
 
     @classmethod
     def from_maxima_years_model_class(cls, maxima, years, model_class,
+                                      temporal_covariate,
                                       ci_method=ConfidenceIntervalMethodFromExtremes.ci_bayes):
         # Load coordinates and dataset
         coordinates, dataset = load_temporal_coordinates_and_dataset(maxima, years)
@@ -53,7 +54,7 @@ class EurocodeConfidenceIntervalFromExtremes(object):
         fitted_estimator = fitted_linear_margin_estimator(model_class, coordinates, dataset, starting_year=1958,
                                                           fit_method=fit_method)
         # Load object from result from extremes
-        return cls.from_estimator_extremes(fitted_estimator, ci_method)
+        return cls.from_estimator_extremes(fitted_estimator, ci_method, temporal_covariate)
 
 
 

extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_bayesian_extremes.py

@@ -53,9 +53,13 @@ class ResultFromBayesianExtremes(AbstractResultFromExtremes):
                 'FUN': "mean",
         }
         res = r.ci(self.result_from_fit, **bayesian_ci_parameters, **common_kwargs)
-        d = self.get_python_dictionary(res)
-        keys = ['Posterior Mean 50-year level', '95% lower CI', '95% upper CI']
-        mean_estimate, lower, upper = [np.array(d[k])[0] for k in keys]
+        if self.gev_param_name_to_dim:
+            a = np.array(res)[0]
+            lower, mean_estimate, upper = a
+        else:
+            d = self.get_python_dictionary(res)
+            keys = ['Posterior Mean 50-year 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)
 
 

spatio_temporal_dataset/coordinates/transformed_coordinates/transformation/abstract_transformation.py

@@ -13,6 +13,9 @@ class AbstractTransformation(object):
     def nb_dimensions(self):
         return self.df_coordinates.shape[1]
 
+    def transform_float(self, coordinate: float):
+        return self.transform_array(np.array([coordinate]))[0]
+
     def transform_array(self, coordinate: np.ndarray):
         assert len(coordinate) == self.nb_dimensions, "coordinate={}, nb_dimensions={}".format(coordinate,
                                                                                                self.nb_dimensions)

test/test_extreme_fit/test_model/test_confidence_interval.py

@@ -16,6 +16,8 @@ from extreme_fit.model.utils import r, set_seed_r, set_seed_for_test
 from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
 from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_coordinates import \
     AbstractTemporalCoordinates
+from spatio_temporal_dataset.coordinates.transformed_coordinates.transformation.abstract_transformation import \
+    CenteredScaledNormalization
 from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
 from spatio_temporal_dataset.spatio_temporal_observations.abstract_spatio_temporal_observations import \
     AbstractSpatioTemporalObservations
@@ -23,7 +25,7 @@ from spatio_temporal_dataset.spatio_temporal_observations.abstract_spatio_tempor
 
 class TestConfidenceInterval(unittest.TestCase):
 
-    def setUp(self) -> None:
+    def load_data(self) -> None:
         set_seed_for_test()
         r("""
         N <- 50
@@ -34,42 +36,69 @@ class TestConfidenceInterval(unittest.TestCase):
         # Compute the stationary temporal margin with isMev
         self.start_year = 0
         df = pd.DataFrame({AbstractCoordinates.COORDINATE_T: range(self.start_year, self.start_year + 50)})
-        self.coordinates = AbstractTemporalCoordinates.from_df(df)
+        self.coordinates = self.load_coordinates(df)
         df2 = pd.DataFrame(data=np.array(r['x_gev']), index=df.index)
         observations = AbstractSpatioTemporalObservations(df_maxima_gev=df2)
         self.dataset = AbstractDataset(observations=observations, coordinates=self.coordinates)
         self.model_classes = [StationaryTemporalModel]
 
+    @staticmethod
+    def load_coordinates(df):
+        return AbstractTemporalCoordinates.from_df(df)
+
     def compute_eurocode_ci(self, model_class):
+        self.load_data()
         estimator = fitted_linear_margin_estimator(model_class, self.coordinates, self.dataset,
                                                    starting_year=0,
                                                    fit_method=self.fit_method)
-        return EurocodeConfidenceIntervalFromExtremes.from_estimator_extremes(estimator, self.ci_method)
+        return EurocodeConfidenceIntervalFromExtremes.from_estimator_extremes(estimator, self.ci_method,
+                                                                              self.start_year)
+
+    @property
+    def bayesian_ci(self):
+        return {
+            StationaryTemporalModel: (6.756903450587758, 10.316338515219085, 15.77861914935531),
+            NonStationaryLocationTemporalModel: (6.588570126641043, 10.055847177064836, 14.332882862817332),
+            NonStationaryLocationAndScaleTemporalModel: (7.836837972383451, 11.162663922795906, 16.171701445841183),
+        }
 
     def test_my_bayes(self):
         self.fit_method = TemporalMarginFitMethod.extremes_fevd_bayesian
         self.ci_method = ConfidenceIntervalMethodFromExtremes.my_bayes
-        self.model_class_to_triplet = {
-            StationaryTemporalModel: (6.756903450587758, 10.316338515219085, 15.77861914935531),
-            NonStationaryLocationTemporalModel: (6.047033481540427, 9.708540966532225, 14.74058551926604),
-            NonStationaryLocationAndScaleTemporalModel: (6.383536224810785, 9.69120774797993, 13.917914357321615),
-        }
+        self.model_class_to_triplet = self.bayesian_ci
 
     def test_ci_bayes(self):
         self.fit_method = TemporalMarginFitMethod.extremes_fevd_bayesian
         self.ci_method = ConfidenceIntervalMethodFromExtremes.ci_bayes
-        self.model_class_to_triplet = {
-            StationaryTemporalModel: (6.756903450587758, 10.316338515219085, 15.77861914935531),
-            # NonStationaryLocationTemporalModel: (6.047033481540427, 9.708540966532225, 14.74058551926604),
-            # NonStationaryLocationAndScaleTemporalModel: (6.383536224810785, 9.69120774797993, 13.917914357321615),
-        }
+        self.model_class_to_triplet = self.bayesian_ci
 
     def tearDown(self) -> None:
         for model_class, expected_triplet in self.model_class_to_triplet.items():
-            eurocode_ci = self.compute_eurocode_ci(StationaryTemporalModel)
+            eurocode_ci = self.compute_eurocode_ci(model_class)
             found_triplet = eurocode_ci.triplet
             for a, b in zip(expected_triplet, found_triplet):
-                self.assertAlmostEqual(a, b, msg="{} {}".format(model_class, found_triplet))
+                self.assertAlmostEqual(a, b, msg="{} \n{}".format(model_class, found_triplet))
+
+
+class TestConfidenceIntervalModifiedCoordinates(TestConfidenceInterval):
+
+    @staticmethod
+    def load_coordinates(df):
+        return AbstractTemporalCoordinates.from_df(df, transformation_class=CenteredScaledNormalization)
+
+    @property
+    def bayesian_ci(self):
+        return {
+            StationaryTemporalModel: (6.756903450587758, 10.316338515219085, 15.77861914935531),
+            NonStationaryLocationTemporalModel: (6.266027110993808, 10.063368195790687, 14.894103640762097),
+            NonStationaryLocationAndScaleTemporalModel: (5.554116722722492, 13.714431163455535, 26.929963957448642),
+        }
+
+    def test_my_bayes(self):
+        super().test_my_bayes()
+
+    def test_ci_bayes(self):
+        super().test_ci_bayes()
 
 
 if __name__ == '__main__':