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

[paper 1] add upper and lower bound into the gev param object. Add script to...

[paper 1] add upper and lower bound into the gev param object. Add script to run fit wihout maximum, and then compare the upper bound with the maximum values that was removed. The test is ok for all altitudes, the upper bound found was always higher than the maximum values
parent 7151bc19
No related merge requests found
Showing with 100 additions and 1 deletion
+100 -1
from typing import Dict
from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal
from experiment.meteo_france_data.scm_models_data.visualization.study_visualization.main_study_visualizer import \
ALL_ALTITUDES_WITHOUT_NAN
from experiment.paper_past_snow_loads.check_mle_convergence_for_trends.without_maximum.study_visualizer_for_fit_witout_maximum import \
StudyVisualizerForFitWithoutMaximum
def fit_without_maximum_value(altitude_to_visualizer: Dict[int, StudyVisualizerForFitWithoutMaximum]):
for v in altitude_to_visualizer.values():
v.maximum_value_test()
if __name__ == '__main__':
altitudes = ALL_ALTITUDES_WITHOUT_NAN[:]
altitude_to_visualizer = {altitude: StudyVisualizerForFitWithoutMaximum(CrocusSnowLoadTotal(altitude=altitude),
multiprocessing=True)
for altitude in altitudes}
fit_without_maximum_value(altitude_to_visualizer)
from typing import Dict, Tuple
import matplotlib
import numpy as np
from cached_property import cached_property
from experiment.meteo_france_data.plot.create_shifted_cmap import get_shifted_map
from experiment.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
from experiment.paper_past_snow_loads.study_visualizer_for_non_stationary_trends import \
StudyVisualizerForNonStationaryTrends
from experiment.trend_analysis.abstract_score import MeanScore
class StudyVisualizerForFitWithoutMaximum(StudyVisualizerForNonStationaryTrends):
def __init__(self, study: AbstractStudy, **kwargs):
super().__init__(study, **kwargs)
@cached_property
def massif_name_to_maximum_index_for_non_null_values(self) -> Tuple[Dict, Dict]:
d = super().massif_name_to_non_null_years_and_maxima
d_without_maximum = {}
d_maximum = {}
for m, (years, maxima) in d.items():
idx = np.argmax(maxima)
maximum = maxima[idx]
maxima = np.delete(maxima, idx)
years = np.delete(years, idx)
d_without_maximum[m] = (years, maxima)
d_maximum[m] = maximum
return d_without_maximum, d_maximum
@property
def massif_name_to_maximum(self) -> Dict:
return self.massif_name_to_maximum_index_for_non_null_values[1]
@cached_property
def massif_name_to_non_null_years_and_maxima(self):
return self.massif_name_to_maximum_index_for_non_null_values[0]
def maximum_value_test(self):
diff = []
for massif_name, maximum in self.massif_name_to_maximum.items():
t = self.massif_name_to_trend_test_that_minimized_aic[massif_name]
msg = '{} {}m'.format(massif_name, self.study.altitude)
upper_bound = t.unconstrained_estimator_gev_params.density_upper_bound
if not np.isinf(upper_bound):
diff.append(upper_bound - maximum)
assert maximum <= upper_bound, msg
if len(diff) > 1:
print('{} mean difference={}'.format(self.study.altitude, min(diff)))
...@@ -130,4 +130,22 @@ class GevParams(AbstractParams): ...@@ -130,4 +130,22 @@ class GevParams(AbstractParams):
if self.shape == 0: if self.shape == 0:
return x return x
else: else:
return np.log(1 + self.shape * x) / self.shape return np.log(1 + self.shape * x) / self.shape
\ No newline at end of file
@property
def bound(self):
return self.location - (self.scale / self.shape)
@property
def density_upper_bound(self):
if self.shape >= 0:
return np.inf
else:
return self.bound
@property
def density_lower_bound(self):
if self.shape <= 0:
return np.inf
else:
return self.bound
\ No newline at end of file
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