Commit ef47242e authored by Le Roux Erwan's avatar Le Roux Erwan
Browse files

[TREND ANALYSIS] add boolean that accept scale 0 parameter that enables to...

[TREND ANALYSIS] add boolean that accept scale 0 parameter that enables to compute easily some mean difference where scale1 is zero
parent d7dd0c85
No related merge requests found
Showing with 28 additions and 14 deletions
+28 -14
...@@ -18,17 +18,19 @@ mpl.rcParams['hatch.linewidth'] = 0.3 ...@@ -18,17 +18,19 @@ mpl.rcParams['hatch.linewidth'] = 0.3
def main_non_stationary_model_comparison(): def main_non_stationary_model_comparison():
stop_loop = False stop_loop = True
for altitude in POSTER_ALTITUDES[:]: for altitude in POSTER_ALTITUDES[:]:
for trend_test_class in [GevLocationTrendTest, GevScaleTrendTest, GevLocationAndScaleTrendTest, ComparisonAgainstMu, ComparisonAgainstSigma][::-1]: for trend_test_class in [GevLocationTrendTest, GevScaleTrendTest, GevLocationAndScaleTrendTest,
ComparisonAgainstMu, ComparisonAgainstSigma][::-1][-1:]:
vizualiser = get_full_altitude_visualizer(Altitude_Hypercube_Year_Visualizer, altitude=altitude, vizualiser = get_full_altitude_visualizer(Altitude_Hypercube_Year_Visualizer, altitude=altitude,
exact_starting_year=1958, reduce_strength_array=False, exact_starting_year=1958, reduce_strength_array=False,
trend_test_class=trend_test_class, trend_test_class=trend_test_class,
) )
# vizualiser.save_to_file = False vizualiser.save_to_file = False
vizualiser.visualize_massif_trend_test_one_altitude(poster_plot=True, write_text_on_massif=False) vizualiser.visualize_massif_trend_test_one_altitude(poster_plot=True, write_text_on_massif=False)
if stop_loop: if stop_loop:
return return
if __name__ == '__main__': if __name__ == '__main__':
main_non_stationary_model_comparison() main_non_stationary_model_comparison()
...@@ -2,6 +2,7 @@ from experiment.trend_analysis.univariate_test.abstract_gev_trend_test import Ab ...@@ -2,6 +2,7 @@ from experiment.trend_analysis.univariate_test.abstract_gev_trend_test import Ab
from experiment.trend_analysis.univariate_test.gev_trend_test_two_parameters import GevLocationAndScaleTrendTest from experiment.trend_analysis.univariate_test.gev_trend_test_two_parameters import GevLocationAndScaleTrendTest
from extreme_estimator.extreme_models.margin_model.temporal_linear_margin_model import \ from extreme_estimator.extreme_models.margin_model.temporal_linear_margin_model import \
NonStationaryLocationAndScaleModel, NonStationaryLocationStationModel, NonStationaryScaleStationModel NonStationaryLocationAndScaleModel, NonStationaryLocationStationModel, NonStationaryScaleStationModel
import numpy as np
class AbstractComparisonNonStationaryModel(AbstractGevTrendTest): class AbstractComparisonNonStationaryModel(AbstractGevTrendTest):
...@@ -17,6 +18,12 @@ class AbstractComparisonNonStationaryModelOneParameter(AbstractComparisonNonStat ...@@ -17,6 +18,12 @@ class AbstractComparisonNonStationaryModelOneParameter(AbstractComparisonNonStat
def degree_freedom_chi2(self) -> int: def degree_freedom_chi2(self) -> int:
return 1 return 1
@property
def test_sign(self) -> int:
# Test sign correspond to the difference between the 2 likelihoods
# Therefore, colors sum up which non stationary model explain best the data
return np.sign(self.likelihood_ratio)
class ComparisonAgainstMu(AbstractComparisonNonStationaryModelOneParameter, GevLocationAndScaleTrendTest): class ComparisonAgainstMu(AbstractComparisonNonStationaryModelOneParameter, GevLocationAndScaleTrendTest):
......
import numpy as np import numpy as np
import pandas as pd import pandas as pd
from cached_property import cached_property
from scipy.stats import chi2 from scipy.stats import chi2
from experiment.trend_analysis.univariate_test.abstract_univariate_test import AbstractUnivariateTest from experiment.trend_analysis.univariate_test.abstract_univariate_test import AbstractUnivariateTest
...@@ -109,11 +110,11 @@ class AbstractGevTrendTest(AbstractUnivariateTest): ...@@ -109,11 +110,11 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
def test_sign(self) -> int: def test_sign(self) -> int:
return np.sign(self.test_trend_slope_strength) return np.sign(self.test_trend_slope_strength)
def get_non_stationary_linear_coef(self, gev_param_name): def get_non_stationary_linear_coef(self, gev_param_name: str):
return self.non_stationary_estimator.margin_function_fitted.get_coef(gev_param_name, return self.non_stationary_estimator.margin_function_fitted.get_coef(gev_param_name,
AbstractCoordinates.COORDINATE_T) AbstractCoordinates.COORDINATE_T)
@property @cached_property
def non_stationary_constant_gev_params(self) -> GevParams: def non_stationary_constant_gev_params(self) -> GevParams:
# Constant parameters correspond to the gev params in 1958 # Constant parameters correspond to the gev params in 1958
return self.non_stationary_estimator.margin_function_fitted.get_gev_params(coordinate=np.array([1958]), return self.non_stationary_estimator.margin_function_fitted.get_gev_params(coordinate=np.array([1958]),
...@@ -145,8 +146,8 @@ class AbstractGevTrendTest(AbstractUnivariateTest): ...@@ -145,8 +146,8 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
def variance_difference_same_sign_as_slope_strenght(self) -> bool: def variance_difference_same_sign_as_slope_strenght(self) -> bool:
return False return False
def mean_difference(self, zeta0, mu1=0.0, sigma1=0.0): def mean_difference(self, zeta0: float, mu1: float = 0.0, sigma1: float = 0.0) -> float:
return GevParams(loc=mu1, scale=sigma1, shape=zeta0).mean return GevParams(loc=mu1, scale=sigma1, shape=zeta0, accept_zero_scale_parameter=True).mean
@property @property
def test_trend_constant_quantile(self): def test_trend_constant_quantile(self):
......
...@@ -14,7 +14,7 @@ class ExtremeParams(AbstractParams, ABC): ...@@ -14,7 +14,7 @@ class ExtremeParams(AbstractParams, ABC):
self.location = loc self.location = loc
self.scale = scale self.scale = scale
self.shape = shape self.shape = shape
# Scale cannot be negative # By default, scale cannot be negative
# (sometimes it happens, when we want to find a quantile for every point of a 2D map # (sometimes it happens, when we want to find a quantile for every point of a 2D map
# then it can happen that a corner point that was not used for fitting correspond to a negative scale, # then it can happen that a corner point that was not used for fitting correspond to a negative scale,
# in the case we set all the parameters as equal to np.nan, and we will not display those points) # in the case we set all the parameters as equal to np.nan, and we will not display those points)
......
from collections import OrderedDict from collections import OrderedDict
from typing import List
from cached_property import cached_property from cached_property import cached_property
...@@ -15,9 +16,11 @@ class GevParams(ExtremeParams): ...@@ -15,9 +16,11 @@ class GevParams(ExtremeParams):
SUMMARY_NAMES = PARAM_NAMES + ExtremeParams.QUANTILE_NAMES SUMMARY_NAMES = PARAM_NAMES + ExtremeParams.QUANTILE_NAMES
NB_SUMMARY_NAMES = len(SUMMARY_NAMES) NB_SUMMARY_NAMES = len(SUMMARY_NAMES)
def __init__(self, loc: float, scale: float, shape: float, block_size: int = None): def __init__(self, loc: float, scale: float, shape: float, block_size: int = None, accept_zero_scale_parameter=False):
super().__init__(loc, scale, shape) super().__init__(loc, scale, shape)
self.block_size = block_size self.block_size = block_size
if accept_zero_scale_parameter and scale == 0.0:
self.has_undefined_parameters = False
def quantile(self, p) -> float: def quantile(self, p) -> float:
if self.has_undefined_parameters: if self.has_undefined_parameters:
...@@ -46,6 +49,7 @@ class GevParams(ExtremeParams): ...@@ -46,6 +49,7 @@ class GevParams(ExtremeParams):
return self.to_dict().__str__() return self.to_dict().__str__()
def quantile_strength_evolution(self, p=0.99, mu1=0.0, sigma1=0.0): def quantile_strength_evolution(self, p=0.99, mu1=0.0, sigma1=0.0):
# todo: chagne the name to time derivative
""" """
Compute the variation of some quantile with respect to time. Compute the variation of some quantile with respect to time.
(when mu1 and sigma1 can be modified with time) (when mu1 and sigma1 can be modified with time)
...@@ -63,13 +67,13 @@ class GevParams(ExtremeParams): ...@@ -63,13 +67,13 @@ class GevParams(ExtremeParams):
# Compute some indicators (such as the mean and the variance) # Compute some indicators (such as the mean and the variance)
def g(self, k): def g(self, k) -> float:
# Compute the g_k parameters as defined in wikipedia # Compute the g_k parameters as defined in wikipedia
# https://fr.wikipedia.org/wiki/Loi_d%27extremum_g%C3%A9n%C3%A9ralis%C3%A9e # https://fr.wikipedia.org/wiki/Loi_d%27extremum_g%C3%A9n%C3%A9ralis%C3%A9e
return gamma(1 - k * self.shape) return gamma(1 - k * self.shape)
@property @property
def mean(self): def mean(self) -> float:
if self.has_undefined_parameters: if self.has_undefined_parameters:
return np.nan return np.nan
elif self.shape >= 1: elif self.shape >= 1:
...@@ -78,7 +82,7 @@ class GevParams(ExtremeParams): ...@@ -78,7 +82,7 @@ class GevParams(ExtremeParams):
return self.location + self.scale * (self.g(k=1) - 1) / self.shape return self.location + self.scale * (self.g(k=1) - 1) / self.shape
@property @property
def variance(self): def variance(self) -> float:
if self.has_undefined_parameters: if self.has_undefined_parameters:
return np.nan return np.nan
elif self.shape >= 0.5: elif self.shape >= 0.5:
...@@ -87,11 +91,11 @@ class GevParams(ExtremeParams): ...@@ -87,11 +91,11 @@ class GevParams(ExtremeParams):
return ((self.scale / self.shape) ** 2) * (self.g(k=2) - self.g(k=1) ** 2) return ((self.scale / self.shape) ** 2) * (self.g(k=2) - self.g(k=1) ** 2)
@property @property
def std(self): def std(self) -> float:
return np.sqrt(self.variance) return np.sqrt(self.variance)
@classmethod @classmethod
def indicator_names(cls): def indicator_names(cls) -> List[str]:
return ['mean', 'std'] + cls.QUANTILE_NAMES[:2] return ['mean', 'std'] + cls.QUANTILE_NAMES[:2]
@cached_property @cached_property
......
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