[Paper drawing] add drawing for the MCMC sampling

experiment/meteo_france_data/scm_models_data/abstract_study.py

@@ -102,6 +102,10 @@ class AbstractStudy(object):
     def observations_annual_maxima(self) -> AnnualMaxima:
         return AnnualMaxima(df_maxima_gev=pd.DataFrame(self.year_to_annual_maxima, index=self.study_massif_names))
 
+    def annual_maxima_and_years(self, massif_name) -> Tuple[np.ndarray, np.ndarray]:
+        df = self.observations_annual_maxima.df_maxima_gev
+        return df.loc[massif_name].values, np.array(df.columns)
+
     @cached_property
     def year_to_annual_maxima(self) -> OrderedDict:
         # Map each year to an array of size nb_massif

experiment/paper1_steps/1 - non stationary model choice/main2_choice_to_not_use_starting_years.py → experiment/paper1/1 - non stationary model choice/main2_choice_to_not_use_starting_years.py

@@ -2,7 +2,7 @@ import time
 
 from experiment.meteo_france_data.scm_models_data.visualization.hypercube_visualization.altitude_year_hypercube_visualizer import \
     Altitude_Hypercube_Year_Visualizer
-from experiment.paper1_steps.utils import get_full_altitude_visualizer, FULL_ALTITUDES
+from experiment.paper1.utils import get_full_altitude_visualizer, FULL_ALTITUDES
 
 
 def main_fast_spatial_repartition():

experiment/paper1_steps/1 - non stationary model choice/main3_non_stationary_strength_evolution.py → experiment/paper1/1 - non stationary model choice/main3_non_stationary_strength_evolution.py

@@ -9,7 +9,7 @@ from experiment.trend_analysis.univariate_test.gev_trend_test_two_parameters imp
 """
 Visualize the 0.99 quantile initial value and its evolution
 """
-from experiment.paper1_steps.utils import get_full_altitude_visualizer, FULL_ALTITUDES
+from experiment.paper1.utils import get_full_altitude_visualizer, FULL_ALTITUDES
 
 
 def main_fast_spatial_risk_evolution():

experiment/paper1_steps/1 - non stationary model choice/main4_common_spatial_altitude_starting_years_impact.py → experiment/paper1/1 - non stationary model choice/main4_common_spatial_altitude_starting_years_impact.py

@@ -9,7 +9,7 @@ from experiment.trend_analysis.univariate_test.gev_trend_test_two_parameters imp
 """
 Visualize the 0.99 quantile initial value and its evolution
 """
-from experiment.paper1_steps.utils import get_full_altitude_visualizer, FULL_ALTITUDES
+from experiment.paper1.utils import get_full_altitude_visualizer, FULL_ALTITUDES
 
 
 def main_fast_spatial_risk_evolution():

experiment/paper1_steps/1 - non stationary model choice/main4_common_spatial_starting_years_impact.py → experiment/paper1/1 - non stationary model choice/main4_common_spatial_starting_years_impact.py

@@ -9,7 +9,7 @@ from experiment.trend_analysis.univariate_test.gev_trend_test_two_parameters imp
 """
 Visualize the 0.99 quantile initial value and its evolution
 """
-from experiment.paper1_steps.utils import get_full_altitude_visualizer, FULL_ALTITUDES
+from experiment.paper1.utils import get_full_altitude_visualizer, FULL_ALTITUDES
 
 
 def main_fast_spatial_risk_evolution():

experiment/paper1_steps/1 - non stationary model choice/main4_individual_starting_years_impact.py → experiment/paper1/1 - non stationary model choice/main4_individual_starting_years_impact.py

@@ -9,7 +9,7 @@ from experiment.trend_analysis.univariate_test.gev_trend_test_two_parameters imp
 """
 Visualize the 0.99 quantile initial value and its evolution
 """
-from experiment.paper1_steps.utils import get_full_altitude_visualizer, FULL_ALTITUDES
+from experiment.paper1.utils import get_full_altitude_visualizer, FULL_ALTITUDES
 
 
 def main_fast_spatial_risk_evolution():

experiment/paper1/main_drawing_bayesian.py

+import seaborn as sns
+
+import matplotlib.pyplot as plt
+from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSwe3Days
+from experiment.meteo_france_data.scm_models_data.visualization.utils import create_adjusted_axes
+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.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
+from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_extract_eurocode_return_level import \
+    ExtractEurocodeReturnLevelFromMyBayesianExtremes
+from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import \
+    ConfidenceIntervalMethodFromExtremes
+
+
+def main_drawing_bayesian():
+    return_level_bayesian = get_return_level_bayesian_example()
+    print(return_level_bayesian.result_from_fit.df_all_samples)
+    print(return_level_bayesian.result_from_fit.df_posterior_samples)
+    print(return_level_bayesian.posterior_eurocode_return_level_samples_for_temporal_covariate)
+    axes = create_adjusted_axes(1, 3)
+    # Plot the trajectories on the first axis
+    ax_trajectories = axes[0]
+    ax_trajectories_inverted = ax_trajectories.twinx()
+    df_all_samples = return_level_bayesian.result_from_fit.df_all_samples
+    iteration_step = [i + 1 for i in df_all_samples.index]
+    lns = []
+    # Last color is for the return level
+    colors = ['r', 'g', 'b', 'tab:orange']
+    gev_param_name_to_color = dict(zip(GevParams.PARAM_NAMES, colors))
+    gev_param_name_to_ax = dict(
+        zip(GevParams.PARAM_NAMES, [ax_trajectories, ax_trajectories, ax_trajectories_inverted]))
+    gev_param_name_to_label = {n: GevParams.greek_letter_from_gev_param_name(n) for n in GevParams.PARAM_NAMES}
+    for j, gev_param_name in enumerate(GevParams.PARAM_NAMES[:]):
+        label = gev_param_name_to_label[gev_param_name]
+        ax = gev_param_name_to_ax[gev_param_name]
+        color = gev_param_name_to_color[gev_param_name]
+        ln = ax.plot(iteration_step, df_all_samples.iloc[:, j].values, label=label, color=color)
+        lns.extend(ln)
+    ax_trajectories_inverted.set_ylim(-0.3, 10)
+    ax_trajectories.set_ylim(-0.3, 90)
+    for ax in [ax_trajectories, ax_trajectories_inverted]:
+        ax.set_xlim(min(iteration_step), max(iteration_step))
+    # labs = [l.get_label() for l in lns]
+    # ax_trajectories.legend(lns, labs, loc=0)
+    ax.axvline(x=return_level_bayesian.result_from_fit.burn_in_nb, color='k', linestyle='--')
+    ax_trajectories_inverted.legend(loc=1)
+    ax_trajectories.legend(loc=2)
+    ax_trajectories.set_xlabel("MCMC iterations")
+
+    # Plot the parameter posterior on axes 1
+    ax_parameter_posterior = axes[1]
+    ax_parameter_posterior_flip = ax_parameter_posterior.twiny()
+    gev_param_name_to_ax = dict(
+        zip(GevParams.PARAM_NAMES, [ax_parameter_posterior, ax_parameter_posterior.twiny(), ax_parameter_posterior_flip]))
+    df_posterior_samples = return_level_bayesian.result_from_fit.df_posterior_samples
+    lns = []
+    for j, gev_param_name in enumerate(GevParams.PARAM_NAMES[:]):
+        label = gev_param_name_to_label[gev_param_name]
+        color = gev_param_name_to_color[gev_param_name]
+        ax = gev_param_name_to_ax[gev_param_name]
+        ln = sns.kdeplot(df_posterior_samples.iloc[:, j], ax=ax, label=label, color=color)
+        lns.append(ln)
+    labs = [l.get_label() for l in lns]
+    ax_parameter_posterior.legend(lns, labs, loc=0)
+
+    # Plot the return level posterior on the last axes
+    ax_return_level_posterior = axes[2]
+    sns.kdeplot(return_level_bayesian.posterior_eurocode_return_level_samples_for_temporal_covariate,
+                ax=ax_return_level_posterior, color=colors[-1])
+    ax_return_level_posterior.set_xlabel("$z_p(\\theta)")
+    ax_return_level_posterior.set_ylabel("$p(z_p(\\theta)|y)")
+
+    plt.show()
+
+
+def get_return_level_bayesian_example():
+    maxima, years = CrocusSwe3Days(altitude=1800).annual_maxima_and_years('Vercors')
+    # plt.plot(years, maxima)
+    # plt.show()
+    model_class = StationaryTemporalModel
+    coordinates, dataset = load_temporal_coordinates_and_dataset(maxima, years)
+    fitted_estimator = fitted_linear_margin_estimator(model_class, coordinates, dataset, starting_year=1958,
+                                                      fit_method=TemporalMarginFitMethod.extremes_fevd_bayesian)
+    return_level_bayesian = ExtractEurocodeReturnLevelFromMyBayesianExtremes(estimator=fitted_estimator,
+                                                                             ci_method=ConfidenceIntervalMethodFromExtremes.my_bayes,
+                                                                             temporal_covariate=2017)
+    return return_level_bayesian
+
+
+if __name__ == '__main__':
+    main_drawing_bayesian()

experiment/paper1_steps/poster_EVAN2019/main_poster_EVAN2019.py → experiment/paper1/poster_EVAN2019/main_poster_EVAN2019.py

@@ -7,7 +7,7 @@ from experiment.meteo_france_data.scm_models_data.visualization.study_visualizat
 from experiment.trend_analysis.univariate_test.gev_trend_test_one_parameter import GevScaleTrendTest, \
     GevLocationTrendTest
 from experiment.trend_analysis.univariate_test.gev_trend_test_two_parameters import GevLocationAndScaleTrendTest
-from experiment.paper1_steps.utils import get_full_altitude_visualizer
+from experiment.paper1.utils import get_full_altitude_visualizer
 
 POSTER_ALTITUDES = [900, 1800, 2700]
 import matplotlib as mpl

experiment/paper1_steps/poster_EVAN2019/shape_prior_check/some_experiment_EVAN.py → experiment/paper1/poster_EVAN2019/shape_prior_check/some_experiment_EVAN.py

@@ -9,7 +9,7 @@ from experiment.trend_analysis.univariate_test.abstract_comparison_non_stationar
 from experiment.trend_analysis.univariate_test.gev_trend_test_one_parameter import GevScaleTrendTest, \
     GevLocationTrendTest
 from experiment.trend_analysis.univariate_test.gev_trend_test_two_parameters import GevLocationAndScaleTrendTest
-from experiment.paper1_steps.utils import get_full_altitude_visualizer
+from experiment.paper1.utils import get_full_altitude_visualizer
 
 POSTER_ALTITUDES = [900, 1800, 2700]
 import matplotlib as mpl