Commit 7596a6dd authored by Le Roux Erwan's avatar Le Roux Erwan
Browse files

[EXTREME ESTIMATOR] Add MLE extreme estimator. refactor test gev temporal...

[EXTREME ESTIMATOR] Add MLE extreme estimator. refactor test gev temporal extremes. move confidence interval results
parent 37558375
No related merge requests found
Showing with 298 additions and 60 deletions
+298 -60
......@@ -2,9 +2,8 @@ from typing import Dict, List
import matplotlib.pyplot as plt
from experiment.eurocode_data.eurocode_return_level_uncertainties import EurocodeLevelUncertaintyFromExtremes
from experiment.eurocode_data.massif_name_to_departement import DEPARTEMENT_TYPES, massif_name_to_eurocode_region
from experiment.eurocode_data.utils import EUROCODE_QUANTILE, EUROCODE_ALTITUDES
from extreme_fit.model.result_from_model_fit.result_from_extremes.eurocode_return_level_uncertainties import EurocodeConfidenceIntervalFromExtremes
from experiment.eurocode_data.massif_name_to_departement import massif_name_to_eurocode_region
from experiment.meteo_france_data.scm_models_data.visualization.utils import create_adjusted_axes
from root_utils import get_display_name_from_object_type
......@@ -55,7 +54,7 @@ def plot_model_name_to_uncertainty_method_to_ordered_dict(d, massif_name, axes):
def plot_label_to_ordered_return_level_uncertainties(ax, massif_name, model_name,
label_to_ordered_return_level_uncertainties:
Dict[str, List[
EurocodeLevelUncertaintyFromExtremes]]):
EurocodeConfidenceIntervalFromExtremes]]):
""" Generic function that might be used by many other more global functions"""
colors = ['tab:blue', 'tab:orange', 'tab:purple', 'tab:olive']
alpha = 0.2
......
import time
from collections import OrderedDict
from experiment.eurocode_data.eurocode_return_level_uncertainties import ConfidenceIntervalMethodFromExtremes
from extreme_fit.model.result_from_model_fit.result_from_extremes.eurocode_return_level_uncertainties import \
ConfidenceIntervalMethodFromExtremes
from experiment.eurocode_data.eurocode_visualizer import \
plot_massif_name_to_model_name_to_uncertainty_method_to_ordered_dict, get_model_name
from experiment.eurocode_data.massif_name_to_departement import DEPARTEMENT_TYPES, MASSIF_NAMES_ALPS
from experiment.eurocode_data.massif_name_to_departement import MASSIF_NAMES_ALPS
from experiment.eurocode_data.utils import EUROCODE_ALTITUDES, LAST_YEAR_FOR_EUROCODE
from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSwe3Days
from experiment.meteo_france_data.scm_models_data.visualization.hypercube_visualization.altitude_hypercube_visualizer import \
......@@ -13,7 +14,6 @@ from experiment.meteo_france_data.scm_models_data.visualization.hypercube_visual
load_altitude_visualizer
from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import StationaryTemporalModel, \
NonStationaryLocationAndScaleTemporalModel
from root_utils import get_display_name_from_object_type
# Model class
import matplotlib as mpl
......@@ -22,7 +22,8 @@ mpl.rcParams['text.usetex'] = True
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):
def massif_name_to_ordered_return_level_uncertainties(model_class, last_year_for_the_data, altitudes, massif_names,
uncertainty_methods):
# Load model name
model_name = get_model_name(model_class)
# Load altitude visualizer
......@@ -35,20 +36,20 @@ def massif_name_to_ordered_return_level_uncertainties(model_class, last_year_for
trend_test_class=None) # type: AltitudeHypercubeVisualizer
# Loop on the data
assert isinstance(altitude_visualizer.tuple_to_study_visualizer, OrderedDict)
massif_name_to_ordered_eurocode_level_uncertainty = {massif_name: {ci_method: [] for ci_method in uncertainty_methods} for massif_name in massif_names}
massif_name_to_ordered_eurocode_level_uncertainty = {
massif_name: {ci_method: [] for ci_method in uncertainty_methods} for massif_name in massif_names}
for altitude, visualizer in altitude_visualizer.tuple_to_study_visualizer.items():
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)
d = visualizer.massif_name_to_altitude_and_eurocode_level_uncertainty(model_class, last_year_for_the_data,
massif_names, ci_method)
# 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(return_level_uncertainty)
massif_name_to_ordered_eurocode_level_uncertainty[massif_name][ci_method].append(
return_level_uncertainty)
return {model_name: massif_name_to_ordered_eurocode_level_uncertainty}
def main_drawing():
fast_plot = [True, False][0]
# Select parameters
......@@ -60,8 +61,10 @@ def main_drawing():
][1:]
altitudes = EUROCODE_ALTITUDES[:]
uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.my_bayes, ConfidenceIntervalMethodFromExtremes.bayes]
show = True
if fast_plot:
show = True
model_class_and_last_year = model_class_and_last_year[:1]
altitudes = altitudes[2:4]
# altitudes = altitudes[:]
......@@ -72,7 +75,8 @@ def main_drawing():
for model_class, last_year_for_the_data in model_class_and_last_year:
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_name_to_ordered_return_level_uncertainties(model_class, last_year_for_the_data, altitudes,
massif_names, uncertainty_methods))
duration = time.time() - start
print(model_class, duration)
# Transform the dictionary into the desired format
......@@ -84,7 +88,7 @@ def main_drawing():
# Plot graph
plot_massif_name_to_model_name_to_uncertainty_method_to_ordered_dict(
massif_name_to_model_name_to_ordered_return_level_uncertainties, nb_massif_names=len(massif_names),
nb_model_names=len(model_class_and_last_year), show=True)
nb_model_names=len(model_class_and_last_year), show=show)
if __name__ == '__main__':
......
......@@ -11,9 +11,8 @@ import numpy as np
import pandas as pd
import seaborn as sns
from experiment.eurocode_data.eurocode_return_level_uncertainties import ExtractEurocodeReturnLevelFromExtremes, \
EurocodeLevelUncertaintyFromExtremes, compute_eurocode_level_uncertainty
from experiment.eurocode_data.massif_name_to_departement import dep_class_to_massif_names
from extreme_fit.model.result_from_model_fit.result_from_extremes.eurocode_return_level_uncertainties import \
EurocodeConfidenceIntervalFromExtremes, compute_eurocode_confidence_interval
from experiment.meteo_france_data.scm_models_data.abstract_extended_study import AbstractExtendedStudy
from experiment.trend_analysis.abstract_score import MeanScore, AbstractTrendScore
from experiment.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
......@@ -355,14 +354,14 @@ 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, EurocodeLevelUncertaintyFromExtremes]]:
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]]:
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]
if self.multiprocessing:
with Pool(NB_CORES) as p:
res = p.starmap(compute_eurocode_level_uncertainty, arguments)
res = p.starmap(compute_eurocode_confidence_interval, arguments)
else:
res = [compute_eurocode_level_uncertainty(*argument) for argument in arguments]
res = [compute_eurocode_confidence_interval(*argument) for argument in arguments]
res_and_altitude = [(self.study.altitude, r) for r in res]
massif_name_to_eurocode_return_level_uncertainty = OrderedDict(zip([massif_name for _, massif_name in massif_ids_and_names], res_and_altitude))
return massif_name_to_eurocode_return_level_uncertainty
......
......@@ -6,9 +6,10 @@ import pandas as pd
from extreme_fit.distribution.gev.gev_params import GevParams
from extreme_fit.model.margin_model.linear_margin_model.linear_margin_model import LinearMarginModel
from extreme_fit.model.result_from_model_fit.abstract_result_from_model_fit import AbstractResultFromModelFit
from extreme_fit.model.result_from_model_fit.result_from_extremes import ResultFromExtremes, ResultFromBayesianExtremes
from extreme_fit.model.result_from_model_fit.result_from_extremes.result_from_bayesian_extremes import AbstractResultFromExtremes, ResultFromBayesianExtremes
from extreme_fit.model.result_from_model_fit.result_from_extremes.result_from_mle_extremes import ResultFromMleExtremes
from extreme_fit.model.result_from_model_fit.result_from_ismev import ResultFromIsmev
from extreme_fit.model.utils import r, ro, get_null, get_margin_formula_extremes, get_coord, get_coord_df
from extreme_fit.model.utils import r, ro, get_null, get_margin_formula_extremes, get_coord_df
from extreme_fit.model.utils import safe_run_r_estimator
from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
......@@ -21,8 +22,6 @@ class TemporalMarginFitMethod(Enum):
class AbstractTemporalLinearMarginModel(LinearMarginModel):
"""Linearity only with respect to the temporal coordinates"""
ISMEV_GEV_FIT_METHOD_STR = 'isMev.gev.fit'
EXTREMES_FEVD_BAYESIAN_FIT_METHOD_STR = 'extRemes.fevd.Bayesian'
def __init__(self, coordinates: AbstractCoordinates, use_start_value=False, params_start_fit=None,
params_sample=None, starting_point=None, fit_method=TemporalMarginFitMethod.is_mev_gev_fit):
......@@ -52,7 +51,7 @@ class AbstractTemporalLinearMarginModel(LinearMarginModel):
# Gev fit with extRemes package
def extremes_fevd_mle_fit(self, x, df_coordinates_temp) -> ResultFromExtremes:
def extremes_fevd_mle_fit(self, x, df_coordinates_temp) -> AbstractResultFromExtremes:
r_type_argument_kwargs, y = self.extreme_arguments(df_coordinates_temp)
res = safe_run_r_estimator(function=r('fevd_fixed'),
x=x,
......@@ -60,9 +59,9 @@ class AbstractTemporalLinearMarginModel(LinearMarginModel):
method='MLE',
**r_type_argument_kwargs
)
return ResultFromExtremes(res, self.margin_function_start_fit.gev_param_name_to_dims)
return ResultFromMleExtremes(res, self.margin_function_start_fit.gev_param_name_to_dims)
def extremes_fevd_bayesian_fit(self, x, df_coordinates_temp) -> ResultFromExtremes:
def extremes_fevd_bayesian_fit(self, x, df_coordinates_temp) -> AbstractResultFromExtremes:
r_type_argument_kwargs, y = self.extreme_arguments(df_coordinates_temp)
# Assert for any non-stationary model that the shape parameter is constant
# (because the prior function considers that the last parameter should be the shape)
......
......@@ -7,10 +7,14 @@ class AbstractResultFromModelFit(object):
def __init__(self, result_from_fit: robjects.ListVector) -> None:
if hasattr(result_from_fit, 'names'):
self.name_to_value = {name: result_from_fit.rx2(name) for name in result_from_fit.names}
self.name_to_value = self.get_python_dictionary(result_from_fit)
else:
self.name_to_value = {}
@staticmethod
def get_python_dictionary(r_dictionary):
return {name: r_dictionary.rx2(name) for name in r_dictionary.names}
@property
def names(self):
return self.name_to_value.keys()
......
from enum import Enum
from typing import List
import numpy as np
from cached_property import cached_property
from experiment.eurocode_data.utils import EUROCODE_QUANTILE, YEAR_OF_INTEREST_FOR_RETURN_LEVEL
from experiment.trend_analysis.univariate_test.utils import load_temporal_coordinates_and_dataset, \
fitted_linear_margin_estimator
from extreme_fit.distribution.gev.gev_params import GevParams
from extreme_fit.estimator.abstract_estimator import AbstractEstimator
from extreme_fit.estimator.margin_estimator.abstract_margin_estimator import LinearMarginEstimator
from extreme_fit.estimator.utils import load_margin_function
from extreme_fit.model.margin_model.linear_margin_model.abstract_temporal_linear_margin_model import \
AbstractTemporalLinearMarginModel, TemporalMarginFitMethod
from extreme_fit.model.margin_model.margin_function.linear_margin_function import LinearMarginFunction
from extreme_fit.model.result_from_model_fit.result_from_extremes import ResultFromExtremes, ResultFromBayesianExtremes
class ConfidenceIntervalMethodFromExtremes(Enum):
# Confidence interval from the ci function
bayes = 0
normal = 1
boot = 2
proflik = 3
# Confidence interval from my functions
my_bayes = 4
def compute_eurocode_level_uncertainty(last_year_for_the_data, smooth_maxima_x_y, model_class, ci_method):
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 EurocodeLevelUncertaintyFromExtremes.from_maxima_years_model_class(smooth_maxima, years, model_class, ci_method)
class EurocodeLevelUncertaintyFromExtremes(object):
YEAR_OF_INTEREST = 2017
def __init__(self, posterior_mean, poster_uncertainty_interval):
self.posterior_mean = posterior_mean
self.poster_uncertainty_interval = poster_uncertainty_interval
@classmethod
def from_estimator_extremes(cls, estimator_extremes: LinearMarginEstimator,
ci_method: ConfidenceIntervalMethodFromExtremes):
extractor = ExtractEurocodeReturnLevelFromExtremes(estimator_extremes, ci_method, cls.YEAR_OF_INTEREST)
return cls(extractor.posterior_mean_eurocode_return_level_for_the_year_of_interest,
extractor.posterior_eurocode_return_level_uncertainty_interval_for_the_year_of_interest)
@classmethod
def from_maxima_years_model_class(cls, maxima, years, model_class,
ci_method=ConfidenceIntervalMethodFromExtremes.bayes):
# Load coordinates and dataset
coordinates, dataset = load_temporal_coordinates_and_dataset(maxima, years)
# Select fit method depending on the ci_method
if ci_method in [ConfidenceIntervalMethodFromExtremes.bayes,
ConfidenceIntervalMethodFromExtremes.my_bayes]:
fit_method = TemporalMarginFitMethod.extremes_fevd_bayesian
else:
fit_method = TemporalMarginFitMethod.extremes_fevd_mle
# Fitted estimator
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)
class ExtractFromExtremes(object):
pass
from extreme_fit.model.result_from_model_fit.result_from_extremes.result_from_bayesian_extremes import \
ResultFromBayesianExtremes
class ExtractEurocodeReturnLevelFromExtremes(object):
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,
alpha_for_confidence_interval: int = ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY,
):
self.ci_method = ci_method
self.estimator = estimator
self.result_from_fit = self.estimator.result_from_model_fit # type: ResultFromBayesianExtremes
assert isinstance(self.result_from_fit, ResultFromBayesianExtremes)
self.result_from_fit = self.estimator.result_from_model_fit
self.year_of_interest = year_of_interest
self.alpha_for_confidence_interval = alpha_for_confidence_interval
self.eurocode_quantile = EUROCODE_QUANTILE
self.alpha_for_confidence_interval = self.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY
class ExtractEurocodeReturnLevelFromCiMethod(AbstractExtractEurocodeReturnLevel):
pass
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)
assert isinstance(self.result_from_fit, ResultFromBayesianExtremes)
@property
def margin_functions_from_fit(self) -> List[LinearMarginFunction]:
......@@ -101,7 +58,8 @@ class ExtractEurocodeReturnLevelFromExtremes(object):
@cached_property
def posterior_eurocode_return_level_samples_for_year_of_interest(self) -> np.ndarray:
"""We divide by 100 to transform the snow water equivalent into snow load"""
return np.array([p.quantile(EUROCODE_QUANTILE) for p in self.gev_params_from_fit_for_year_of_interest]) / 100
return np.array(
[p.quantile(self.eurocode_quantile) for p in self.gev_params_from_fit_for_year_of_interest]) / 100
@property
def posterior_mean_eurocode_return_level_for_the_year_of_interest(self) -> np.ndarray:
......@@ -113,4 +71,4 @@ class ExtractEurocodeReturnLevelFromExtremes(object):
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)
for q in bottom_and_upper_quantile]
for q in bottom_and_upper_quantile]
\ No newline at end of file
import numpy as np
import pandas as pd
from rpy2 import robjects
from extreme_fit.model.result_from_model_fit.abstract_result_from_model_fit import \
AbstractResultFromModelFit
from extreme_fit.model.utils import r
class AbstractResultFromExtremes(AbstractResultFromModelFit):
def __init__(self, result_from_fit: robjects.ListVector, gev_param_name_to_dim=None) -> None:
super().__init__(result_from_fit)
self.gev_param_name_to_dim = gev_param_name_to_dim
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))
from enum import Enum
class ConfidenceIntervalMethodFromExtremes(Enum):
# Confidence interval from the ci function
bayes = 0
normal = 1
boot = 2
proflik = 3
# Confidence interval from my functions
my_bayes = 4
from enum import Enum
from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_extract_eurocode_return_level import \
ExtractEurocodeReturnLevelFromMyBayesianExtremes, ExtractEurocodeReturnLevelFromCiMethod
from experiment.trend_analysis.univariate_test.utils import load_temporal_coordinates_and_dataset, \
fitted_linear_margin_estimator
from extreme_fit.estimator.margin_estimator.abstract_margin_estimator import LinearMarginEstimator
from extreme_fit.model.margin_model.linear_margin_model.abstract_temporal_linear_margin_model import \
TemporalMarginFitMethod
from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import \
ConfidenceIntervalMethodFromExtremes
def compute_eurocode_confidence_interval(last_year_for_the_data, smooth_maxima_x_y, model_class, ci_method):
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)
class EurocodeConfidenceIntervalFromExtremes(object):
YEAR_OF_INTEREST = 2017
def __init__(self, posterior_mean, poster_uncertainty_interval):
self.posterior_mean = posterior_mean
self.poster_uncertainty_interval = poster_uncertainty_interval
@classmethod
def from_estimator_extremes(cls, estimator_extremes: LinearMarginEstimator,
ci_method: ConfidenceIntervalMethodFromExtremes):
if ci_method == ConfidenceIntervalMethodFromExtremes.my_bayes:
extractor = ExtractEurocodeReturnLevelFromMyBayesianExtremes(estimator_extremes, ci_method, cls.YEAR_OF_INTEREST)
else:
extractor = ExtractEurocodeReturnLevelFromCiMethod(estimator_extremes, ci_method, cls.YEAR_OF_INTEREST)
return cls(extractor.posterior_mean_eurocode_return_level_for_the_year_of_interest,
extractor.posterior_eurocode_return_level_uncertainty_interval_for_the_year_of_interest)
@classmethod
def from_maxima_years_model_class(cls, maxima, years, model_class,
ci_method=ConfidenceIntervalMethodFromExtremes.bayes):
# Load coordinates and dataset
coordinates, dataset = load_temporal_coordinates_and_dataset(maxima, years)
# Select fit method depending on the ci_method
if ci_method in [ConfidenceIntervalMethodFromExtremes.bayes,
ConfidenceIntervalMethodFromExtremes.my_bayes]:
fit_method = TemporalMarginFitMethod.extremes_fevd_bayesian
else:
fit_method = TemporalMarginFitMethod.extremes_fevd_mle
# Fitted estimator
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)
import numpy as np
import pandas as pd
from rpy2 import robjects
from extreme_fit.distribution.gev.gev_params import GevParams
from extreme_fit.model.result_from_model_fit.abstract_result_from_model_fit import \
AbstractResultFromModelFit
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.utils import get_margin_coef_ordered_dict
from extreme_fit.model.utils import r
class ResultFromExtremes(AbstractResultFromModelFit):
def __init__(self, result_from_fit: robjects.ListVector, gev_param_name_to_dim=None) -> None:
super().__init__(result_from_fit)
self.gev_param_name_to_dim = gev_param_name_to_dim
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))
class ResultFromBayesianExtremes(ResultFromExtremes):
class ResultFromBayesianExtremes(AbstractResultFromExtremes):
def __init__(self, result_from_fit: robjects.ListVector, gev_param_name_to_dim=None,
burn_in_percentage=0.5) -> None:
......
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.utils import get_margin_coef_ordered_dict
class ResultFromMleExtremes(AbstractResultFromExtremes):
@property
def margin_coef_ordered_dict(self):
values = self.name_to_value['results']
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)
......@@ -4,13 +4,10 @@ import numpy as np
import pandas as pd
from experiment.trend_analysis.univariate_test.abstract_gev_trend_test import fitted_linear_margin_estimator
from extreme_fit.distribution.gev.gev_params import GevParams
from extreme_fit.estimator.margin_estimator.abstract_margin_estimator import LinearMarginEstimator
from extreme_fit.model.margin_model.linear_margin_model.abstract_temporal_linear_margin_model import \
AbstractTemporalLinearMarginModel, TemporalMarginFitMethod
TemporalMarginFitMethod
from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import StationaryTemporalModel, \
NonStationaryLocationTemporalModel, NonStationaryLocationAndScaleTemporalModel
from extreme_fit.model.result_from_model_fit.result_from_extremes import ResultFromExtremes
from extreme_fit.model.utils import r, set_seed_r
from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_coordinates import \
......@@ -18,10 +15,9 @@ from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_
from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
from spatio_temporal_dataset.spatio_temporal_observations.abstract_spatio_temporal_observations import \
AbstractSpatioTemporalObservations
from test.test_utils import load_non_stationary_temporal_margin_models
class TestGevTemporalBayesian(unittest.TestCase):
class TestGevTemporalExtremesBayesian(unittest.TestCase):
def setUp(self) -> None:
set_seed_r()
......@@ -42,12 +38,12 @@ class TestGevTemporalBayesian(unittest.TestCase):
def test_gev_temporal_margin_fit_stationary(self):
# Create estimator
estimator_fitted = fitted_linear_margin_estimator(StationaryTemporalModel, self.coordinates, self.dataset,
estimator = fitted_linear_margin_estimator(StationaryTemporalModel, self.coordinates, self.dataset,
starting_year=0,
fit_method=self.fit_method)
ref = {'loc': 0.34272436381693616, 'scale': 1.3222588712831973, 'shape': 0.30491484962825105}
for year in range(1, 3):
mle_params_estimated = estimator_fitted.margin_function_from_fit.get_gev_params(np.array([year])).to_dict()
mle_params_estimated = estimator.margin_function_from_fit.get_gev_params(np.array([year])).to_dict()
for key in ref.keys():
self.assertAlmostEqual(ref[key], mle_params_estimated[key], places=3)
......
import unittest
import numpy as np
import pandas as pd
from experiment.trend_analysis.univariate_test.abstract_gev_trend_test import fitted_linear_margin_estimator
from extreme_fit.model.margin_model.linear_margin_model.abstract_temporal_linear_margin_model import \
TemporalMarginFitMethod
from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import StationaryTemporalModel, \
NonStationaryLocationTemporalModel, NonStationaryLocationAndScaleTemporalModel
from extreme_fit.model.utils import r, set_seed_r
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.dataset.abstract_dataset import AbstractDataset
from spatio_temporal_dataset.spatio_temporal_observations.abstract_spatio_temporal_observations import \
AbstractSpatioTemporalObservations
class TestGevTemporalExtremesMle(unittest.TestCase):
def setUp(self) -> None:
set_seed_r()
r("""
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
""")
# 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)
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.fit_method = TemporalMarginFitMethod.extremes_fevd_mle
def test_gev_temporal_margin_fit_stationary(self):
# Create estimator
estimator = fitted_linear_margin_estimator(StationaryTemporalModel, self.coordinates, self.dataset,
starting_year=0,
fit_method=self.fit_method)
ref = {'loc': 0.02191974259369493, 'scale': 1.0347946062900268, 'shape': 0.829052520147379}
for year in range(1, 3):
mle_params_estimated = estimator.margin_function_from_fit.get_gev_params(np.array([year])).to_dict()
for key in ref.keys():
self.assertAlmostEqual(ref[key], mle_params_estimated[key], places=3)
def test_gev_temporal_margin_fit_non_stationary_location(self):
# Create estimator
estimator = fitted_linear_margin_estimator(NonStationaryLocationTemporalModel, self.coordinates, self.dataset,
starting_year=0,
fit_method=self.fit_method)
# Checks that parameters returned are indeed different
mle_params_estimated_year1 = estimator.margin_function_from_fit.get_gev_params(np.array([1])).to_dict()
mle_params_estimated_year3 = estimator.margin_function_from_fit.get_gev_params(np.array([3])).to_dict()
self.assertNotEqual(mle_params_estimated_year1, mle_params_estimated_year3)
def test_gev_temporal_margin_fit_non_stationary_location_and_scale(self):
# Create estimator
estimator = fitted_linear_margin_estimator(NonStationaryLocationAndScaleTemporalModel, self.coordinates,
self.dataset,
starting_year=0,
fit_method=self.fit_method)
# Checks that parameters returned are indeed different
mle_params_estimated_year1 = estimator.margin_function_from_fit.get_gev_params(np.array([1])).to_dict()
mle_params_estimated_year3 = estimator.margin_function_from_fit.get_gev_params(np.array([3])).to_dict()
self.assertNotEqual(mle_params_estimated_year1, mle_params_estimated_year3)
if __name__ == '__main__':
unittest.main()
......@@ -3,7 +3,10 @@ import unittest
import numpy as np
import pandas as pd
from experiment.trend_analysis.univariate_test.utils import fitted_linear_margin_estimator
from extreme_fit.estimator.margin_estimator.abstract_margin_estimator import LinearMarginEstimator
from extreme_fit.model.margin_model.linear_margin_model.abstract_temporal_linear_margin_model import \
TemporalMarginFitMethod
from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import StationaryTemporalModel, \
NonStationaryLocationTemporalModel
from extreme_fit.model.utils import r, set_seed_r
......@@ -33,12 +36,13 @@ class TestGevTemporal(unittest.TestCase):
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.fit_method = TemporalMarginFitMethod.is_mev_gev_fit
def test_gev_temporal_margin_fit_stationary(self):
# Create estimator
margin_model = StationaryTemporalModel(self.coordinates)
estimator = LinearMarginEstimator(self.dataset, margin_model)
estimator.fit()
estimator = fitted_linear_margin_estimator(StationaryTemporalModel, self.coordinates, self.dataset,
starting_year=0,
fit_method=self.fit_method)
ref = {'loc': 0.04309190816463247, 'scale': 2.0688696961628437, 'shape': 0.8291528207825063}
for year in range(1, 3):
mle_params_estimated = estimator.margin_function_from_fit.get_gev_params(np.array([year])).to_dict()
......@@ -49,7 +53,6 @@ class TestGevTemporal(unittest.TestCase):
# Create estimator
margin_models = load_non_stationary_temporal_margin_models(self.coordinates)
for margin_model in margin_models:
# margin_model = NonStationaryLocationStationModel(self.coordinates)
estimator = LinearMarginEstimator(self.dataset, margin_model)
estimator.fit()
# Checks that parameters returned are indeed different
......@@ -71,7 +74,8 @@ class TestGevTemporal(unittest.TestCase):
self.assertNotEqual(mle_params_estimated_year5, mle_params_estimated_year3)
def fit_non_stationary_estimator(self, starting_point):
margin_model = NonStationaryLocationTemporalModel(self.coordinates, starting_point=starting_point + self.start_year)
margin_model = NonStationaryLocationTemporalModel(self.coordinates,
starting_point=starting_point + self.start_year)
estimator = LinearMarginEstimator(self.dataset, margin_model)
estimator.fit()
return estimator
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment