diff --git a/extreme_fit/function/param_function/param_function.py b/extreme_fit/function/param_function/param_function.py
index 0f236d27d72a10f10986af2f0184bab7651bd7ac..9b8f5aafec4a01a31c8fdc2763da25572e1734fc 100644
--- a/extreme_fit/function/param_function/param_function.py
+++ b/extreme_fit/function/param_function/param_function.py
@@ -89,6 +89,6 @@ class SplineParamFunction(AbstractParamFunction):
for i, (dim, max_degree) in enumerate(self.dim_and_degree):
assert isinstance(dim, int)
spline_coef = self.coef.dim_to_spline_coef[dim]
- spl = BSpline(t=spline_coef.knots, c=spline_coef.coefficients, k=max_degree)
- gev_param_value += spl(coordinate)[0]
+ spl = BSpline(t=spline_coef.knots, c=spline_coef.coefficients, k=max_degree, extrapolate=False)
+ gev_param_value += spl(coordinate[dim])
return gev_param_value
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 df1542dc6de41e8d9aeb7e13b86e859494bc316b..a2d54e4972eb7c0ba59ed4da75f5a6e2b4ed1114 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
@@ -15,9 +15,9 @@ from extreme_fit.model.utils import r
class AbstractResultFromExtremes(AbstractResultFromModelFit):
- def __init__(self, result_from_fit: robjects.ListVector, param_name_to_dim=None, dim_to_coordinate=None) -> None:
+ def __init__(self, result_from_fit: robjects.ListVector, param_name_to_dims=None, dim_to_coordinate=None) -> None:
super().__init__(result_from_fit)
- self.param_name_to_dim = param_name_to_dim
+ self.param_name_to_dims = param_name_to_dims
self.dim_to_coordinate = dim_to_coordinate
@property
@@ -46,7 +46,7 @@ class AbstractResultFromExtremes(AbstractResultFromModelFit):
@property
def is_non_stationary(self):
- return len(self.param_name_to_dim) > 0
+ return len(self.param_name_to_dims) > 0
def load_dataframe_from_r_matrix(self, name):
r_matrix = self.name_to_value[name]
@@ -61,21 +61,21 @@ class AbstractResultFromExtremes(AbstractResultFromModelFit):
'tscale': False,
'type': r.c("return.level") if return_level_interval else r.c("parameter")
}
- if self.param_name_to_dim:
+ if self.param_name_to_dims:
if isinstance(transformed_temporal_covariate, (int, float, np.int, np.float, np.int64)):
d = {GevParams.greek_letter_from_param_name_confidence_interval(param_name) + '1': r.c(transformed_temporal_covariate) for
- param_name in self.param_name_to_dim.keys()}
+ param_name in self.param_name_to_dims.keys()}
elif isinstance(transformed_temporal_covariate, np.ndarray):
d = OrderedDict()
- linearity_in_shape = GevParams.SHAPE in self.param_name_to_dim
+ linearity_in_shape = GevParams.SHAPE in self.param_name_to_dims
nb_calls = 4 # or 4 (1 and 3 did not work for the test)
for param_name in GevParams.PARAM_NAMES:
- suffix = '0' if param_name in self.param_name_to_dim else ''
+ suffix = '0' if param_name in self.param_name_to_dims else ''
covariate = np.array([1] * nb_calls)
d2 = {GevParams.greek_letter_from_param_name_confidence_interval(param_name, linearity_in_shape) + suffix: r.c(covariate)}
d.update(d2)
- if param_name in self.param_name_to_dim:
- for coordinate_idx, _ in self.param_name_to_dim[param_name]:
+ if param_name in self.param_name_to_dims:
+ for coordinate_idx, _ in self.param_name_to_dims[param_name]:
idx_str = str(coordinate_idx + 1)
covariate = float(transformed_temporal_covariate.copy()[coordinate_idx])
covariate = np.array([covariate] * nb_calls)
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 e2176c23875bd15e5ab7d1c74da8bb1bb80195da..4ebd71101100ee9111f588e67678ad873b00caa2 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
@@ -51,7 +51,7 @@ class ResultFromBayesianExtremes(AbstractResultFromExtremes):
def get_coef_dict_from_posterior_sample(self, s: pd.Series):
assert len(s) >= 3
values = {i: v for i, v in enumerate(s)}
- return get_margin_coef_ordered_dict(self.param_name_to_dim, values)
+ return get_margin_coef_ordered_dict(self.param_name_to_dims, values)
@property
def margin_coef_ordered_dict(self):
@@ -65,7 +65,7 @@ class ResultFromBayesianExtremes(AbstractResultFromExtremes):
'FUN': "mean",
}
res = r.ci(self.result_from_fit, **bayesian_ci_parameters, **common_kwargs)
- if self.param_name_to_dim:
+ if self.param_name_to_dims:
a = np.array(res)[0]
lower, mean_estimate, upper = a
else:
diff --git a/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_evgam.py b/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_evgam.py
index 592066be4c5b3cdbb77aaccceb542746d43a1580..9ff502683cd815e1cb36996f9b267664dcc17463 100644
--- a/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_evgam.py
+++ b/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_evgam.py
@@ -32,12 +32,19 @@ class ResultFromEvgam(AbstractResultFromExtremes):
assert len(p) == len(self.maxima)
# Compute nllh
nllh = 0
- for loc, log_scale, shape, maximum in zip(*parameters):
- gev_params = GevParams(loc, np.exp(log_scale), shape)
+ for j, (loc, log_scale, shape, maximum) in enumerate(zip(*parameters)):
+ scale = np.exp(log_scale)
+ gev_params = GevParams(loc, scale, shape)
p = gev_params.density(maximum)
nllh += -np.log(p)
return nllh
+ def get_gev_params(self, idx):
+ param_names = ['location', 'logscale', 'shape']
+ loc, log_scale, shape = [np.array(self.get_python_dictionary(self.name_to_value[param_name])['fitted'])[idx]
+ for param_name in param_names]
+ return GevParams(loc, np.exp(log_scale), shape)
+
@property
def maxima(self):
return np.array(self.get_python_dictionary(self.name_to_value['location'])['model'][0])
@@ -61,39 +68,55 @@ class ResultFromEvgam(AbstractResultFromExtremes):
@property
def margin_coef_ordered_dict(self):
coefficients = np.array(self.name_to_value['coefficients'])
- param_name_to_str_formula = {k: str(v) for k, v in self.get_python_dictionary(self.name_to_value['formula']).items()}
+ param_name_to_str_formula = {k: str(v) for k, v in
+ self.get_python_dictionary(self.name_to_value['formula']).items()}
r_param_names_with_spline = [k for k, v in param_name_to_str_formula.items() if "s(" in v]
+ r_param_names_with_spline = [k if k != "scale" else "logscale" for k in r_param_names_with_spline]
if len(r_param_names_with_spline) == 0:
- return get_margin_coef_ordered_dict(self.param_name_to_dim, coefficients, self.type_for_mle,
+ return get_margin_coef_ordered_dict(self.param_name_to_dims, coefficients, self.type_for_mle,
dim_to_coordinate_name=self.dim_to_coordinate)
else:
- # todo we might need to delete some coefficient for the spline param so that it does not create some assertion error
- coef_dict = get_margin_coef_ordered_dict(self.param_name_to_dim, coefficients, self.type_for_mle,
- dim_to_coordinate_name=self.dim_to_coordinate)
+ # Compute spline param_name to dim_to_knots_and_coefficients
spline_param_name_to_dim_knots_and_coefficient = {}
- for r_param_name in r_param_names_with_spline:
- print('here')
- param_name = self.r_param_name_to_param_name[r_param_name]
- dim_knots_and_coefficient = {}
- dims = self.param_name_to_dim[param_name]
- if len(dims) > 1:
- raise NotImplementedError
- else:
- dim, max_degree = dims[0]
- d = self.get_python_dictionary(self.name_to_value[r_param_name])
- coefficients = np.array(d["coefficients"])
- smooth = list(d['smooth'])[0]
- knots = np.array(self.get_python_dictionary(smooth)['knots'])
- assert len(knots) == len(coefficients) + 1 + max_degree
- dim_knots_and_coefficient[dim] = (knots, coefficients)
+ param_names_with_spline = [self.r_param_name_to_param_name[r_param_name] for r_param_name in
+ r_param_names_with_spline]
+ for param_name, r_param_name in zip(param_names_with_spline, r_param_names_with_spline):
+ dim_knots_and_coefficient = self.compute_dim_to_knots_and_coefficient(param_name, r_param_name)
spline_param_name_to_dim_knots_and_coefficient[param_name] = dim_knots_and_coefficient
-
+ # Modify the param_name_to_dim for the spline variables (to not miss any variable)
+ param_name_to_dims = self.param_name_to_dims.copy()
+ for param_name in param_names_with_spline:
+ new_dims = []
+ for dim, _ in param_name_to_dims[param_name]:
+ nb_coefficients_for_param_name = len(
+ spline_param_name_to_dim_knots_and_coefficient[param_name][dim][1])
+ new_dim = (dim, nb_coefficients_for_param_name - 1)
+ new_dims.append(new_dim)
+ param_name_to_dims[param_name] = new_dims
+ # Extract the coef list
+ coef_dict = get_margin_coef_ordered_dict(param_name_to_dims, coefficients, self.type_for_mle,
+ dim_to_coordinate_name=self.dim_to_coordinate)
return coef_dict, spline_param_name_to_dim_knots_and_coefficient
+ def compute_dim_to_knots_and_coefficient(self, param_name, r_param_name):
+ dim_knots_and_coefficient = {}
+ dims = self.param_name_to_dims[param_name]
+ if len(dims) > 1:
+ raise NotImplementedError
+ else:
+ dim, max_degree = dims[0]
+ d = self.get_python_dictionary(self.name_to_value[r_param_name])
+ coefficients = np.array(d["coefficients"])
+ smooth = list(d['smooth'])[0]
+ knots = np.array(self.get_python_dictionary(smooth)['knots'])
+ assert len(knots) == len(coefficients) + 1 + max_degree
+ dim_knots_and_coefficient[dim] = (knots, coefficients)
+ return dim_knots_and_coefficient
+
@property
def r_param_name_to_param_name(self):
return {
'location': GevParams.LOC,
- 'scale': GevParams.SCALE,
+ 'logscale': GevParams.SCALE,
'shape': GevParams.SHAPE,
- }
\ No newline at end of file
+ }
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 e24978c3413a43424e80a5f7c7e6cdd0b58af55b..dd7553c366a5fb77d2c2f98f100e2f80b91ce17b 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
@@ -25,7 +25,7 @@ class ResultFromMleExtremes(AbstractResultFromExtremes):
values = {i: param for i, param in enumerate(np.array(d['par']))}
else:
values = {i: np.array(v)[0] for i, v in enumerate(d.values())}
- return get_margin_coef_ordered_dict(self.param_name_to_dim, values, self.type_for_mle,
+ return get_margin_coef_ordered_dict(self.param_name_to_dims, values, self.type_for_mle,
dim_to_coordinate_name=self.dim_to_coordinate)
def _confidence_interval_method(self, common_kwargs, ci_method, return_period):
diff --git a/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_spline.py b/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_spline.py
index 90360e3ba5b22ab6aab16ca9347b216638ddbd10..101205677e32c0d445d4426eb5bdf05822eb119f 100644
--- a/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_spline.py
+++ b/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_spline.py
@@ -4,7 +4,8 @@ import numpy as np
import pandas as pd
from extreme_fit.distribution.gev.gev_params import GevParams
-from extreme_fit.model.margin_model.polynomial_margin_model.polynomial_margin_model import NonStationaryQuadraticLocationModel, \
+from extreme_fit.model.margin_model.polynomial_margin_model.polynomial_margin_model import \
+ NonStationaryQuadraticLocationModel, \
NonStationaryQuadraticScaleModel, NonStationaryQuadraticLocationGumbelModel, NonStationaryQuadraticScaleGumbelModel
from extreme_fit.model.margin_model.spline_margin_model.temporal_spline_model_degree_1 import \
NonStationaryTwoLinearLocationModel, NonStationaryTwoLinearScaleModel
@@ -48,8 +49,7 @@ class TestGevTemporalSpline(unittest.TestCase):
starting_year=0,
fit_method=self.fit_method)
# Checks that parameters returned are indeed different
- mle_params_estimated_year1 = estimator.function_from_fit.get_params(np.array([1])).to_dict()
- print(mle_params_estimated_year1, type(mle_params_estimated_year1))
+ mle_params_estimated_year1 = estimator.function_from_fit.get_params(np.array([0])).to_dict()
mle_params_estimated_year3 = estimator.function_from_fit.get_params(np.array([21])).to_dict()
mle_params_estimated_year5 = estimator.function_from_fit.get_params(np.array([self.last_year])).to_dict()
self.assertNotEqual(mle_params_estimated_year1, mle_params_estimated_year3)
@@ -58,18 +58,26 @@ class TestGevTemporalSpline(unittest.TestCase):
diff1 = mle_params_estimated_year1[param_to_test] - mle_params_estimated_year3[param_to_test]
diff2 = mle_params_estimated_year3[param_to_test] - mle_params_estimated_year5[param_to_test]
self.assertNotAlmostEqual(diff1, diff2)
+ # Assert that the computed parameters are the same
+ # for year in range(self.start_year, self.start_year+5):
+ # gev_params_from_result = estimator.result_from_model_fit.get_gev_params(year).to_dict()
+ # my_gev_params = estimator.function_from_fit.get_params(np.array([year])).to_dict()
+ # for param_name in GevParams.PARAM_NAMES:
+ # self.assertAlmostEqual(gev_params_from_result[param_name], my_gev_params[param_name],
+ # msg='for the {} parameter at year={}'.format(param_name, year))
# Assert that indicators are correctly computed
- self.assertAlmostEqual(estimator.result_from_model_fit.nllh, estimator.nllh())
+ self.assertAlmostEqual(estimator.result_from_model_fit.nllh, estimator.nllh(), places=0)
# self.assertAlmostEqual(estimator.result_from_model_fit.aic, estimator.aic())
# self.assertAlmostEqual(estimator.result_from_model_fit.bic, estimator.bic())
- # def test_gev_temporal_margin_fit_spline_two_linear_location(self):
- # self.function_test_gev_temporal_margin_fit_non_stationary_spline(NonStationaryTwoLinearLocationModel,
- # param_to_test=GevParams.LOC)
+ def test_gev_temporal_margin_fit_spline_two_linear_location(self):
+ self.function_test_gev_temporal_margin_fit_non_stationary_spline(NonStationaryTwoLinearLocationModel,
+ param_to_test=GevParams.LOC)
+
#
- # def test_gev_temporal_margin_fit_spline_two_linear_scale(self):
- # self.function_test_gev_temporal_margin_fit_non_stationary_spline(NonStationaryTwoLinearScaleModel,
- # param_to_test=GevParams.SCALE)
+ def test_gev_temporal_margin_fit_spline_two_linear_scale(self):
+ self.function_test_gev_temporal_margin_fit_non_stationary_spline(NonStationaryTwoLinearScaleModel,
+ param_to_test=GevParams.SCALE)
if __name__ == '__main__':