import numpy as np
#np.bool = np.bool_
#import seaborn as sns
#import pandas as pd
import matplotlib.pyplot as plt
import pymc as pm
import arviz as az
import pytensor.tensor as at
import pymc_experimental.distributions as pmx
from scipy.stats import genextreme as gev

WARNING (pytensor.tensor.blas): Using NumPy C-API based implementation for BLAS functions.

from CleanGeslaBM import BlockMaxima
from ModelSetup import model_setup
from HilbertSpaceLatLon import HSGP_LatLon
from CovarianceSphericalMetric import Matern52Chordal # ExponentialChordal, Matern32Chordal
from BaselineModel import Baseline
from DataAnalysisUtil import train_test_split_with_info, LOO_interpolation_validation, lat_lon_in_BM

from HilbertSpaceLatLon import calculate_K, calculate_Kapprox, approx_hsgp_hyperparams_xyz

from DataAnalysisPlt import  BHM_GeoScatter_train_and_new_sites_plot, BHM_mean_GeoScatter_train_test_split,Baseline_GEV_param_KDE_comparison_plot, Baseline_GEV_param_resampeling_KDE_plots

from DataAnalysisPlt import (BHM_boxplot_models_comparison,BHM_subset_pair_plot, BHM_supp_GEV_KDE_plot, BHM_geographical_spread_GEV_param_posterior_plot) # BHM_violinplot_models_comparison, BHM_subset_denisty_comparsion_plot,

from DataAnalysisPlt import BHM_hdpi_and_MLM_ci_return_lvl_multiple_sites_plot, BHM_violinplot_models_comparison_old # BHM_hpdi_and_MLM_ci_return_lvl_Loo_multiple_models_plot, BHM_GEV_param_LOO_posterior_plot,LOO_pair_plot,

from DataAnalysisPlt import BHM_comparative_mean_and_hpdi_var_name_surface_plots #,BHM_interpolated_train_test_var_name_surface_mean_or_std, BHM_comparative_mean_and_hpdi_var_name_surface_shared_cbar_plots

RANDOM_SEED = 2024
rng = np.random.default_rng(RANDOM_SEED)
#GeslaDataset object, where datapath contains all GESLA3 data with respect to some region 
metafile = "/home/sm_seuba/Documents/GESLA3_ALL_2.csv"
region_datapath = "/home/sm_seuba/Documents/GESLA3 DATA BALTIC NORTH CMEMS/"
#region_datapath = "/home/sm_seuba/Documents/GESLA3 DATA BALTIC CMEMS/"
#region_datapath = "/home/sm_seuba/Documents/GESLA3 DATA SWE CMEMS/"
#region_datapath = "/home/sm_seuba/Documents/GESLA3 DATA FIN SWE CMEMS/"

region_BM, region_sites =  BlockMaxima(metafile, region_datapath, start=1850, detrend=20, window=True, center=True)

region_base_BM = region_BM.dropna(axis=1, thresh= 20)

#We remove one of the four stations in the top4 subset
#Loo_column, XY_Loo = LOO_interpolation_validation(region_base_BM, max_noNaNs_idx=1) #0) #3) #2)
#region_base_BM.drop(columns=[Loo_column], inplace=True)

#train, train_info, test,test_info = train_test_split_with_info(region_base_BM,rows=55, thresh=35, case=0)
#__, region_base_info, _, _ = train_test_split_with_info(region_BM, rows=0, thresh=10, case=1)

xy_base, xs_base, ys_base = lat_lon_in_BM(region_base_BM)

max(ys_base), max(ys_base)

(65.696899, 65.696899)

np.arange(53.80, 66.70,  0.99, dtype="float32").shape, np.arange(9.60,30.00,  1.5, dtype="float32").shape

((14,), (14,))

x_new = np.arange(53.80, 66.70, 0.99, dtype="float32")
y_new = np.arange(9.60, 30.00,  1.5, dtype="float32")
X, Y = np.meshgrid(x_new, y_new, indexing="ij")

X.shape,Y.shape

((14, 14), (14, 14))

Xs, Ys = X.flatten(), Y.flatten()
XY_test_longer = [(Xs[i], Ys[i]) for i in range(Xs.shape[0])]

#Loo_column, XY_Loo = LOO_interpolation_validation(region_base_BM, max_noNaNs_idx=0) #1) #2) #3)
#region_base_BM.drop(columns=[Loo_column], inplace=True)

train, train_info, test, test_info = train_test_split_with_info(region_base_BM,rows=55, thresh=35, case=0) #45 30
base_case_BM, region_base_info, _, _ = train_test_split_with_info(region_base_BM, rows=0, thresh=0, case=1) #train_test_split_with_info(region_BM, rows=0, thresh=10, case=1)
#__, region_base_info, _, _ = train_test_split_with_info(region_BM, rows=0, thresh=10, case=1)

train_info

{'Case ': 0,
 'Row': 55,
 'Tresh': 35,
 'BM shape': (79, 36),
 '#Entrys': 2844,
 '#NaNs': 864,
 'NaNs/Entrys': 0.3038}

test_info

{'Case ': 0,
 'Row': 0,
 'Tresh': 0,
 'BM shape': (134, 15),
 '#Entrys': 2010,
 '#NaNs': 1545,
 'NaNs/Entrys': 0.7687}

region_base_info

{'Case ': 1,
 'Row': 0,
 'Tresh': 0,
 'BM shape': (134, 51),
 '#Entrys': 6834,
 '#NaNs': 3947,
 'NaNs/Entrys': 0.5776}

BHM_GeoScatter_train_and_new_sites_plot(train, XY_test_longer,static=True)

BHM_mean_GeoScatter_train_test_split(region_base_BM, train, test, train_test_split="train", static=True)

BHM_mean_GeoScatter_train_test_split(region_base_BM, train, test, train_test_split="both", static=True)

XY_train, coord_train, dim_train, data_train = model_setup(train)
#XY_test, coord_test, dim_test, data_test = model_setup(region_BM[LOO_interpolation(region_base_BM, max_noNaNs=)])
#train_sites, test_sites = [i for i in range(dim_train[1])],  [dim_train[1] +1 ]# for i in range(dim_test[1])]
XY_test_longer_array = np.array(XY_test_longer)
train_sites, LOO_sites, test_longer_sites = [i for i in range(dim_train[1])], [dim_train[1] +1 ], [i + dim_train[1] for i in range(len(XY_test_longer)) ]

Baseline_data_train = Baseline(data_train)
Baseline_region_BM = Baseline(region_base_BM)

Baseline_data_train_stations = list(Baseline_data_train.index)

Baseline_GEV_param_KDE_comparison_plot(Baseline_data_train, Baseline_region_BM,MLE_consistent=False)

Baseline_GEV_param_KDE_comparison_plot(Baseline_data_train, Baseline_region_BM,MLE_consistent=True,CI_consistent=True)

Baseline_GEV_param_resampeling_KDE_plots(data_train, 10000,Baseline_region_BM)

not_consistent_CI_MLE = list(Baseline_data_train[(Baseline_data_train["xi_i"] <= -0.50)].index)
consistent_CI_MLE = list(Baseline_data_train[(Baseline_data_train["xi_i"] >-0.50) &(Baseline_data_train["notNaNs"] > 60)].index)
Baseline_consistent_CI_MLE = Baseline_data_train[(Baseline_data_train["xi_i"] > -0.50)]
train_columns = list(train.columns)

not_consistent_CI_MLE_region = list(Baseline_region_BM[(Baseline_region_BM["xi_i"] <= -0.50)].index)
consistent_CI_MLE_region = list(Baseline_region_BM[(Baseline_region_BM["xi_i"] >-0.50) & (Baseline_region_BM["notNaNs"] > 100)].index)
Baseline_consistent_CI_MLE_region = Baseline_region_BM[(Baseline_region_BM["xi_i"] > -0.50)]
train_columns_region = list(region_base_BM)

fig, axes = plt.subplots(nrows=2, ncols=4, figsize=(26,10))
row_count = 0
column_count = 0
for coord in consistent_CI_MLE:
    if column_count == 4:
        row_count +=1
        column_count = 0
    mu_param = Baseline_data_train["mu_i"].loc[[coord]].to_numpy(dtype="float")[0]
    sigma_param = Baseline_data_train["sigma_i"].loc[[coord]].to_numpy(dtype="float")[0]
    xi_param = np.round(Baseline_data_train["xi_i"].loc[[coord]].to_numpy(dtype="float"),5)[0]
    supp_GEV = Baseline_data_train["supp_GEV"].loc[[coord]].to_numpy(dtype="float")[0]
    years, values = list(train[coord].dropna().index), train[coord].dropna().to_numpy(dtype="float")
    values = train[coord].dropna().to_numpy(dtype="float")
    n_obs = Baseline_data_train["notNaNs"].loc[[coord]].to_numpy(dtype="float")[0]
    train[coord].dropna().plot(ax=axes[row_count, column_count]).scatter(y=values, x=years, color="C1")
    axes[row_count, column_count].axhline(mu_param, color="black",  label="MLE mu")
    axes[row_count, column_count].axhline(supp_GEV, color="black",  linestyle="dashed",label="MLE supp GEV")
    axes[row_count, column_count].axhline(mu_param+sigma_param, linestyle="dotted", color="red",  label="MLE mu +- sigma")
    axes[row_count, column_count].axhline(mu_param-sigma_param, linestyle="dotted",color="red")
    axes[row_count, column_count].set_title(f"{region_sites[coord]} #obs: {n_obs} GEV xi: {xi_param}", fontsize=15)
    axes[row_count, column_count].legend(loc='upper left')
    column_count +=1 
for axes in axes.flatten():
        if not axes.has_data():
            fig.delaxes(axes)
fig.suptitle("Stations with consistent MLE and CI", fontsize=30)
fig.show()

mu_param = Baseline_data_train["mu_i"].loc[[coord]].to_numpy(dtype="float")[0]
sigma_param = Baseline_data_train["sigma_i"].loc[[coord]].to_numpy(dtype="float")[0]
xi_param = np.round(Baseline_data_train["xi_i"].loc[[coord]].to_numpy(dtype="float"),5)[0]

fig, axes = plt.subplots(nrows=2, ncols=4, figsize=(26,10))
row_count = 0
column_count = 0
for coord in consistent_CI_MLE:
    if column_count == 4:
        row_count +=1
        column_count = 0
    xi_param = np.round(Baseline_data_train["xi_i"].loc[[coord]].to_numpy(dtype="float"),5)[0]
    mu_param = Baseline_data_train["mu_i"].loc[[coord]].to_numpy(dtype="float")[0]
    sigma_param = Baseline_data_train["sigma_i"].loc[[coord]].to_numpy(dtype="float")[0]
    supp_GEV = Baseline_data_train["supp_GEV"].loc[[coord]].to_numpy(dtype="float")[0]
    x_samples = np.arange(0, stop=supp_GEV, step=0.1)
    n_obs = Baseline_data_train["notNaNs"].loc[[coord]].to_numpy(dtype="float")[0]
    train[coord].dropna().plot(kind="kde",ax=axes[row_count, column_count], label=f"MLE xi {xi_param} ")
    axes[row_count, column_count].scatter(x=mu_param, y=0, c="C1")
    axes[row_count, column_count].axvline(supp_GEV, color="black",  linestyle="dashed",label="MLE supp GEV")
    axes[row_count, column_count].plot(x_samples,gev.pdf(x_samples, -1*xi_param, loc=mu_param, scale=sigma_param), color="C1", linestyle="dashed",label="pdf by MLE" )
    axes[row_count, column_count].legend(loc='upper left', fontsize=13)
    axes[row_count, column_count].set_title(f"{region_sites[coord]} station  #obs {n_obs}", fontsize=15) #({train_columns.index(coord)}), obs {n_obs}", fontsize=15)
    axes[row_count, column_count].tick_params(axis="both", labelsize=13)
    column_count +=1
for axes in axes.flatten():
        if not axes.has_data():
            fig.delaxes(axes)
fig.suptitle("Stations with consistent MLE and CI (training data)", fontsize=30)
fig.show()

fig, axes = plt.subplots(nrows=2, ncols=4,figsize=(26,10))
row_count = 0
column_count = 0
for coord in consistent_CI_MLE_region:
    if column_count == 4:
        row_count +=1
        column_count = 0
    mu_param = Baseline_region_BM["mu_i"].loc[[coord]].to_numpy(dtype="float")[0]
    sigma_param = Baseline_region_BM["sigma_i"].loc[[coord]].to_numpy(dtype="float")[0]
    xi_param = np.round(Baseline_region_BM["xi_i"].loc[[coord]].to_numpy(dtype="float"),5)[0]
    supp_GEV = Baseline_region_BM["supp_GEV"].loc[[coord]].to_numpy(dtype="float")[0]
    years, values = list(region_base_BM.dropna(axis=1, thresh=10)[coord].dropna().index), region_base_BM.dropna(axis=1, thresh=10)[coord].dropna().to_numpy(dtype="float")
    values = region_base_BM.dropna(axis=1, thresh=10)[coord].dropna().to_numpy(dtype="float")
    n_obs = Baseline_region_BM["notNaNs"].loc[[coord]].to_numpy(dtype="float")[0]
    region_base_BM.dropna(axis=1, thresh=10)[coord].dropna().plot(ax=axes[row_count, column_count]).scatter(y=values, x=years, color="C1")
    axes[row_count, column_count].axhline(mu_param, color="black",  label="MLE mu")
    axes[row_count, column_count].axhline(supp_GEV, color="black",  linestyle="dashed",label="MLE supp GEV")
    axes[row_count, column_count].axhline(mu_param+sigma_param, linestyle="dotted", color="red",  label="MLE mu +- sigma")
    axes[row_count, column_count].axhline(mu_param-sigma_param, linestyle="dotted",color="red")
    axes[row_count, column_count].set_title(f"{region_sites[coord]} #obs: {n_obs} GEV xi: {xi_param}", fontsize=15)
    axes[row_count, column_count].legend(loc='upper left')
    column_count +=1 
for axes in axes.flatten():
        if not axes.has_data():
            fig.delaxes(axes)
fig.suptitle("Stations with consistent MLE and or CI region", fontsize=30)
fig.show()

fig, axes = plt.subplots(nrows=2, ncols=4, figsize=(26,10))
row_count = 0
column_count = 0
for coord in consistent_CI_MLE_region:
    if column_count == 4:
        row_count +=1
        column_count = 0
    xi_param = np.round(Baseline_region_BM["xi_i"].loc[[coord]].to_numpy(dtype="float"),5)[0]
    sigma_param = Baseline_region_BM["sigma_i"].loc[[coord]].to_numpy(dtype="float")[0]
    mu_param = Baseline_region_BM["mu_i"].loc[[coord]].to_numpy(dtype="float")[0]
    supp_GEV = Baseline_region_BM["supp_GEV"].loc[[coord]].to_numpy(dtype="float")[0]
    n_obs = Baseline_region_BM["notNaNs"].loc[[coord]].to_numpy(dtype="float")[0]
    x_samples = np.arange(0, stop=supp_GEV, step=0.1)
    region_base_BM[coord].dropna().plot(kind="kde",ax=axes[row_count, column_count], label=f"MLE xi {xi_param} ")
    axes[row_count, column_count].scatter(x=mu_param, y=0, c="C1")
    axes[row_count, column_count].axvline(supp_GEV, color="black",  linestyle="dashed",label="MLE supp GEV")
    axes[row_count, column_count].plot(x_samples,gev.pdf(x_samples, -1*xi_param, loc=mu_param, scale=sigma_param), color="C1", linestyle="dashed",label="pdf by MLE" )
    axes[row_count, column_count].legend(loc='upper left', fontsize=13)
    axes[row_count, column_count].set_title(f"{region_sites[coord]} station #obs {n_obs}", fontsize=15) #({train_columns_region.index(coord)}), #obs {n_obs}", fontsize=15)
    axes[row_count, column_count].tick_params(axis="both", labelsize=13)
    column_count +=1
for axes in axes.flatten():
        if not axes.has_data():
            fig.delaxes(axes)
fig.suptitle("Stations with consistent MLE and CI (region)", fontsize=30)
fig.show()

XY_base,_ ,__,____ = model_setup(region_base_BM)

approx_hsgp_hyperparams_xyz(XY_base,  lengthscale=350, variable_name="xi")

For the variable xi
Recommended smallest number of basis vectors [Xs,Ys,Zs] (m): [11, 11, 11]
Recommended smallest scaling factor [Xs,Ys,Zs] (c): [2.4, 6.6, 1.6]

approx_hsgp_hyperparams_xyz(XY_base,  lengthscale=350, variable_name="mu")

For the variable mu
Recommended smallest number of basis vectors [Xs,Ys,Zs] (m): [11, 11, 11]
Recommended smallest scaling factor [Xs,Ys,Zs] (c): [2.4, 6.6, 1.6]

approx_hsgp_hyperparams_xyz(XY_base,  lengthscale=350, variable_name="sigma")

For the variable sigma
Recommended smallest number of basis vectors [Xs,Ys,Zs] (m): [11, 11, 11]
Recommended smallest scaling factor [Xs,Ys,Zs] (c): [2.4, 6.6, 1.6]

chosen_ell = 350
fig, axs = plt.subplots(2, 4, figsize=(14, 7), sharey=True)
K = calculate_K(XY_base, chosen_ell)
axs[0, 0].imshow(K, cmap="inferno", vmin=0, vmax=1)
axs[0, 0].set(xlabel="x1", ylabel="x2", title=f"True Gram matrix\nTrue $\\ell$ = {chosen_ell}")
axs[1, 0].axis("off")
im_kwargs = {
    "cmap": "inferno",
    "vmin": 0,
    "vmax": 1,
    "interpolation": "none",
}

## column 1

m, c = [11, 11, 11],  2.2
m_star = np.prod(m)
K_approx = calculate_Kapprox(XY_base, c, m,chosen_ell)
axs[0, 1].imshow(K_approx, **im_kwargs)
axs[0, 1].set_title(f"m = {m}, c = {c}, m* = {m_star}")

mm, c = [12, 11, 20], 2.8
m_star = np.prod(m)
K_approx = calculate_Kapprox(XY_base, c, m, chosen_ell)
axs[1, 1].imshow(K_approx, **im_kwargs)
axs[1, 1].set(xlabel="x1", ylabel="x2", title=f"m = {m}, c = {c},  m* = {m_star}")

## column 2

#m, c = [40, 11, 40], 6.6
m, c = [12, 11, 26], 3.0
m_star = np.prod(m)
K_approx = calculate_Kapprox(XY_base, c, m, chosen_ell)
axs[0, 2].imshow(K_approx, **im_kwargs)
axs[0, 2].set_title(f"m = {m}, c = {c}, m* = {m_star}")



m, c = [14, 11, 26], 3.0
m_star = np.prod(m)
K_approx = calculate_Kapprox(XY_base, c, m,chosen_ell)
axs[1, 2].imshow(K_approx, **im_kwargs)
axs[1, 2].set_title(f"m = {m},  c = {c}, m* = {m_star} ")

## column 3
m, c = [35, 11, 35], 6.6
m_star = np.prod(m)
K_approx = calculate_Kapprox(XY_base, c, m,chosen_ell)
axs[0, 3].imshow(K_approx, **im_kwargs)
axs[0, 3].set_title(f"m = {m},  c = {c}, m* = {m_star}")

m, c = [12, 11, 20], 3.1
m_star = np.prod(m)
K_approx = calculate_Kapprox(XY_base, c, m,chosen_ell)
axs[1, 3].imshow(K_approx, **im_kwargs)
axs[1, 3].set_title(f"m = {m},  c = {c}, m* = {m_star}")

for ax in axs.flatten():
    ax.grid(False)

fig.suptitle(f'The true Gram maxtrix for cov_xi: Matern52Chordal with ls={chosen_ell} and the HS approximated Gram, where n={train_info["#Entrys"]} and n^2={train_info["#Entrys"]**2} ', fontsize=16)
plt.tight_layout();

chosen_ell = 350
fig, axs = plt.subplots(2, 4, figsize=(14, 7), sharey=True)
K = calculate_K(XY_base, chosen_ell)
axs[0, 0].imshow(K, cmap="inferno", vmin=0, vmax=1)
axs[0, 0].set(xlabel="x1", ylabel="x2", title=f"True Gram matrix\nTrue $\\ell$ = {chosen_ell}")
axs[1, 0].axis("off")
im_kwargs = {
    "cmap": "inferno",
    "vmin": 0,
    "vmax": 1,
    "interpolation": "none",
}

## column 1

m, c = [11, 11, 11],  2.2
m_star = np.prod(m)
K_approx = calculate_Kapprox(XY_base, c, m,chosen_ell)
axs[0, 1].imshow(K_approx, **im_kwargs)
axs[0, 1].set_title(f"m = {m}, c = {c}, m* = {m_star}")

mm, c = [12, 11, 20], 2.8
m_star = np.prod(m)
K_approx = calculate_Kapprox(XY_base, c, m, chosen_ell)
axs[1, 1].imshow(K_approx, **im_kwargs)
axs[1, 1].set(xlabel="x1", ylabel="x2", title=f"m = {m}, c = {c},  m* = {m_star}")

## column 2

#m, c = [40, 11, 40], 6.6
m, c = [12, 11, 26], 3.0
m_star = np.prod(m)
K_approx = calculate_Kapprox(XY_base, c, m, chosen_ell)
axs[0, 2].imshow(K_approx, **im_kwargs)
axs[0, 2].set_title(f"m = {m}, c = {c}, m* = {m_star}")



m, c = [14, 11, 26], 3.0
m_star = np.prod(m)
K_approx = calculate_Kapprox(XY_base, c, m,chosen_ell)
axs[1, 2].imshow(K_approx, **im_kwargs)
axs[1, 2].set_title(f"m = {m},  c = {c}, m* = {m_star} ")

## column 3
m, c = [35, 11, 35], 6.6
m_star = np.prod(m)
K_approx = calculate_Kapprox(XY_base, c, m,chosen_ell)
axs[0, 3].imshow(K_approx, **im_kwargs)
axs[0, 3].set_title(f"m = {m},  c = {c}, m* = {m_star}")

m, c = [12, 11, 20], 3.1
m_star = np.prod(m)
K_approx = calculate_Kapprox(XY_base, c, m,chosen_ell)
axs[1, 3].imshow(K_approx, **im_kwargs)
axs[1, 3].set_title(f"m = {m},  c = {c}, m* = {m_star}")

for ax in axs.flatten():
    ax.grid(False)

fig.suptitle(f'The true Gram maxtrix for cov_mu: Matern52Chordal with ls={chosen_ell} and the HS approximated Gram, where n={train_info["#Entrys"]} and n^2={train_info["#Entrys"]**2} ', fontsize=16)
plt.tight_layout();

chosen_ell= 350
fig, axs = plt.subplots(2, 4, figsize=(14, 7), sharey=True)
K = calculate_K(XY_base,chosen_ell)
axs[0, 0].imshow(K, cmap="inferno", vmin=0, vmax=1)
axs[0, 0].set(xlabel="x1", ylabel="x2", title=f"True Gram matrix\nTrue $\\ell$ = {chosen_ell}")
axs[1, 0].axis("off")
im_kwargs = {
    "cmap": "inferno",
    "vmin": 0,
    "vmax": 1,
    "interpolation": "none",
}

## column 1


m, c = [11, 11, 11],  2.2
m_star = np.prod(m)
K_approx = calculate_Kapprox(XY_base, c, m,chosen_ell)
axs[0, 1].imshow(K_approx, **im_kwargs)
axs[0, 1].set_title(f"m = {m}, c = {c}, m* = {m_star}")


m, c = [12, 11, 20], 2.8
m_star = np.prod(m)
K_approx = calculate_Kapprox(XY_base, c, m, chosen_ell)
axs[1, 1].imshow(K_approx, **im_kwargs)
axs[1, 1].set(xlabel="x1", ylabel="x2", title=f"m = {m}, c = {c},  m* = {m_star}")

## column 2


m, c = [12, 11, 22], 3.2
m_star = np.prod(m)
K_approx = calculate_Kapprox(XY_base, c, m, chosen_ell)
axs[0, 2].imshow(K_approx, **im_kwargs)
axs[0, 2].set_title(f"m = {m}, c = {c}, m* = {m_star}")


m, c = [14, 11, 22], 3.2
m_star = np.prod(m)
K_approx = calculate_Kapprox(XY_base, c, m,chosen_ell)
axs[1, 2].imshow(K_approx, **im_kwargs)
axs[1, 2].set_title(f"m = {m},  c = {c}, m* = {m_star} ")


## column 3
m, c = [35, 11,50], 6.2
m_star = np.prod(m)
K_approx = calculate_Kapprox(XY_base, c, m,chosen_ell)
axs[0, 3].imshow(K_approx, **im_kwargs)
axs[0, 3].set_title(f"m = {m},  c = {c}, m* = {m_star}")

m, c = [12, 11, 20], 3.1
m_star = np.prod(m)
K_approx = calculate_Kapprox(XY_base, c, m,chosen_ell)
axs[1, 3].imshow(K_approx, **im_kwargs)
axs[1, 3].set_title(f"m = {m},  c = {c}, m* = {m_star}")

for ax in axs.flatten():
    ax.grid(False)

fig.suptitle(f'The true Gram maxtrix for cov_sigma: Matern52Chordal with ls={chosen_ell} and HS approximated ', fontsize=16)
plt.tight_layout();

compare_column_0,_ =  LOO_interpolation_validation(region_base_BM, max_noNaNs_idx=0) #(56.105301, 15.5894) #max_noNaNs_idx=0
compare_column_1,_ =  LOO_interpolation_validation(region_base_BM, max_noNaNs_idx=1) #(59.3242, 18.081699) #max_noNaNs_idx=1
compare_column_2,_ =  LOO_interpolation_validation(region_base_BM, max_noNaNs_idx=2)
compare_column_3,_ =  LOO_interpolation_validation(region_base_BM, max_noNaNs_idx=3)
compare_column_top4_list = [compare_column_0, compare_column_1, compare_column_2, compare_column_3]
columns_list = list(train.columns)
compare_column_idx_0, compare_column_idx_1, compare_column_idx_2, compare_column_idx_3 = columns_list.index(compare_column_0),columns_list.index(compare_column_1),columns_list.index(compare_column_2), columns_list.index(compare_column_3)
compare_stations_in_top_4_dict = {"station": [compare_column_idx_0,compare_column_idx_1, compare_column_idx_2, compare_column_idx_3], "location": "Top 4 (fewest NaNs) stations"}

geographical_spread = {"station": [0,1,2,6,9,18,13,24,22,29,35]}
geographical_spread_Skagerrak = {"station": [13, 18, 22, 24, 29, 35], "location": "Skagerrak"}
geographical_spread_CoastOfSkane = {"station": [0, 9, 16, 21, 28], "location": "The coastline of Scania County"}
geographical_spread_CostOfSthlmToSmaland = {"station": [3, 10, 11, 15, 31, 33], "location": "The costline from Stockholm to Smaland County" }
geographical_spread_GulfOfFinland = {"station": [1, 2, 6], "location": "Gulf of Finland"}
geographical_spread_BothniaBaySouth = {"station": [5, 7, 12, 14, 19, 26, 27, 32],"location": "Bothnia Bay South"}
geographical_spread_BothniaBayNorth = {"station": [4, 8, 17, 20, 23,25, 30, 34], "location": "Bothnia Bay North"}
geographical_spread_TheTopOneOfEachZone = {"station": [0, 1, 15, 17, 26, 35],  "location": "Top one (fewest NaNs) station from each zone"}

geographical_spread_list = [geographical_spread_Skagerrak,geographical_spread_CoastOfSkane, geographical_spread_CostOfSthlmToSmaland, 
                            geographical_spread_GulfOfFinland, geographical_spread_BothniaBaySouth, geographical_spread_BothniaBayNorth] #geographical_spread_TheTopOneOfEachZone]

p5, p10, p20, p30, p50, p100, p200, p500, p1000,p5000, p10000, p50000, p100000 = 1/5, 1/10, 1/20, 1/30, 1/50, 1/100,1/200, 1/500, 1/1000,1/5000, 1/10000, 1/50000, 1/100000




ls_mu = 350
ls_sigma = ls_mu #325
ls_xi = ls_mu# 300

m, c = [12, 11, 20], 2.8

c_mu = c# 2.8   #3.0 #7.2 #1.2 #7.2 #3.4 #3.4 #5.60
c_sigma = c#3.2      #7.2  #1.2 # 7.2 #3.60 #5.40
c_xi = c# 3.0   #7.2 #1.2 #7.2 #2.4 #3.40 #5.20
#m_1 = [20, 20]
m_mu = m #[12, 11, 20]  #[12, 11, 22] #[20, 20, 20] 
m_sigma = m#[12, 11, 20]    #[12, 11, 22] #[20, 20, 20] 
m_xi =m#[12, 11, 20]    #[12, 11,22]#[20, 20, 20] # [30,30,30] #[20, 20, 20]
#m_1 = [40, 40]
#m_1 = [200,200]
#m_1 = [400, 400]
#m_1 = [100,100]

gev_var_names = ["xi_i", "mu_i", "sigma_i"]

model_C = "Common"
model_H = "Hilbert"
#model_K = "Kernel"
model_L = "Latent"
model_S = "Separate"

with pm.Model(coords=coord_train) as model_C:
    #Observations data
    sites =pm.Data("sites", train_sites)
    #Hyperparameters
    mu_mu = pm.Normal("mu_mu",mu=.7, sigma=.2)
    sigma_mu = pm.HalfNormal("sigma_mu", sigma=.25)
    #Sigma
    mu_sigma = pm.TruncatedNormal("mu_sigma", mu=.18, sigma=0.015, lower=.01, upper=.25) 
    sigma_sigma= pm.HalfNormal("sigma_sigma", sigma=.2)
    #Xi
    #mu_xi = pm.Normal("mu_xi", mu=-.3, sigma=.2)
    sigma_xi = pm.HalfNormal("sigma_xi", sigma=1)

#   Site-specific parameters of the GEV
    mu_i = pm.Normal("mu_i",mu=mu_mu, sigma=sigma_mu, shape=sites.shape ,dims="station")
    sigma_i = pm.TruncatedNormal('sigma_i',mu=mu_sigma, sigma=sigma_sigma, lower=0.01, upper=.3, shape=sites.shape, dims="station")
    xi_i = pm.TruncatedNormal('xi_i', mu=0, sigma=sigma_xi,lower=-0.4, upper=0.1,  shape=sites.shape,dims="station")


    #The GEV for each site
    Ys = pmx.GenExtreme("Ys", mu=mu_i, sigma=sigma_i, xi=xi_i, scipy=False, observed=data_train, dims=("years","station"))
    #Return level
    Zs_5 = pm.Deterministic("Zs_5", mu_i- sigma_i / xi_i *  (1 -(-np.log(1-p5))**(-xi_i)), dims= "station")
    Zs_10 = pm.Deterministic("Zs_10", mu_i - sigma_i / xi_i *  (1 -(-np.log(1-p10))**(-xi_i)), dims="station")
    Zs_20 = pm.Deterministic("Zs_20", mu_i - sigma_i / xi_i *  (1 -(-np.log(1-p20))**(-xi_i)), dims="station")
    Zs_30 = pm.Deterministic("Zs_30", mu_i - sigma_i / xi_i *  (1 -(-np.log(1-p30))**(-xi_i)), dims="station")
    Zs_50 = pm.Deterministic("Zs_50", mu_i - sigma_i / xi_i *  (1 -(-np.log(1-p50))**(-xi_i)), dims="station")
    Zs_100 = pm.Deterministic("Zs_100", mu_i - sigma_i / xi_i * (1 -(-np.log(1-p100))**(-xi_i)), dims="station")
    Zs_200 = pm.Deterministic("Zs_200", mu_i - sigma_i / xi_i *  (1 -(-np.log(1-p200))**(-xi_i)), dims="station")
    Zs_500 = pm.Deterministic("Zs_500", mu_i - sigma_i / xi_i * (1 -(-np.log(1-p500))**(-xi_i)), dims="station")
    Zs_1000 = pm.Deterministic("Zs_1000", mu_i - sigma_i / xi_i * (1 -(-np.log(1-p1000))**(-xi_i)), dims="station")
    Zs_5000 = pm.Deterministic("Zs_5000", mu_i - sigma_i / xi_i * (1 -(-np.log(1-p5000))**(-xi_i)), dims="station")
    Zs_10000 = pm.Deterministic("Zs_10000", mu_i - sigma_i / xi_i * (1 -(-np.log(1-p10000))**(-xi_i)), dims="station")
    Zs_50000 = pm.Deterministic("Zs_50000", mu_i - sigma_i / xi_i * (1 -(-np.log(1-p50000))**(-xi_i)), dims="station")
    Zs_100000 = pm.Deterministic("Zs_100000", mu_i - sigma_i / xi_i * (1 -(-np.log(1-p100000))**(-xi_i)), dims="station")

/home/sm_seuba/.conda/envs/pymc_env/lib/python3.11/site-packages/pymc/model/core.py:1359: ImputationWarning:

Data in Ys contains missing values and will be automatically imputed from the sampling distribution.

with pm.Model(coords=coord_train) as model_H:
    #XY tuples needed for the covariance 
    sites =pm.Data("sites", train_sites)
    _XY = pm.Data("XY_train", XY_train, dims=("station", "feature"))    
    #HSGP Mu
    eta_mu = pm.HalfNormal("eta_mu", sigma=0.5)#sigma=0.05)
    cov_mu = eta_mu**2* Matern52Chordal(2, ls=ls_mu)
    beta_mu_inter = pm.Normal("beta_mu_0", mu=0.45, sigma=.1)
    beta_mu = pm.Normal("beta_mu", mu=.001, sigma=0.0025, shape=2)
    mean_mu =pm.gp.mean.Linear(coeffs=beta_mu, intercept=beta_mu_inter)
    hsgp_mu = HSGP_LatLon(m_mu,  c=c_mu, ls=ls_mu, eps=eta_mu, cov_func_name="Matern52", cov_func=cov_mu, mean_func=mean_mu)
    f_mu = hsgp_mu.prior("f_mu", X=_XY, dims="station")
    #HSGP Sigma
    eta_sigma = pm.HalfNormal("eta_sigma", sigma=0.5)#sigma=0.05)
    cov_sigma = eta_sigma** 2 * Matern52Chordal(2, ls=ls_sigma)
    beta_sigma_inter =pm.TruncatedNormal("beta_sigma_0",mu=.18, sigma=.03, lower=0.1, upper=.25)  #beta_sigma_inter =pm.TruncatedNormal("beta_sigma_0",mu=.16, sigma=.03, lower=0.05, upper=.25) 
    beta_sigma = pm.TruncatedNormal("beta_sigma",mu=.00025, sigma=.0002, lower=0.001, upper=.005, shape=2) #beta_sigma = pm.TruncatedNormal("beta_sigma",mu=.00025, sigma=.0001, lower=0.001, upper=.005, shape=2) #, shape=sites.shape)
    mean_sigma = pm.gp.mean.Linear(coeffs=beta_sigma, intercept=beta_sigma_inter)
    hsgp_sigma = HSGP_LatLon(m=m_sigma,  c=c_sigma, ls=ls_sigma,eps=eta_sigma,cov_func_name="Matern52", cov_func=cov_sigma, mean_func=mean_sigma)
    f_sigma = hsgp_sigma.prior("f_sigma", X=_XY)
    #HGSP Xi
    eta_xi = pm.HalfNormal("eta_xi", sigma=0.5)#sigma=0.05)
    cov_xi = eta_xi** 2 * Matern52Chordal(2, ls=ls_xi)
    beta_xi = pm.TruncatedNormal("beta_xi", mu=-0.002, sigma=0.003, lower=-.05, upper=0.01, shape=2) #, shape=sites.shape)
    mean_xi = pm.gp.mean.Linear(coeffs=beta_xi)
    hsgp_xi = HSGP_LatLon(m=m_xi, c=c_xi, ls=ls_xi, eps=eta_xi,cov_func_name="Matern52",cov_func=cov_xi, mean_func=mean_xi)
    f_xi = hsgp_xi.prior("f_xi", X=_XY)
    #Coefficents
    eps_mu = pm.Normal("mu_eps_mu", mu=0, sigma=0.0001) #eps_mu = pm.Normal("mu_eps_mu", mu=0.05, sigma=0.01) #, shape=sites.shape)
    eps_sigma = pm.HalfNormal("eps_sigma", sigma=0.0001 ) #eps_sigma = pm.HalfNormal("eps_sigma", sigma=0.01)
    eps_xi = pm.Normal("eps_xi", mu=0, sigma=.0001) #eps_xi = pm.Normal("eps_xi", mu=0, sigma=.001)
    #GEV parameters
    xi_i= pm.Deterministic("xi_i", f_xi+eps_xi, dims="station")
    sigma_i = pm.Deterministic("sigma_i", at.abs(f_sigma)+eps_sigma, dims="station")
    mu_i = pm.Deterministic("mu_i",  at.abs(f_mu) +eps_mu, dims="station")
    #GEV
    Ys = pmx.GenExtreme("Ys", mu=mu_i, sigma=sigma_i, xi=xi_i, scipy=False, observed=data_train, dims=("years","station"))
    #Return level
    Zs_5 = pm.Deterministic("Zs_5", mu_i- sigma_i / xi_i *  (1 -(-np.log(1-p5))**(-xi_i)), dims= "station")
    Zs_10 = pm.Deterministic("Zs_10", mu_i - sigma_i / xi_i *  (1 -(-np.log(1-p10))**(-xi_i)), dims="station")
    Zs_20 = pm.Deterministic("Zs_20", mu_i - sigma_i / xi_i *  (1 -(-np.log(1-p20))**(-xi_i)), dims="station")
    Zs_30 = pm.Deterministic("Zs_30", mu_i - sigma_i / xi_i *  (1 -(-np.log(1-p30))**(-xi_i)), dims="station")
    Zs_50 = pm.Deterministic("Zs_50", mu_i - sigma_i / xi_i *  (1 -(-np.log(1-p50))**(-xi_i)), dims="station")
    Zs_100 = pm.Deterministic("Zs_100", mu_i - sigma_i / xi_i * (1 -(-np.log(1-p100))**(-xi_i)), dims="station")
    Zs_200 = pm.Deterministic("Zs_200", mu_i - sigma_i / xi_i *  (1 -(-np.log(1-p200))**(-xi_i)), dims="station")
    Zs_500 = pm.Deterministic("Zs_500", mu_i - sigma_i / xi_i * (1 -(-np.log(1-p500))**(-xi_i)), dims="station")
    Zs_1000 = pm.Deterministic("Zs_1000", mu_i - sigma_i / xi_i * (1 -(-np.log(1-p1000))**(-xi_i)), dims="station")
    Zs_5000 = pm.Deterministic("Zs_5000", mu_i - sigma_i / xi_i * (1 -(-np.log(1-p5000))**(-xi_i)), dims="station")
    Zs_10000 = pm.Deterministic("Zs_10000", mu_i - sigma_i / xi_i * (1 -(-np.log(1-p10000))**(-xi_i)), dims="station")
    Zs_50000 = pm.Deterministic("Zs_50000", mu_i - sigma_i / xi_i * (1 -(-np.log(1-p50000))**(-xi_i)), dims="station")
    Zs_100000 = pm.Deterministic("Zs_100000", mu_i - sigma_i / xi_i * (1 -(-np.log(1-p100000))**(-xi_i)), dims="station")

/home/sm_seuba/.conda/envs/pymc_env/lib/python3.11/site-packages/pymc/model/core.py:1359: ImputationWarning:

Data in Ys contains missing values and will be automatically imputed from the sampling distribution.

with pm.Model(coords=coord_train) as model_L:
    #XY tuples needed for the covariance 
    sites =pm.Data("sites", train_sites)
    _XY = pm.Data("XY_train", XY_train, dims=("station", "feature"))    
    #GP Mu
    eta_mu = pm.HalfNormal("eta_mu", sigma=0.05)
    cov_mu = eta_mu**2* Matern52Chordal(2, ls=ls_mu)
    beta_mu_inter = pm.Normal("beta_mu_0", mu=0.45, sigma=.1)
    beta_mu = pm.Normal("beta_mu", mu=.001, sigma=0.0025, shape=2)
    mean_mu =pm.gp.mean.Linear(coeffs=beta_mu, intercept=beta_mu_inter)
    gp_mu = pm.gp.Latent(mean_func=mean_mu, cov_func= cov_mu)
    f_mu = gp_mu.prior("f_mu", X=_XY)
    #GP Sigma
    eta_sigma =  pm.HalfNormal("eta_sigma", sigma=0.05)
    cov_sigma = eta_sigma** 2 * Matern52Chordal(2, ls=ls_sigma)
    beta_sigma_inter =pm.TruncatedNormal("beta_sigma_0",mu=.18, sigma=.03, lower=0.1, upper=.25) 
    beta_sigma = pm.TruncatedNormal("beta_sigma",mu=.00025, sigma=.0002, lower=0.001, upper=.005, shape=2) #, shape=sites.shape)
    mean_sigma = pm.gp.mean.Linear(coeffs=beta_sigma, intercept=beta_sigma_inter)
    gp_sigma = pm.gp.Latent(mean_func=mean_sigma, cov_func= cov_sigma)
    f_sigma = gp_sigma.prior("f_sigma", X=_XY)
    #GP Xi
    eta_xi =  pm.HalfNormal("eta_xi", sigma=0.05)
    cov_xi = eta_xi** 2 * Matern52Chordal(2, ls=ls_xi)
    beta_xi = pm.TruncatedNormal("beta_xi", mu=-0.0025, sigma=0.003, lower=-.05, upper=0.01, shape=2) #, shape=sites.shape)
    mean_xi = pm.gp.mean.Linear(coeffs=beta_xi)
    gp_xi = pm.gp.Latent(mean_func=mean_xi, cov_func= cov_xi)
    f_xi = gp_xi.prior("f_xi", X=_XY)
    #Error
    eps_mu = pm.Normal("mu_eps_mu", mu=0, sigma=0.0001) #eps_mu = pm.Normal("mu_eps_mu", mu=0.05, sigma=0.01) #, shape=sites.shape)
    eps_sigma = pm.HalfNormal("eps_sigma",sigma=0.0001 ) #eps_sigma = pm.HalfNormal("eps_sigma", sigma=0.01)
    eps_xi = pm.Normal("eps_xi", mu=0, sigma=.0001) #eps_xi = pm.Normal("eps_xi", mu=0, sigma=.001)
    #GEV parameters
    xi_i= pm.Deterministic("xi_i", f_xi+eps_xi, dims="station")
    sigma_i = pm.Deterministic("sigma_i", at.abs(f_sigma) + eps_sigma, dims="station")
    mu_i = pm.Deterministic("mu_i",  at.abs(f_mu) +eps_mu, dims="station")
    #GEV
    Ys = pmx.GenExtreme("Ys", mu=mu_i, sigma=sigma_i, xi=xi_i, scipy=False, observed=data_train, dims=("years","station"))
    #Return level
    Zs_5 = pm.Deterministic("Zs_5", mu_i- sigma_i / xi_i *  (1 -(-np.log(1-p5))**(-xi_i)), dims= "station")
    Zs_10 = pm.Deterministic("Zs_10", mu_i - sigma_i / xi_i *  (1 -(-np.log(1-p10))**(-xi_i)), dims="station")
    Zs_20 = pm.Deterministic("Zs_20", mu_i - sigma_i / xi_i *  (1 -(-np.log(1-p20))**(-xi_i)), dims="station")
    Zs_30 = pm.Deterministic("Zs_30", mu_i - sigma_i / xi_i *  (1 -(-np.log(1-p30))**(-xi_i)), dims="station")
    Zs_50 = pm.Deterministic("Zs_50", mu_i - sigma_i / xi_i *  (1 -(-np.log(1-p50))**(-xi_i)), dims="station")
    Zs_100 = pm.Deterministic("Zs_100", mu_i - sigma_i / xi_i * (1 -(-np.log(1-p100))**(-xi_i)), dims="station")
    Zs_200 = pm.Deterministic("Zs_200", mu_i - sigma_i / xi_i *  (1 -(-np.log(1-p200))**(-xi_i)), dims="station")
    Zs_500 = pm.Deterministic("Zs_500", mu_i - sigma_i / xi_i * (1 -(-np.log(1-p500))**(-xi_i)), dims="station")
    Zs_1000 = pm.Deterministic("Zs_1000", mu_i - sigma_i / xi_i * (1 -(-np.log(1-p1000))**(-xi_i)), dims="station")
    Zs_5000 = pm.Deterministic("Zs_5000", mu_i - sigma_i / xi_i * (1 -(-np.log(1-p5000))**(-xi_i)), dims="station")
    Zs_10000 = pm.Deterministic("Zs_10000", mu_i - sigma_i / xi_i * (1 -(-np.log(1-p10000))**(-xi_i)), dims="station")
    Zs_50000 = pm.Deterministic("Zs_50000", mu_i - sigma_i / xi_i * (1 -(-np.log(1-p50000))**(-xi_i)), dims="station")
    Zs_100000 = pm.Deterministic("Zs_100000", mu_i - sigma_i / xi_i * (1 -(-np.log(1-p100000))**(-xi_i)), dims="station")

/home/sm_seuba/.conda/envs/pymc_env/lib/python3.11/site-packages/pymc/model/core.py:1359: ImputationWarning:

Data in Ys contains missing values and will be automatically imputed from the sampling distribution.

with pm.Model(coords=coord_train) as model_S:
    s = pm.Data("s", train_sites)
    #y_obs = pm.MutableData("y_obs", BM, dims=("years","station"))
    #GEV paramters
    mu_i = pm.Normal("mu_i",mu=.74, sigma=0.3,dims="station")
    sigma_i= pm.TruncatedNormal("sigma_i", mu=.18, sigma=0.2, lower=.1, upper=.25,dims="station")
    xi_i = pm.TruncatedNormal("xi_i", mu=0, sigma=.5, lower=-.4, upper=.1, dims="station")
    #GEV
    Ys = pmx.GenExtreme("Ys", mu=mu_i, sigma=sigma_i, xi=xi_i, scipy=False, observed=data_train, dims=("years","station"))
    #Return level
    Zs_5 = pm.Deterministic("Zs_5", mu_i- sigma_i / xi_i *  (1 -(-np.log(1-p5))**(-xi_i)), dims= "station")
    Zs_10 = pm.Deterministic("Zs_10", mu_i - sigma_i / xi_i *  (1 -(-np.log(1-p10))**(-xi_i)), dims="station")
    Zs_20 = pm.Deterministic("Zs_20", mu_i - sigma_i / xi_i *  (1 -(-np.log(1-p20))**(-xi_i)), dims="station")
    Zs_30 = pm.Deterministic("Zs_30", mu_i - sigma_i / xi_i *  (1 -(-np.log(1-p30))**(-xi_i)), dims="station")
    Zs_50 = pm.Deterministic("Zs_50", mu_i - sigma_i / xi_i *  (1 -(-np.log(1-p50))**(-xi_i)), dims="station")
    Zs_100 = pm.Deterministic("Zs_100", mu_i - sigma_i / xi_i * (1 -(-np.log(1-p100))**(-xi_i)), dims="station")
    Zs_200 = pm.Deterministic("Zs_200", mu_i - sigma_i / xi_i *  (1 -(-np.log(1-p200))**(-xi_i)), dims="station")
    Zs_500 = pm.Deterministic("Zs_500", mu_i - sigma_i / xi_i * (1 -(-np.log(1-p500))**(-xi_i)), dims="station")
    Zs_1000 = pm.Deterministic("Zs_1000", mu_i - sigma_i / xi_i * (1 -(-np.log(1-p1000))**(-xi_i)), dims="station")
    Zs_5000 = pm.Deterministic("Zs_5000", mu_i - sigma_i / xi_i * (1 -(-np.log(1-p5000))**(-xi_i)), dims="station")
    Zs_10000 = pm.Deterministic("Zs_10000", mu_i - sigma_i / xi_i * (1 -(-np.log(1-p10000))**(-xi_i)), dims="station")
    Zs_50000 = pm.Deterministic("Zs_50000", mu_i - sigma_i / xi_i * (1 -(-np.log(1-p50000))**(-xi_i)), dims="station")
    Zs_100000 = pm.Deterministic("Zs_100000", mu_i - sigma_i / xi_i * (1 -(-np.log(1-p100000))**(-xi_i)), dims="station")

/home/sm_seuba/.conda/envs/pymc_env/lib/python3.11/site-packages/pymc/model/core.py:1359: ImputationWarning:

Data in Ys contains missing values and will be automatically imputed from the sampling distribution.

pm.model_to_graphviz(model_C)

pm.model_to_graphviz(model_H)

pm.model_to_graphviz(model_L)

pm.model_to_graphviz(model_S)

with model_C:
    idataC = pm.sample_prior_predictive(samples=10000, random_seed=rng)

Sampling: [Ys_observed, Ys_unobserved, mu_i, mu_mu, mu_sigma, sigma_i, sigma_mu, sigma_sigma, sigma_xi, xi_i]

with model_H:
    idataH = pm.sample_prior_predictive(samples=10000, random_seed=rng)

Sampling: [Ys_observed, Ys_unobserved, beta_mu, beta_mu_0, beta_sigma, beta_sigma_0, beta_xi, eps_sigma, eps_xi, eta_mu, eta_sigma, eta_xi, f_mu_hsgp_coeffs_, f_sigma_hsgp_coeffs_, f_xi_hsgp_coeffs_, mu_eps_mu]

with model_L:
    idataL = pm.sample_prior_predictive(samples=10000, random_seed=rng)

Sampling: [Ys_observed, Ys_unobserved, beta_mu, beta_mu_0, beta_sigma, beta_sigma_0, beta_xi, eps_sigma, eps_xi, eta_mu, eta_sigma, eta_xi, f_mu_rotated_, f_sigma_rotated_, f_xi_rotated_, mu_eps_mu]

with model_S:
    idataS = pm.sample_prior_predictive(samples=10000, random_seed=rng)

Sampling: [Ys_observed, Ys_unobserved, mu_i, sigma_i, xi_i]

all_idata = [idataC, idataH, idataL, idataS]
models_names = ["Common", "Hilbert", "Latent", "Separate"]

fig, axes = plt.subplots(nrows=4, ncols=2, figsize=(20,20))
row_titles = ["PPC plot (KDE)", "PPC plot (ECDF)"]
for row, idata in enumerate(all_idata):
    az.plot_ppc(idata, observed=True, observed_rug=True, num_pp_samples=5000, group="prior", figsize=(12, 6), ax=axes[row,0]);
    az.plot_ppc(idata, observed=True, observed_rug=True,  kind='cumulative',num_pp_samples=5000, group="prior", figsize=(12, 6), ax=axes[row,1]);
    axes[row,0].set_ylabel(models_names[row], fontsize=20)
for col, row_title in enumerate(row_titles):
    axes[0,col].set_title(row_title,  fontdict={"fontsize":20})
fig.suptitle("Visualization of Prior Predictive Checks", fontsize=30);
fig.show();

with model_C:
    steps = pm.NUTS(target_accept=0.995)
    traceC =pm.sample(draws=10000, tune=5000, chains=4, cores=4, random_seed=rng, step=steps,initvals={ "xi_i":np.full(dim_train[1],0)}, return_inferencedata=True)
    #traceC =pm.sample(draws=5000, tune=5000, chains=4, cores=4, random_seed=rng, step=steps,initvals={ "xi_i":np.full(dim_train[1],0)}, return_inferencedata=True)
    #traceC =pm.sample(draws=3000, tune=3000, chains=4, cores=4, random_seed=rng, step=steps,initvals={ "xi_i":np.full(dim_train[1],0)}, return_inferencedata=True)
    #traceC =pm.sample(draws=1000, tune=1000, chains=4, cores=4, random_seed=rng, step=steps,initvals={ "xi_i":np.full(dim_train[1],0)}, return_inferencedata=True)

Multiprocess sampling (4 chains in 4 jobs)

NUTS: [mu_mu, sigma_mu, mu_sigma, sigma_sigma, sigma_xi, mu_i, sigma_i, xi_i, Ys_unobserved]

with model_H:
    steps = pm.NUTS(target_accept=0.99) #target_accept=0.98) 
    traceH =pm.sample(draws=10000, tune=5000, chains=4, cores=4, random_seed=rng, step=steps,initvals={"beta_xi":np.full(2,0), "eps_xi":0}, return_inferencedata=True)
    #traceH =pm.sample(draws=5000, tune=5000, chains=4, cores=4, random_seed=rng, step=steps,initvals={"beta_xi":np.full(2,0), "eps_xi":0}, return_inferencedata=True)
    #traceH =pm.sample(draws=3000, tune=3000, chains=4, cores=4, random_seed=rng, step=steps,initvals={"beta_xi":np.full(2,0), "eps_xi":0}, return_inferencedata=True)
    #traceH =pm.sample(draws=1000, tune=1000, chains=4, cores=4, random_seed=rng, step=steps,initvals={"beta_xi":np.full(2,0), "eps_xi":0}, return_inferencedata=True)

/home/sm_seuba/.conda/envs/pymc_env/lib/python3.11/site-packages/pytensor/tensor/elemwise.py:754: RuntimeWarning:

invalid value encountered in power

/home/sm_seuba/.conda/envs/pymc_env/lib/python3.11/site-packages/pytensor/tensor/elemwise.py:754: RuntimeWarning:

invalid value encountered in power

/home/sm_seuba/.conda/envs/pymc_env/lib/python3.11/site-packages/pytensor/tensor/elemwise.py:754: RuntimeWarning:

invalid value encountered in power

/home/sm_seuba/.conda/envs/pymc_env/lib/python3.11/site-packages/pytensor/tensor/elemwise.py:754: RuntimeWarning:

invalid value encountered in power

/home/sm_seuba/.conda/envs/pymc_env/lib/python3.11/site-packages/pytensor/tensor/elemwise.py:754: RuntimeWarning:

invalid value encountered in power

/home/sm_seuba/.conda/envs/pymc_env/lib/python3.11/site-packages/pytensor/tensor/elemwise.py:754: RuntimeWarning:

invalid value encountered in power

Multiprocess sampling (4 chains in 4 jobs)

NUTS: [eta_mu, beta_mu_0, beta_mu, f_mu_hsgp_coeffs_, eta_sigma, beta_sigma_0, beta_sigma, f_sigma_hsgp_coeffs_, eta_xi, beta_xi, f_xi_hsgp_coeffs_, mu_eps_mu, eps_sigma, eps_xi, Ys_unobserved]

with model_L:
    steps = pm.NUTS(target_accept=0.99) #0.98
    traceL =pm.sample(draws=10000, tune=5000, chains=4, cores=4, random_seed=rng, step=steps,initvals={"beta_xi":np.full(2,0), "eps_xi":0}, return_inferencedata=True)
    #traceL =pm.sample(draws=5000, tune=5000, chains=4, cores=4, random_seed=rng, step=steps,initvals={"beta_xi":np.full(2,0), "eps_xi":0}, return_inferencedata=True)
    #traceL =pm.sample(draws=3000, tune=3000, chains=4, cores=4, random_seed=rng, step=steps,initvals={"beta_xi":np.full(2,0), "eps_xi":0}, return_inferencedata=True)
    #traceL =pm.sample(draws=1000, tune=1000, chains=4, cores=4, random_seed=rng, step=steps,initvals={"beta_xi":np.full(2,0), "eps_xi":0}, return_inferencedata=True)

/home/sm_seuba/.conda/envs/pymc_env/lib/python3.11/site-packages/pytensor/tensor/elemwise.py:754: RuntimeWarning:

invalid value encountered in power

/home/sm_seuba/.conda/envs/pymc_env/lib/python3.11/site-packages/pytensor/tensor/elemwise.py:754: RuntimeWarning:

invalid value encountered in power

/home/sm_seuba/.conda/envs/pymc_env/lib/python3.11/site-packages/pytensor/tensor/elemwise.py:754: RuntimeWarning:

invalid value encountered in power

/home/sm_seuba/.conda/envs/pymc_env/lib/python3.11/site-packages/pytensor/tensor/elemwise.py:754: RuntimeWarning:

invalid value encountered in power

/home/sm_seuba/.conda/envs/pymc_env/lib/python3.11/site-packages/pytensor/tensor/elemwise.py:754: RuntimeWarning:

invalid value encountered in power

/home/sm_seuba/.conda/envs/pymc_env/lib/python3.11/site-packages/pytensor/tensor/elemwise.py:754: RuntimeWarning:

invalid value encountered in power

Multiprocess sampling (4 chains in 4 jobs)

NUTS: [eta_mu, beta_mu_0, beta_mu, f_mu_rotated_, eta_sigma, beta_sigma_0, beta_sigma, f_sigma_rotated_, eta_xi, beta_xi, f_xi_rotated_, mu_eps_mu, eps_sigma, eps_xi, Ys_unobserved]

with model_S:
    steps = pm.NUTS(target_accept=0.995)
    traceS =pm.sample(draws=10000, tune=5000, chains=4, cores=4, random_seed=rng, step=steps, initvals={ "xi_i":np.full(dim_train[1],0)},return_inferencedata=True)
    #traceS =pm.sample(draws=5000, tune=5000, chains=4, cores=4, random_seed=rng, step=steps, initvals={ "xi_i":np.full(dim_train[1],0)},return_inferencedata=True)
    #traceS =pm.sample(draws=3000, tune=3000, chains=4, cores=4, random_seed=rng, step=steps, initvals={ "xi_i":np.full(dim_train[1],0)},return_inferencedata=True)
    #traceS =pm.sample(draws=1000, tune=1000, chains=4, cores=4, random_seed=rng, step=steps, initvals={ "xi_i":np.full(dim_train[1],0)},return_inferencedata=True)

Multiprocess sampling (4 chains in 4 jobs)

NUTS: [mu_i, sigma_i, xi_i, Ys_unobserved]

with model_C:
    pm.compute_log_likelihood(traceC)

with model_H:
    pm.compute_log_likelihood(traceH)

with model_L:
    pm.compute_log_likelihood(traceL)

with model_S:
    pm.compute_log_likelihood(traceS)

fig, axes = plt.subplots(nrows=1, ncols=4, figsize=(15,5))
all_trace = [traceC, traceH, traceL, traceS]
for col, trace in enumerate(all_trace):
    az.plot_energy(trace, ax=axes[col]);
    axes[col].set_title(models_names[col],  fontdict={"fontsize":20});
fig.suptitle("Energy transition distribution and marginal energy distribution plots", fontsize=30);
fig.tight_layout();
fig.show();

looC, looH,  looL, looS = az.loo(traceC),az.loo(traceH),  az.loo(traceL), az.loo(traceS)

/home/sm_seuba/.conda/envs/pymc_env/lib/python3.11/site-packages/arviz/stats/stats.py:795: UserWarning:

Estimated shape parameter of Pareto distribution is greater than 0.70 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.

/home/sm_seuba/.conda/envs/pymc_env/lib/python3.11/site-packages/arviz/stats/stats.py:795: UserWarning:

Estimated shape parameter of Pareto distribution is greater than 0.70 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.

looC

Computed from 40000 posterior samples and 1980 observations log-likelihood matrix.

         Estimate       SE
elpd_loo   310.59    30.69
p_loo       78.15        -

There has been a warning during the calculation. Please check the results.
------

Pareto k diagnostic values:
                         Count   Pct.
(-Inf, 0.70]   (good)     1976   99.8%
   (0.70, 1]   (bad)         4    0.2%
   (1, Inf)   (very bad)    0    0.0%

looH

Computed from 40000 posterior samples and 1980 observations log-likelihood matrix.

         Estimate       SE
elpd_loo   319.04    30.05
p_loo       39.27        -
------

Pareto k diagnostic values:
                         Count   Pct.
(-Inf, 0.70]   (good)     1980  100.0%
   (0.70, 1]   (bad)         0    0.0%
   (1, Inf)   (very bad)    0    0.0%

looL

Computed from 40000 posterior samples and 1980 observations log-likelihood matrix.

         Estimate       SE
elpd_loo   318.17    30.20
p_loo       35.40        -
------

Pareto k diagnostic values:
                         Count   Pct.
(-Inf, 0.70]   (good)     1980  100.0%
   (0.70, 1]   (bad)         0    0.0%
   (1, Inf)   (very bad)    0    0.0%

looS

Computed from 40000 posterior samples and 1980 observations log-likelihood matrix.

         Estimate       SE
elpd_loo   309.92    31.55
p_loo       80.55        -

There has been a warning during the calculation. Please check the results.
------

Pareto k diagnostic values:
                         Count   Pct.
(-Inf, 0.70]   (good)     1976   99.8%
   (0.70, 1]   (bad)         4    0.2%
   (1, Inf)   (very bad)    0    0.0%

fig, axes = plt.subplots(nrows=1, ncols=4, figsize=(20,7))
all_the_loo = [looC, looH, looL, looS]
for col, loo in enumerate(all_the_loo):
    az.plot_khat(loo, show_bins=True, threshold=.7, hover_label=True, ax=axes[col]);
    axes[col].axhline(0.70, color="red",  linestyle="dashed",label="0.7");
    axes[col].axhline(0.50, color="black",  linestyle="dashed",label="0.5");
    axes[col].set_title(models_names[col],  fontdict={"fontsize":20});
fig.suptitle("Plot Pareto tail indices for diagnosing convergence in PSIS-LOO.", fontsize=30);
fig.tight_layout();
fig.show();

/home/sm_seuba/.conda/envs/pymc_env/lib/python3.11/site-packages/arviz/plots/backends/matplotlib/khatplot.py:44: UserWarning:

hover labels are only available with interactive backends. To switch to an interactive backend from ipython or jupyter, use `%matplotlib` there should be no need to restart the kernel. For other cases, see https://matplotlib.org/3.1.0/tutorials/introductory/usage.html#backends

fig, axes = plt.subplots(nrows=4, ncols=3, figsize=(30,20))
all_trace = [traceC, traceH, traceL, traceS]
var_names = ["xi_i", "mu_i", "sigma_i"]
for row, trace in enumerate(all_trace):
    for col,var in enumerate(var_names):
        if row == 0:
            axes[row, col].set_title(var,  fontdict={"fontsize":20})
        az.plot_parallel(trace, var_names=var, ax=axes[row, col]);
    axes[row,0].set_ylabel(models_names[row], fontsize=20)
fig.suptitle("Plot parallel coordinates plot showing posterior points with and without divergences.", fontsize=30)
fig.tight_layout();
fig.show();

/tmp/ipykernel_802820/3784022070.py:11: UserWarning:

Creating legend with loc="best" can be slow with large amounts of data.

/home/sm_seuba/.conda/envs/pymc_env/lib/python3.11/site-packages/IPython/core/events.py:93: UserWarning:

Creating legend with loc="best" can be slow with large amounts of data.

/home/sm_seuba/.conda/envs/pymc_env/lib/python3.11/site-packages/IPython/core/pylabtools.py:152: UserWarning:

Creating legend with loc="best" can be slow with large amounts of data.

az.plot_rank(traceC, var_names="xi_i",figsize=(29,25));
plt.suptitle("Common rank plot for the GEV shape param", fontsize=20);

az.plot_rank(traceC, var_names="mu_i",figsize=(29,25));
plt.suptitle("Common rank plot for the GEV location param", fontsize=20);

az.plot_rank(traceC, var_names="sigma_i",figsize=(29,25));
plt.suptitle("Common rank plot for the GEV scale param", fontsize=20);

az.plot_rank(traceH, var_names="xi_i",figsize=(29,25));
plt.suptitle("Hilbert rank plot for the GEV shape param", fontsize=20);

az.plot_rank(traceH, var_names="mu_i",figsize=(29,25));
plt.suptitle("Hilbert rank plot for the GEV location param", fontsize=20);

az.plot_rank(traceH, var_names="sigma_i",figsize=(29,25));
plt.suptitle("Hilbert rank plot for the GEV scale param", fontsize=20);

az.plot_rank(traceL, var_names="xi_i",figsize=(29,25));
plt.suptitle("Latent rank plot for the GEV shape param", fontsize=20);

az.plot_rank(traceL, var_names="mu_i",figsize=(29,25));
plt.suptitle("Latent rank plot for the GEV location param", fontsize=20);

az.plot_rank(traceL, var_names="sigma_i",figsize=(29,25));
plt.suptitle("Latent rank plot for the GEV scale param", fontsize=20);

az.plot_rank(traceS, var_names="xi_i",figsize=(29,25));
plt.suptitle("Separate rank plot for the GEV shape param", fontsize=20);

az.plot_rank(traceS, var_names="mu_i",figsize=(29,25));
plt.suptitle("Separate rank plot for the GEV location param", fontsize=20);

az.plot_rank(traceS, var_names="sigma_i",figsize=(29,25));
plt.suptitle("Separate rank plot for the GEV scale param", fontsize=20);

az.plot_ess(traceC, kind="local", var_names=["xi_i"],  figsize=(29,25),drawstyle="steps-mid", color="k",
    linestyle="-", marker=None, rug=True,  rug_kwargs={"color": "r"});
plt.suptitle("Common local ESS plot for the GEV shape param", fontsize=20);

az.plot_ess(traceC, kind="local", var_names=["mu_i"],  figsize=(29,25),drawstyle="steps-mid", color="k",
    linestyle="-", marker=None, rug=True,  rug_kwargs={"color": "r"});
plt.suptitle("Common local ESS plot for the GEV location param", fontsize=20);

az.plot_ess(traceC, kind="local", var_names=["sigma_i"],  figsize=(29,25),drawstyle="steps-mid", color="k",
    linestyle="-", marker=None, rug=True,  rug_kwargs={"color": "r"});
plt.suptitle("Common local ESS plot for the GEV scale param", fontsize=20);

az.plot_ess(traceH, kind="local", var_names=["xi_i"],  figsize=(29,25),drawstyle="steps-mid", color="k",
    linestyle="-", marker=None, rug=True,  rug_kwargs={"color": "r"});
plt.suptitle("Hilbert local ESS plot for the GEV shape param", fontsize=20);

az.plot_ess(traceH, kind="local", var_names=["mu_i"],  figsize=(29,25),drawstyle="steps-mid", color="k",
    linestyle="-", marker=None, rug=True,  rug_kwargs={"color": "r"});
plt.suptitle("Hilbert local ESS plot for the GEV location param", fontsize=20);

az.plot_ess(traceH, kind="local", var_names=["sigma_i"],  figsize=(29,25),drawstyle="steps-mid", color="k",
    linestyle="-", marker=None, rug=True,  rug_kwargs={"color": "r"});
plt.suptitle("Hilbert local ESS plot for the GEV scale param", fontsize=20);

az.plot_ess(traceL, kind="local", var_names=["xi_i"],  figsize=(29,25),drawstyle="steps-mid", color="k",
    linestyle="-", marker=None, rug=True,  rug_kwargs={"color": "r"});
plt.suptitle("Latent local ESS plot for the GEV shape param", fontsize=20);

az.plot_ess(traceL, kind="local", var_names=["mu_i"],  figsize=(29,25),drawstyle="steps-mid", color="k",
    linestyle="-", marker=None, rug=True,  rug_kwargs={"color": "r"});
plt.suptitle("Latent local ESS plot for the GEV location param", fontsize=20);

az.plot_ess(traceL, kind="local", var_names=["sigma_i"],  figsize=(29,25),drawstyle="steps-mid", color="k",
    linestyle="-", marker=None, rug=True,  rug_kwargs={"color": "r"});
plt.suptitle("Latent local ESS plot for the GEV scale param", fontsize=20);

az.plot_ess(traceS, kind="local", var_names=["xi_i"],  figsize=(29,25),drawstyle="steps-mid", color="k",
    linestyle="-", marker=None, rug=True,  rug_kwargs={"color": "r"});
plt.suptitle("Separate local ESS plot for the GEV shape param", fontsize=20);

az.plot_ess(traceS, kind="local", var_names=["mu_i"],  figsize=(29,25),drawstyle="steps-mid", color="k",
    linestyle="-", marker=None, rug=True,  rug_kwargs={"color": "r"});
plt.suptitle("Separate local ESS plot for the GEV location param", fontsize=20);

az.plot_ess(traceS, kind="local", var_names=["sigma_i"],  figsize=(29,25),drawstyle="steps-mid", color="k",
    linestyle="-", marker=None, rug=True,  rug_kwargs={"color": "r"});
plt.suptitle("Separate local ESS plot for the GEV scale param", fontsize=20);

az.plot_ess(traceC, kind="quantile", var_names=["xi_i"], figsize=(29,25), drawstyle="steps-mid", color="k",
    linestyle="-", marker=None, rug=True, rug_kwargs={"color": "r"});
plt.suptitle("Common quantile ESS plot for the GEV shape param", fontsize=20);

az.plot_ess(traceC, kind="quantile", var_names=["mu_i"], figsize=(29,25), drawstyle="steps-mid", color="k",
    linestyle="-", marker=None, rug=True, rug_kwargs={"color": "r"});
plt.suptitle("Common quantile ESS plot for the GEV location param", fontsize=20);

az.plot_ess(traceC, kind="quantile", var_names=["sigma_i"], figsize=(29,25), drawstyle="steps-mid", color="k",
    linestyle="-", marker=None, rug=True, rug_kwargs={"color": "r"});
plt.suptitle("Common quantile ESS plot for the GEV scale param", fontsize=20);

az.plot_ess(traceH, kind="quantile", var_names=["xi_i"], figsize=(29,25), drawstyle="steps-mid", color="k",
    linestyle="-", marker=None, rug=True, rug_kwargs={"color": "r"});
plt.suptitle("Hilbert quantile ESS plot for the GEV shape param", fontsize=20);

az.plot_ess(traceH, kind="quantile", var_names=["mu_i"], figsize=(29,25), drawstyle="steps-mid", color="k",
    linestyle="-", marker=None, rug=True, rug_kwargs={"color": "r"});
plt.suptitle("Hilbert quantile ESS plot for the GEV location param", fontsize=20);

az.plot_ess(traceH, kind="quantile", var_names=["sigma_i"], figsize=(29,25), drawstyle="steps-mid", color="k",
    linestyle="-", marker=None, rug=True, rug_kwargs={"color": "r"});
plt.suptitle("Hilbert quantile ESS plot for the GEV scale param", fontsize=20);

az.plot_ess(traceL, kind="quantile", var_names=["xi_i"], figsize=(29,25), drawstyle="steps-mid", color="k",
    linestyle="-", marker=None, rug=True, rug_kwargs={"color": "r"});
plt.suptitle("Latent quantile ESS plot for the GEV shape param", fontsize=20);

az.plot_ess(traceL, kind="quantile", var_names=["mu_i"], figsize=(29,25), drawstyle="steps-mid", color="k",
    linestyle="-", marker=None, rug=True, rug_kwargs={"color": "r"});
plt.suptitle("Latent quantile ESS plot for the GEV location param", fontsize=20);

az.plot_ess(traceL, kind="quantile", var_names=["sigma_i"], figsize=(29,25), drawstyle="steps-mid", color="k",
    linestyle="-", marker=None, rug=True, rug_kwargs={"color": "r"});
plt.suptitle("Latent quantile ESS plot for the GEV scale param", fontsize=20);

az.plot_ess(traceS, kind="quantile", var_names=["xi_i"], figsize=(20,20), drawstyle="steps-mid", color="k",
    linestyle="-", marker=None, rug=True, rug_kwargs={"color": "r"});
plt.suptitle("Separate quantile ESS plot for the GEV shape param", fontsize=20);

az.plot_ess(traceS, kind="quantile", var_names=["mu_i"], figsize=(29,25), drawstyle="steps-mid", color="k",
    linestyle="-", marker=None, rug=True, rug_kwargs={"color": "r"});
plt.suptitle("Separate quantile ESS plot for the GEV location param", fontsize=20);

az.plot_ess(traceS, kind="quantile", var_names=["sigma_i"], figsize=(29,25), drawstyle="steps-mid", color="k",
    linestyle="-", marker=None, rug=True, rug_kwargs={"color": "r"});
plt.suptitle("Separate quantile ESS plot for the GEV scale param", fontsize=20);

az.plot_ess(traceC, kind="evolution", var_names=["xi_i"], figsize=(20,20));
plt.suptitle("Common evolution ESS plot for the GEV shape param", fontsize=20);

az.plot_ess(traceC, kind="evolution", var_names=["mu_i"], figsize=(29,25));
plt.suptitle("Common evolution ESS plot for the GEV location param", fontsize=20);

az.plot_ess(traceC, kind="evolution", var_names=["sigma_i"], figsize=(20,20));
plt.suptitle("Common evolution ESS plot for the GEV scale param", fontsize=20);

az.plot_ess(traceH, kind="evolution", var_names=["xi_i"], figsize=(29,25));
plt.suptitle("Hilbert evolution ESS plot for the GEV shape param", fontsize=20);

az.plot_ess(traceH, kind="evolution", var_names=["mu_i"], figsize=(20,20));
plt.suptitle("Hilbert evolution ESS plot for the GEV location param", fontsize=20);

az.plot_ess(traceH, kind="evolution", var_names=["sigma_i"], figsize=(29,25));
plt.suptitle("Hilbert evolution ESS plot for the GEV scale param", fontsize=20);

az.plot_ess(traceL, kind="evolution", var_names=["xi_i"], figsize=(20,20));
plt.suptitle("Latent evolution ESS plot for the GEV shape param", fontsize=20);

az.plot_ess(traceL, kind="evolution", var_names=["mu_i"], figsize=(29,25));
plt.suptitle("Latent evolution ESS plot for the GEV location param", fontsize=20);

az.plot_ess(traceL, kind="evolution", var_names=["sigma_i"], figsize=(20,20));
plt.suptitle("Latent evolution ESS plot for the GEV scale param", fontsize=20);

az.plot_ess(traceS, kind="evolution", var_names=["xi_i"], figsize=(29,25));
plt.suptitle("Separate evolution ESS plot for the GEV shape param", fontsize=20);

az.plot_ess(traceS, kind="evolution", var_names=["mu_i"], figsize=(20,20));
plt.suptitle("Separate evolution ESS plot for the GEV location param", fontsize=20);

az.plot_ess(traceS, kind="evolution", var_names=["sigma_i"], figsize=(29,25));
plt.suptitle("Separate evolution ESS plot for the GEV scale param", fontsize=20);

az.plot_autocorr(traceC, var_names="xi_i", max_lag=20, combined=True, figsize=(20,20));
plt.suptitle("Common ACF plots for the GEV shape param", fontsize=20);

az.plot_autocorr(traceC, var_names="mu_i", max_lag=20, combined=True, figsize=(29,25));
plt.suptitle("Common ACF plots for the GEV location param", fontsize=20);

az.plot_autocorr(traceC, var_names="sigma_i", max_lag=20, combined=True, figsize=(20,20));
plt.suptitle("Common ACF plots for the GEV scale param", fontsize=20);

az.plot_autocorr(traceH, var_names="xi_i", max_lag=20, combined=True, figsize=(29,25));
plt.suptitle("Hilbert ACF plots for the GEV shape param", fontsize=20);

az.plot_autocorr(traceH, var_names="mu_i", max_lag=20, combined=True, figsize=(20,20));
plt.suptitle("Hilbert ACF plots for the GEV location param", fontsize=20);

az.plot_autocorr(traceH, var_names="sigma_i", max_lag=20, combined=True, figsize=(29,25));
plt.suptitle("Hilbert ACF plots for the GEV scale param", fontsize=20);

az.plot_autocorr(traceL, var_names="xi_i", max_lag=20, combined=True, figsize=(20,20));
plt.suptitle("Latent ACF plots for the GEV shape param", fontsize=20);

az.plot_autocorr(traceL, var_names="mu_i", max_lag=20, combined=True, figsize=(29,25));
plt.suptitle("Latent ACF plots for the GEV location param", fontsize=20);

az.plot_autocorr(traceL, var_names="sigma_i", max_lag=20, combined=True, figsize=(20,20));
plt.suptitle("Latent ACF plots for the GEV scale param", fontsize=20);

az.plot_autocorr(traceS, var_names="xi_i", max_lag=20, combined=True, figsize=(29,25));
plt.suptitle("Separate ACF plots for the GEV shape param", fontsize=20);

az.plot_autocorr(traceS, var_names="mu_i", max_lag=20, combined=True, figsize=(20,20));
plt.suptitle("Separate ACF plots for the GEV location param", fontsize=20);

az.plot_autocorr(traceS, var_names="sigma_i", max_lag=20, combined=True, figsize=(29,25));
plt.suptitle("Separate ACF plots for the GEV scale param", fontsize=20);

az.plot_mcse(traceC, var_names=["xi_i"], figsize=(20,20), drawstyle="steps-mid", color="k",
    linestyle="-", marker=None, rug=True,  rug_kwargs={"color": "r"});
plt.suptitle("Common Monte Carlo SE plots for the GEV shape param", fontsize=20);

az.plot_mcse(traceC, var_names=["mu_i"], figsize=(29,25), drawstyle="steps-mid", color="k",
    linestyle="-", marker=None, rug=True,  rug_kwargs={"color": "r"});
plt.suptitle("Common Monte Carlo SE plots for the GEV location param", fontsize=20);

az.plot_mcse(traceC, var_names=["sigma_i"], figsize=(20,20), drawstyle="steps-mid", color="k",
    linestyle="-", marker=None, rug=True,  rug_kwargs={"color": "r"});
plt.suptitle("Common Monte Carlo SE plots for the GEV scale param", fontsize=20);

az.plot_mcse(traceH, var_names=["xi_i"], figsize=(29,25), drawstyle="steps-mid", color="k",
    linestyle="-", marker=None, rug=True,  rug_kwargs={"color": "r"});
plt.suptitle("Hilbert Monte Carlo SE plots for the GEV shape param", fontsize=20);

az.plot_mcse(traceH, var_names=["mu_i"], figsize=(20,20), drawstyle="steps-mid", color="k",
    linestyle="-", marker=None, rug=True,  rug_kwargs={"color": "r"});
plt.suptitle("Hilbert Monte Carlo SE plots for the GEV location param", fontsize=20);

az.plot_mcse(traceH, var_names=["sigma_i"], figsize=(25,25), drawstyle="steps-mid", color="k",
    linestyle="-", marker=None, rug=True,  rug_kwargs={"color": "r"});
plt.suptitle("Hilbert Monte Carlo SE plots for the GEV scale param", fontsize=20);

az.plot_mcse(traceL, var_names=["xi_i"], figsize=(29,25), drawstyle="steps-mid", color="k",
    linestyle="-", marker=None, rug=True,  rug_kwargs={"color": "r"});
plt.suptitle("Latent Monte Carlo SE plots for the GEV shape param", fontsize=20);

az.plot_mcse(traceL, var_names=["mu_i"], figsize=(20,20), drawstyle="steps-mid", color="k",
    linestyle="-", marker=None, rug=True,  rug_kwargs={"color": "r"});
plt.suptitle("Latent Monte Carlo SE plots for the GEV location param", fontsize=20);

az.plot_mcse(traceL, var_names=["sigma_i"], figsize=(29,25), drawstyle="steps-mid", color="k",
    linestyle="-", marker=None, rug=True,  rug_kwargs={"color": "r"});
plt.suptitle("Latent Monte Carlo SE plots for the GEV scale param", fontsize=20);

az.plot_mcse(traceS, var_names=["xi_i"], figsize=(20,20), drawstyle="steps-mid", color="k",
    linestyle="-", marker=None, rug=True,  rug_kwargs={"color": "r"});
plt.suptitle("Separate Monte Carlo SE plots for the GEV shape param", fontsize=20);

az.plot_mcse(traceS, var_names=["mu_i"], figsize=(29,25), drawstyle="steps-mid", color="k",
    linestyle="-", marker=None, rug=True,  rug_kwargs={"color": "r"});
plt.suptitle("Separate Monte Carlo SE plots for the GEV location param", fontsize=20);

az.plot_mcse(traceS, var_names=["sigma_i"], figsize=(20,20), drawstyle="steps-mid", color="k",
    linestyle="-", marker=None, rug=True,  rug_kwargs={"color": "r"});
plt.suptitle("Separate Monte Carlo SE plots for the GEV scale param", fontsize=20);

comparison_dict = {"Common":traceC, "Hilbert":traceH,  "Latent":traceL, "Separate":traceS}
comp_int = {"Hilbert": traceH, "Latent":traceL}

model_comparison = az.compare(comparison_dict,ic="loo", method="BB-pseudo-BMA")
model_comparison

/home/sm_seuba/.conda/envs/pymc_env/lib/python3.11/site-packages/arviz/stats/stats.py:795: UserWarning:

Estimated shape parameter of Pareto distribution is greater than 0.70 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.

/home/sm_seuba/.conda/envs/pymc_env/lib/python3.11/site-packages/arviz/stats/stats.py:795: UserWarning:

Estimated shape parameter of Pareto distribution is greater than 0.70 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.

model_comparison = az.compare(comparison_dict,ic="loo", method="stacking")
model_comparison

/home/sm_seuba/.conda/envs/pymc_env/lib/python3.11/site-packages/arviz/stats/stats.py:795: UserWarning:

Estimated shape parameter of Pareto distribution is greater than 0.70 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.

/home/sm_seuba/.conda/envs/pymc_env/lib/python3.11/site-packages/arviz/stats/stats.py:795: UserWarning:

Estimated shape parameter of Pareto distribution is greater than 0.70 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.

model_int_comparison = az.compare(comp_int,ic="loo") # method="BB-pseudo-BMA")
model_int_comparison

model_int_comparison = az.compare(comp_int,ic="loo", method="stacking")
model_int_comparison

az.plot_compare(model_comparison, insample_dev=True);

idataC.extend(traceC)
with model_C:
    pm.sample_posterior_predictive(idataC, extend_inferencedata=True, random_seed=rng)

Sampling: [Ys_observed]

idataH.extend(traceH)
with model_H:
    pm.sample_posterior_predictive(idataH, extend_inferencedata=True, random_seed=rng)

Sampling: [Ys_observed]

idataL.extend(traceL)
with model_L:
    pm.sample_posterior_predictive(idataL, extend_inferencedata=True, random_seed=rng)

Sampling: [Ys_observed]

idataS.extend(traceS)
with model_S:
    pm.sample_posterior_predictive(idataS, extend_inferencedata=True, random_seed=rng)

Sampling: [Ys_observed]

all_idata = [idataC, idataH, idataL, idataS]
models_names = ["Common", "Hilbert", "Latent", "Separate"]
all_idata_models_dict = {models_names[i]:all_idata[i] for i in range(len(models_names))}

fig, axes = plt.subplots(nrows=4, ncols=2, figsize=(20,20))
row_titles = ["PPC plot (KDE)", "PPC plot (ECDF)"]
for row, idata in enumerate(all_idata):
    az.plot_ppc(idata, observed=True, observed_rug=True, figsize=(12, 6), ax=axes[row,0]);
    az.plot_ppc(idata, observed=True, observed_rug=True,  kind='cumulative', figsize=(12, 6), ax=axes[row,1]);
    axes[row,0].set_ylabel(models_names[row], fontsize=20);
for col, row_title in enumerate(row_titles):
    axes[0,col].title.set_text(row_title);
fig.suptitle("Visualization of Posterior Predictive Checks 1:2", fontsize=30);
fig.tight_layout();
fig.show();

fig, axes = plt.subplots(nrows=4, ncols=3, figsize=(20,20))
row_titles = ["LOO-PIT", "LOO-PIT-ECDF","Bayesian p-value"]
for row, idata in enumerate(all_idata):
    az.plot_loo_pit(idata=idata, y="Ys_observed", ax=axes[row,0]);
    az.plot_loo_pit(idata=idata, y="Ys_observed",y_hat="Ys_observed", ecdf=True, hdi_prob=.95,  ax=axes[row,1]);
    az.plot_bpv(idata, kind="p_value", ax=axes[row,2]);
    axes[row,0].set_ylabel(models_names[row], fontsize=20)
for col, row_title in enumerate(row_titles):
    axes[0,col].title.set_text(row_title);
fig.suptitle("Visualization of Posterior Predictive Checks 2:2", fontsize=30);
fig.tight_layout();
fig.show();

az.rhat(idataC, var_names=gev_var_names, method="rank")

<xarray.Dataset>
Dimensions:  (station: 36)
Coordinates:
  * station  (station) int64 0 1 2 3 4 5 6 7 8 9 ... 27 28 29 30 31 32 33 34 35
Data variables:
    xi_i     (station) float64 1.0 1.0 1.0 1.0 1.0 1.0 ... 1.0 1.0 1.0 1.0 1.0
    mu_i     (station) float64 1.0 1.0 1.0 1.0 1.0 ... 1.0 1.0 0.9999 1.0 1.0
    sigma_i  (station) float64 1.0 1.0 1.0 1.0 1.0 1.0 ... 1.0 1.0 1.0 1.0 1.0

array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35])

array([1.00012348, 1.000118  , 0.99999912, 1.00003813, 1.00006023,
       1.00044718, 1.00007157, 1.00011735, 1.00025   , 1.00001348,
       1.00012625, 1.00014164, 1.00001506, 1.00010566, 1.00034514,
       1.00010828, 0.99999159, 1.00024789, 1.00017961, 1.00028053,
       1.00024899, 1.00002209, 1.00002071, 1.00016833, 1.0001346 ,
       1.00010401, 1.00017937, 1.00016297, 1.00002132, 1.00001207,
       1.00021816, 1.00020154, 1.00010145, 1.00007503, 1.00008113,
       1.00019083])

array([1.00010822, 1.0001588 , 1.00003128, 1.00008048, 1.00008232,
       1.00012643, 0.99999267, 1.00002591, 1.00000141, 1.00002357,
       1.00029664, 1.00004124, 1.00015723, 1.00006441, 0.99998331,
       1.00009504, 1.00000443, 0.99997786, 1.00008963, 1.00005259,
       1.00022517, 0.99999681, 1.00003091, 0.99999009, 1.00003945,
       1.00028005, 1.00025709, 1.00005993, 1.0000313 , 1.00023374,
       1.00013701, 0.99997854, 1.00010976, 0.99994792, 1.00001314,
       1.00008804])

array([0.99998881, 1.00008464, 1.00001915, 1.00002624, 1.00007878,
       1.00013609, 1.00009968, 1.00008019, 0.9999941 , 1.00022457,
       1.00006224, 1.00001573, 1.00001477, 1.0000793 , 1.00014071,
       1.00011242, 1.0003145 , 1.0002158 , 1.00023461, 1.00015322,
       1.00023628, 1.00009116, 1.00010531, 1.00007147, 0.99997905,
       1.00005232, 0.99998972, 1.00007853, 1.00009501, 0.99998551,
       1.00012318, 1.00008699, 1.00015634, 1.00007187, 1.00003861,
       1.0000979 ])

PandasIndex(Index([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35],
      dtype='int64', name='station'))

az.rhat(idataH, var_names=gev_var_names, method="rank")

<xarray.Dataset>
Dimensions:  (station: 36)
Coordinates:
  * station  (station) int64 0 1 2 3 4 5 6 7 8 9 ... 27 28 29 30 31 32 33 34 35
Data variables:
    xi_i     (station) float64 1.0 1.0 1.0 1.0 1.0 1.0 ... 1.0 1.0 1.0 1.0 1.0
    mu_i     (station) float64 1.0 1.0 1.0 1.0 1.0 1.0 ... 1.0 1.0 1.0 1.0 1.0
    sigma_i  (station) float64 1.0 1.0 1.0 1.0 1.0 1.0 ... 1.0 1.0 1.0 1.0 1.0

array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35])

array([0.99996565, 0.99995291, 0.99998121, 1.00007448, 1.00016995,
       1.00008598, 0.99997321, 0.99997215, 1.00004475, 1.00000786,
       0.9999975 , 1.00009733, 1.0001826 , 1.00003557, 0.99999101,
       0.99995952, 0.99995823, 1.00001115, 1.00006038, 0.99999417,
       0.99994726, 0.99997036, 1.00009194, 1.00007038, 1.00003876,
       1.00003622, 1.00006607, 1.00005562, 0.9999562 , 0.99997654,
       0.99994813, 1.00002265, 0.99995277, 1.00002878, 1.00006476,
       1.00009616])

array([0.99999756, 0.99999726, 1.00005905, 1.00003629, 0.99999375,
       0.99999225, 1.00003083, 1.00003932, 1.00010581, 0.99999427,
       1.00000245, 1.00006661, 1.00002612, 1.00013666, 1.00003612,
       0.99997523, 1.00022296, 1.00006578, 1.00002643, 0.99997986,
       1.00009579, 1.00011116, 1.00008181, 1.00001002, 1.00016267,
       0.99997235, 0.99999983, 0.99995795, 0.99997604, 1.00007861,
       1.00006592, 1.00003477, 1.00018802, 0.99996606, 0.99998001,
       1.00008464])

array([1.00004106, 1.00020855, 1.00012484, 1.00009909, 1.00001494,
       1.00013115, 1.00025154, 1.0000722 , 1.00000125, 0.99996845,
       1.00011668, 1.00007451, 1.00005348, 1.00000404, 1.00001398,
       1.00003009, 0.99999633, 1.00001202, 1.00007883, 1.00001686,
       1.00001517, 1.00004338, 1.00001609, 1.0000979 , 0.9999972 ,
       0.99996282, 1.0001573 , 1.00019808, 0.99999833, 1.00011803,
       1.00002002, 1.00001337, 1.000108  , 1.00000628, 1.00007357,
       1.00005598])

PandasIndex(Index([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35],
      dtype='int64', name='station'))

az.rhat(idataL, var_names=gev_var_names, method="rank")

<xarray.Dataset>
Dimensions:  (station: 36)
Coordinates:
  * station  (station) int64 0 1 2 3 4 5 6 7 8 9 ... 27 28 29 30 31 32 33 34 35
Data variables:
    xi_i     (station) float64 1.0 1.0 1.0 1.0 1.0 1.0 ... 1.0 1.0 1.0 1.0 1.0
    mu_i     (station) float64 1.0 1.0 1.0 1.0 1.0 1.0 ... 1.0 1.0 1.0 1.0 1.0
    sigma_i  (station) float64 1.0 1.0 1.0 1.0 1.0 ... 1.0 1.0 1.0 0.9999 1.0

array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35])

array([1.00001646, 1.00008482, 0.99998408, 1.00003951, 1.00011252,
       1.00005202, 1.00002725, 1.00003746, 1.00001264, 0.99999372,
       1.00002093, 1.00001924, 1.00005128, 0.99999322, 1.00009668,
       1.00019316, 1.00002442, 0.99999419, 1.00000237, 1.00002759,
       1.00003363, 1.00001633, 1.00000158, 1.00004905, 0.99999046,
       1.00000321, 0.99994909, 1.00009936, 1.00002461, 1.00006563,
       1.00015065, 1.000019  , 0.99998436, 0.99996221, 1.00000866,
       0.99999758])

array([1.00003618, 1.00001879, 1.00011999, 1.00002691, 0.99997701,
       0.99995492, 1.00007443, 1.00000592, 1.00003018, 1.00009648,
       1.00000688, 1.00006867, 1.00018335, 1.0000421 , 1.00011954,
       1.00000712, 1.00013084, 0.99997595, 1.00003114, 1.00004817,
       1.00010185, 1.00002662, 0.99995703, 1.00009873, 0.99996347,
       1.00005468, 0.99996686, 1.0000649 , 0.99997368, 1.00007717,
       0.99998574, 0.99998343, 1.00014441, 1.00011272, 0.99996952,
       1.00021822])

array([1.00002306, 1.0000127 , 1.00015749, 1.00004753, 1.00013266,
       1.00003044, 1.00002231, 1.00006432, 1.00007808, 1.00005451,
       1.00014174, 1.00004655, 0.99999253, 1.00019849, 1.0000274 ,
       1.00004269, 1.00035742, 1.00011933, 1.00003926, 1.00002438,
       1.0000178 , 1.0003335 , 1.0000234 , 1.00009796, 1.00017614,
       0.99992886, 0.99999566, 0.99998598, 1.00007676, 1.00008915,
       1.0000783 , 1.00001965, 1.00024657, 1.00015741, 0.99994536,
       1.00007012])

PandasIndex(Index([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35],
      dtype='int64', name='station'))

az.rhat(idataS, var_names=gev_var_names, method="rank")

<xarray.Dataset>
Dimensions:  (station: 36)
Coordinates:
  * station  (station) int64 0 1 2 3 4 5 6 7 8 9 ... 27 28 29 30 31 32 33 34 35
Data variables:
    xi_i     (station) float64 1.0 1.0 1.0 1.0 1.0 1.0 ... 1.0 1.0 1.0 1.0 1.0
    mu_i     (station) float64 1.0 1.0 1.0 1.0 1.0 1.0 ... 1.0 1.0 1.0 1.0 1.0
    sigma_i  (station) float64 1.0 1.0 1.0 1.0 1.0 1.0 ... 1.0 1.0 1.0 1.0 1.0

array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35])

array([1.00009834, 1.00016629, 0.99999356, 1.0000242 , 1.00001685,
       1.00005038, 1.00001791, 1.00021836, 1.00003317, 1.00004126,
       1.00001791, 1.00002098, 1.00014191, 1.00010258, 0.99997807,
       1.00003284, 1.00004011, 1.00001637, 1.00010735, 1.00010482,
       0.99997848, 1.00008759, 1.00011422, 1.00023423, 1.0000079 ,
       1.00008115, 1.00008332, 0.99997139, 1.00005982, 1.00027301,
       1.00015921, 1.00006201, 1.00012065, 1.00010622, 1.0000929 ,
       1.00007199])

array([1.00013463, 1.00025058, 1.00004315, 1.00011194, 1.00001635,
       0.99998484, 1.00004586, 1.00016766, 0.99996769, 1.00002586,
       1.00010328, 1.00002735, 1.00004061, 1.00002047, 1.00007661,
       1.00024793, 1.00014618, 1.00015954, 1.00000595, 1.00006391,
       1.00011509, 1.00001143, 1.00013153, 1.00004378, 0.99997426,
       1.0001581 , 0.99999221, 1.00005708, 1.00016907, 1.00006322,
       0.99999407, 1.00008705, 1.00021583, 1.00001639, 1.00003834,
       1.00007977])

array([1.00004143, 1.00018836, 0.9999781 , 1.00004998, 1.00003313,
       1.00009757, 0.99999888, 0.99997236, 1.00004672, 1.00007534,
       1.00013137, 1.00008288, 0.99997581, 1.0000322 , 0.99997683,
       1.0000025 , 1.00007862, 1.00075558, 1.00020537, 0.99995969,
       1.00007359, 0.9999594 , 1.0000753 , 1.00025827, 1.00002026,
       1.00006716, 1.00023756, 1.00002792, 1.00005758, 1.00028237,
       0.99998545, 1.00013602, 1.00032281, 1.00001382, 0.99998904,
       1.00011238])

PandasIndex(Index([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35],
      dtype='int64', name='station'))

az.rhat(idataC, var_names=gev_var_names, method="folded")

<xarray.Dataset>
Dimensions:  (station: 36)
Coordinates:
  * station  (station) int64 0 1 2 3 4 5 6 7 8 9 ... 27 28 29 30 31 32 33 34 35
Data variables:
    xi_i     (station) float64 1.0 1.0 1.0 1.0 1.0 1.0 ... 1.0 1.0 1.0 1.0 1.0
    mu_i     (station) float64 1.0 1.0 1.0 1.0 1.0 ... 1.0 1.0 0.9999 1.0 1.0
    sigma_i  (station) float64 1.0 1.0 1.0 1.0 1.0 1.0 ... 1.0 1.0 1.0 1.0 1.0

array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35])

array([1.00012348, 1.00009489, 0.99995983, 0.99998236, 1.00006023,
       1.00007256, 1.00007157, 1.00005313, 1.00025   , 1.00001348,
       1.00012625, 1.00010986, 0.99999801, 1.00010566, 1.00034514,
       1.00010828, 0.99996784, 1.00024789, 1.00012751, 1.00013102,
       1.00024899, 1.00002209, 1.00002071, 0.99998303, 1.00006359,
       1.00010401, 1.00017937, 1.00000111, 1.00002132, 0.9999958 ,
       1.00021816, 1.00020154, 0.99998562, 1.00007503, 1.00000438,
       1.00019083])

array([1.00010822, 1.0001588 , 1.00003128, 1.00008048, 1.00008232,
       1.00012643, 0.99999267, 1.00002591, 0.99996073, 0.99999053,
       1.00029664, 1.00004124, 1.00015723, 1.00006441, 0.99998331,
       1.00009504, 1.00000443, 0.99997786, 0.99996289, 1.00005259,
       1.0000722 , 0.99999681, 1.00003091, 0.99995027, 1.00003945,
       1.00028005, 1.00025709, 0.99993896, 1.0000313 , 1.00023374,
       1.00013701, 0.99996608, 1.00010976, 0.99991642, 1.00001165,
       1.00008804])

array([0.99998881, 1.00008464, 0.99995513, 1.00002624, 1.00007878,
       1.00011781, 1.00009968, 1.00008019, 0.9999941 , 1.00022457,
       1.00006224, 1.00001573, 1.00001477, 1.0000793 , 1.00014071,
       1.00011242, 0.99993834, 1.0002158 , 1.00023461, 0.99995107,
       1.00014876, 1.00009116, 0.9999961 , 1.00007147, 0.99997905,
       1.00005232, 0.99998972, 1.00007853, 1.00009501, 0.99998551,
       1.00012318, 1.00008699, 0.99999698, 1.00007187, 1.00003861,
       1.0000979 ])

PandasIndex(Index([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35],
      dtype='int64', name='station'))

az.rhat(idataH, var_names=gev_var_names, method="folded")

<xarray.Dataset>
Dimensions:  (station: 36)
Coordinates:
  * station  (station) int64 0 1 2 3 4 5 6 7 8 9 ... 27 28 29 30 31 32 33 34 35
Data variables:
    xi_i     (station) float64 1.0 1.0 0.9999 1.0 1.0 ... 1.0 1.0 1.0 1.0 1.0
    mu_i     (station) float64 1.0 1.0 1.0 1.0 1.0 ... 1.0 1.0 1.0 0.9999 1.0
    sigma_i  (station) float64 0.9999 1.0 1.0 1.0 1.0 ... 1.0 1.0 0.9999 1.0 1.0

array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35])

array([0.99996565, 0.99995291, 0.99993775, 1.00007448, 1.00016995,
       1.00008598, 0.99997321, 0.99997215, 1.00004475, 1.00000786,
       0.9999975 , 1.00009733, 1.0001826 , 1.00003557, 0.99999101,
       0.99995152, 0.99994963, 1.00001115, 1.00006038, 0.99999417,
       0.99994399, 0.99997036, 1.00009194, 1.00007038, 1.00003876,
       1.00003622, 1.00006607, 1.00005562, 0.9999562 , 0.99997654,
       0.99994813, 1.00002265, 0.99995277, 0.99995231, 1.00006476,
       1.00009616])

array([0.99999756, 0.99999726, 0.99997395, 1.00003629, 0.99999375,
       0.99999225, 1.00003083, 1.00003932, 1.00010581, 0.99999427,
       0.99995886, 1.00006661, 1.00002612, 1.00013666, 1.00003612,
       0.99997523, 1.00022296, 1.00006578, 1.00000378, 0.99994391,
       0.99995965, 1.00011116, 1.00008181, 0.99999804, 1.00016267,
       0.99996662, 0.99997196, 0.99995795, 0.99996569, 1.00007861,
       1.00006592, 1.00003477, 1.00018802, 0.99996606, 0.9999423 ,
       1.00005289])

array([0.99993303, 1.00003808, 1.00012484, 1.00009909, 0.99999201,
       1.00000274, 0.99992416, 1.0000722 , 0.99992418, 0.99994208,
       1.00003404, 1.00007451, 1.00002523, 1.00000404, 0.99997714,
       1.00003009, 0.99998806, 0.99993997, 1.00007883, 1.00000277,
       1.00000894, 0.99997516, 1.00001609, 1.00002109, 0.9999972 ,
       0.99996282, 1.0001573 , 0.99995661, 0.99998013, 1.00011803,
       0.99996151, 1.00001337, 1.00004414, 0.99994986, 1.00001472,
       1.00005598])

PandasIndex(Index([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35],
      dtype='int64', name='station'))

az.rhat(idataL, var_names=gev_var_names, method="folded")

<xarray.Dataset>
Dimensions:  (station: 36)
Coordinates:
  * station  (station) int64 0 1 2 3 4 5 6 7 8 9 ... 27 28 29 30 31 32 33 34 35
Data variables:
    xi_i     (station) float64 1.0 1.0 1.0 1.0 1.0 ... 0.9999 0.9999 0.9999 1.0
    mu_i     (station) float64 1.0 1.0 1.0 1.0 1.0 1.0 ... 1.0 1.0 1.0 1.0 1.0
    sigma_i  (station) float64 1.0 1.0 1.0 1.0 1.0 ... 1.0 1.0 1.0 0.9999 1.0

array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35])

array([1.00001312, 1.00007896, 0.99998408, 1.00003951, 0.99997932,
       1.00005202, 1.00002725, 1.00003746, 1.00001264, 0.99998515,
       1.00002093, 1.00001924, 1.00005128, 0.99994387, 1.00009668,
       1.00019316, 0.9999855 , 0.99995742, 0.99991767, 1.00000983,
       0.99990704, 0.99995375, 0.99996944, 1.00002588, 0.99994119,
       0.99998502, 0.99994909, 1.00009936, 0.99998024, 1.00006563,
       0.99997609, 1.000019  , 0.99993436, 0.99992191, 0.99993706,
       0.99997986])

array([1.00003249, 1.00001879, 1.00011999, 0.99998947, 0.99997701,
       0.99995492, 1.00007443, 1.00000592, 0.9999586 , 1.00009648,
       1.00000688, 0.99994145, 1.00018335, 1.0000421 , 1.00011954,
       1.00000712, 1.00013084, 0.99997595, 1.00003114, 1.00004817,
       1.00010185, 1.00002082, 0.99993283, 1.00009873, 0.99996347,
       1.00005468, 0.99993888, 1.0000649 , 0.9999692 , 0.99993007,
       0.9999491 , 0.99998343, 1.00014441, 1.00011272, 0.99996952,
       1.00021822])

array([0.99998585, 1.0000127 , 1.00004617, 0.99997009, 1.00013266,
       0.99997977, 0.99996265, 1.00006432, 1.00007808, 0.99992073,
       1.00014174, 0.99999771, 0.99999253, 1.00019849, 1.0000274 ,
       1.00001812, 1.00005775, 1.00011933, 1.00003866, 1.00001879,
       0.9999998 , 1.00003271, 0.9999784 , 1.00009796, 1.00017614,
       0.99992261, 0.99999566, 0.99993119, 1.00007676, 0.9999677 ,
       1.0000783 , 0.99998316, 1.00012807, 0.99998363, 0.99994536,
       1.0000305 ])

PandasIndex(Index([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35],
      dtype='int64', name='station'))

az.rhat(idataS, var_names=gev_var_names, method="folded")

<xarray.Dataset>
Dimensions:  (station: 36)
Coordinates:
  * station  (station) int64 0 1 2 3 4 5 6 7 8 9 ... 27 28 29 30 31 32 33 34 35
Data variables:
    xi_i     (station) float64 1.0 1.0 1.0 1.0 1.0 1.0 ... 1.0 1.0 1.0 1.0 1.0
    mu_i     (station) float64 1.0 1.0 1.0 1.0 1.0 1.0 ... 1.0 1.0 1.0 1.0 1.0
    sigma_i  (station) float64 1.0 1.0 1.0 1.0 1.0 1.0 ... 1.0 1.0 1.0 1.0 1.0

array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35])

array([1.00009834, 1.00016629, 0.99999356, 1.00001827, 1.00001685,
       1.00001422, 1.00000337, 1.00021787, 0.99995334, 1.00004126,
       1.00001713, 0.9999667 , 1.00014191, 1.00008954, 0.99997807,
       1.00003284, 1.00004011, 1.00000352, 1.00001251, 1.00010291,
       0.99997848, 1.00008759, 1.00011422, 1.00023423, 0.99998156,
       0.99995757, 1.00002147, 0.99997139, 1.00005982, 1.0001156 ,
       1.00015921, 1.00006201, 1.00012065, 1.00010622, 1.0000929 ,
       1.00007199])

array([1.00013463, 0.99999586, 0.99999966, 1.00011194, 0.9999952 ,
       0.99998484, 0.99998177, 0.99994874, 0.99996769, 1.00002586,
       1.00010328, 1.00002735, 1.00001019, 0.99999145, 0.99996814,
       1.00024793, 1.00001278, 1.00015954, 1.00000595, 1.00003017,
       1.00011509, 1.00001143, 1.00013153, 1.00004378, 0.99997426,
       1.0001581 , 0.99999221, 0.99994316, 1.00016907, 1.00006322,
       0.99999407, 0.99995744, 1.00021583, 1.00001639, 1.00003834,
       1.00007977])

array([1.00003503, 1.00018836, 0.99997459, 1.00004998, 1.00003313,
       1.00009757, 0.99999848, 0.99991643, 0.9999557 , 1.00000272,
       1.00007972, 1.00008288, 0.99997581, 1.0000322 , 0.99997683,
       1.0000025 , 1.00007862, 1.00075558, 1.00008925, 0.99994296,
       1.00007359, 0.9999594 , 1.00006006, 1.00010102, 0.99997778,
       1.00006716, 1.00023756, 0.99999473, 1.00005758, 1.00000074,
       0.99998545, 0.99997385, 1.00032281, 1.00001382, 0.99998904,
       1.00011238])

PandasIndex(Index([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35],
      dtype='int64', name='station'))

az.ess(idataC, var_names=gev_var_names, method="local",prob=[0.05,0.95])

<xarray.Dataset>
Dimensions:  (station: 36)
Coordinates:
  * station  (station) int64 0 1 2 3 4 5 6 7 8 9 ... 27 28 29 30 31 32 33 34 35
Data variables:
    xi_i     (station) float64 2.264e+04 2.142e+04 ... 2.305e+04 2.104e+04
    mu_i     (station) float64 2.393e+04 2.912e+04 ... 2.917e+04 2.215e+04
    sigma_i  (station) float64 2.21e+04 2.584e+04 ... 2.026e+04 2.52e+04

array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35])

array([22641.157996  , 21416.48053858, 21180.71486932, 20230.72077683,
       22210.87743817, 22165.60447367, 20575.85677374, 21336.37353302,
       21938.22418251, 22186.76207976, 20805.49524227, 20139.33747108,
       22343.38564443, 19654.95686382, 22754.63697214, 20947.65362865,
       20872.99203282, 21566.22480574, 19960.20269526, 21545.36613071,
       22264.00488413, 21778.33180117, 23646.92191463, 20942.38173764,
       19620.04661486, 22512.84102225, 21669.3349779 , 21871.54811114,
       20170.54828959, 21038.68387809, 20385.58982785, 19556.36438671,
       19989.30744701, 20196.17442734, 23050.75363066, 21037.55553918])

array([23934.31966223, 29124.8951694 , 30102.8716542 , 27149.04824335,
       29456.57803258, 28823.04569316, 28510.40956184, 31848.20361761,
       29787.80346786, 30051.21334922, 27159.15182642, 28963.45262292,
       29660.44172903, 28447.80962796, 30166.60015492, 23711.32505274,
       30683.56146392, 23017.80340094, 28507.24462347, 29945.58149498,
       29283.51067649, 29543.97274993, 28418.75343467, 29516.6011866 ,
       29369.97301698, 22407.34746591, 28641.13352549, 30358.17452491,
       24167.64515677, 29409.52141486, 29796.98914399, 25777.44521356,
       29243.37760336, 27041.90714704, 29173.67158555, 22153.3675854 ])

array([22101.72196422, 25842.03896008, 19265.13153399, 27285.65190626,
       21815.91197592, 26940.26944838, 19858.35588811, 25347.22529319,
       20028.22015709, 19380.12501519, 29852.60734017, 27257.3829933 ,
       25525.96711876, 22417.61508904, 25463.94894581, 23833.83935582,
       28518.77534848, 23289.39746176, 26127.08833051, 26654.96518309,
       16258.24822571, 27423.23950413, 26313.28937032, 24123.77851281,
       29976.00698359, 23242.5987712 , 27779.75094518, 25643.3973813 ,
       23424.85077045, 27300.85078198, 20339.79338397, 26683.94565946,
       25780.55233137, 27680.92702935, 20260.57980832, 25204.11613993])

PandasIndex(Index([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35],
      dtype='int64', name='station'))

az.ess(idataH, var_names=gev_var_names, method="local",prob=[0.05,0.95])

<xarray.Dataset>
Dimensions:  (station: 36)
Coordinates:
  * station  (station) int64 0 1 2 3 4 5 6 7 8 9 ... 27 28 29 30 31 32 33 34 35
Data variables:
    xi_i     (station) float64 3.415e+04 3.464e+04 ... 3.441e+04 3.247e+04
    mu_i     (station) float64 3.592e+04 3.434e+04 ... 3.408e+04 3.423e+04
    sigma_i  (station) float64 3.878e+04 3.585e+04 ... 3.621e+04 3.562e+04

array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35])

array([34148.69967526, 34638.16113392, 30807.56836704, 32314.78408559,
       31609.97647253, 35124.12986736, 30572.01700728, 34929.1064713 ,
       31801.98754987, 34552.22390754, 32382.21535921, 33238.20836785,
       34919.91983671, 33438.7471558 , 35275.04082406, 36142.03411342,
       34073.05710117, 35448.54476834, 32652.8845928 , 34303.66927069,
       32498.9191409 , 35823.83314199, 31687.06741089, 32785.00069424,
       34055.18605685, 32913.08108033, 30158.82942785, 35157.00245204,
       31123.03927704, 32004.48878225, 32400.15061088, 35382.24935998,
       30066.06521603, 29704.77385699, 34406.3740515 , 32466.77247856])

array([35915.39543947, 34338.96413865, 32413.90094417, 34686.73388654,
       32035.49985591, 32655.53596208, 33315.85026177, 33444.59960671,
       33000.9121897 , 34685.44954236, 34298.87773739, 33815.18480721,
       33552.94886662, 34980.79512642, 33654.69903396, 37212.77482638,
       34169.41820632, 36643.6042728 , 35336.25788315, 33396.20982807,
       35005.34082362, 34550.02676361, 33255.04010307, 33322.59542088,
       34224.63666746, 37086.63373096, 33635.53130572, 33690.25289893,
       36827.85180133, 32901.26060981, 31981.46182694, 36166.1358692 ,
       32953.3341611 , 35773.76735049, 34084.16843306, 34234.34010587])

array([38782.07249481, 35846.80667828, 33691.15454974, 37521.69551581,
       33414.54888382, 35378.6497867 , 34573.5852952 , 35272.06438496,
       34535.02258821, 37656.40768635, 37487.17198862, 38740.5269059 ,
       36702.88570329, 36744.12107959, 35355.51287658, 36942.18401024,
       35737.80865303, 38712.46183112, 36535.58011383, 37985.68677858,
       35508.65371874, 34431.07611097, 37402.55239807, 38517.96004526,
       37090.60912346, 37141.33575815, 33131.92102448, 36416.37046323,
       38042.88184748, 37170.11512601, 37024.22881819, 38254.53133275,
       34837.47074749, 36151.66046507, 36205.15881115, 35618.28069416])

PandasIndex(Index([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35],
      dtype='int64', name='station'))

az.ess(idataL, var_names=gev_var_names, method="local",prob=[0.05,0.95])

<xarray.Dataset>
Dimensions:  (station: 36)
Coordinates:
  * station  (station) int64 0 1 2 3 4 5 6 7 8 9 ... 27 28 29 30 31 32 33 34 35
Data variables:
    xi_i     (station) float64 3.687e+04 3.454e+04 ... 3.433e+04 3.383e+04
    mu_i     (station) float64 3.314e+04 3.015e+04 ... 3.118e+04 3.129e+04
    sigma_i  (station) float64 3.505e+04 3.474e+04 ... 3.724e+04 3.403e+04

array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35])

array([36867.67838453, 34542.8972876 , 32414.6986446 , 35077.52212479,
       31301.21839629, 34770.34102101, 34261.72501781, 35717.06695418,
       34226.29801615, 37985.80854893, 34507.02547788, 34893.49383392,
       34788.94353323, 36391.91926539, 35182.4617199 , 35923.45345817,
       35336.74896675, 34685.43596484, 35955.38545761, 34685.40648918,
       34316.24154987, 35910.26500683, 35914.0506371 , 35341.68201564,
       35965.53379032, 32842.94140154, 31472.94665665, 35073.50296186,
       34922.11048451, 31922.98663132, 35466.78101277, 35858.18984575,
       31100.62268928, 32586.88463874, 34333.0208774 , 33827.15465196])

array([33137.74056503, 30145.14680874, 30978.76798731, 30811.67364614,
       27278.96403038, 29348.28731842, 30538.110594  , 30289.04355905,
       29910.16048753, 33602.77247075, 31901.40208537, 30593.75667728,
       29123.07214952, 32286.94708306, 32430.69153466, 35258.78866867,
       31425.24601224, 32014.77526894, 33945.83703367, 32283.27961048,
       29653.76068022, 33702.42663344, 30860.80379839, 32520.95473806,
       31299.57020934, 33447.7772953 , 32409.8442517 , 30408.28385167,
       35631.1463479 , 30571.55319721, 31176.5081747 , 33421.88078321,
       32537.95869761, 34325.69823792, 31179.3344071 , 31286.28272806])

array([35046.70139394, 34738.11164429, 35041.50387877, 37584.63205304,
       32295.75087363, 38242.88126501, 35135.85026409, 33939.61755702,
       34547.08287348, 38600.39884224, 34858.22416373, 37315.2491152 ,
       37668.4570319 , 38263.25214444, 33849.39513044, 36150.42218138,
       36345.80060146, 36043.37250296, 37086.73896282, 35487.31814766,
       35915.02984541, 36025.58854062, 36694.26255509, 36048.82465231,
       37450.79853666, 35982.18403303, 34933.08983432, 35283.50318812,
       34156.15090882, 37439.81015818, 34659.21292056, 36805.70254913,
       35679.53503751, 34873.35392196, 37236.62390024, 34034.89361956])

PandasIndex(Index([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35],
      dtype='int64', name='station'))

az.ess(idataS, var_names=gev_var_names, method="local",prob=[0.05,0.95])

<xarray.Dataset>
Dimensions:  (station: 36)
Coordinates:
  * station  (station) int64 0 1 2 3 4 5 6 7 8 9 ... 27 28 29 30 31 32 33 34 35
Data variables:
    xi_i     (station) float64 2.38e+04 2.22e+04 ... 2.491e+04 2.066e+04
    mu_i     (station) float64 2.37e+04 3.028e+04 ... 3.05e+04 2.337e+04
    sigma_i  (station) float64 2.061e+04 2.211e+04 ... 2.113e+04 2.324e+04

array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35])

array([23804.0707776 , 22204.8204408 , 23305.17597248, 21277.5665275 ,
       20976.67566811, 23996.63696841, 22583.70301444, 20064.10870684,
       24924.7595785 , 25329.79047552, 21846.59140534, 22760.22867932,
       22359.91176523, 20705.15508586, 21742.25948543, 21164.89032048,
       23617.67339801, 21981.09906813, 23873.21360266, 22979.52495888,
       22695.53159703, 20890.7631399 , 22322.35075719, 23897.11194813,
       23214.53373463, 21913.68899835, 23608.92854497, 21697.69317474,
       20498.60421575, 22502.35495488, 23908.80035869, 23591.94166334,
       23414.89414342, 23523.55343152, 24909.55851466, 20664.13432228])

array([23704.81002932, 30277.27684072, 30458.14476477, 27100.76251824,
       29305.94424352, 29271.74377522, 30564.46611158, 32488.70132303,
       31891.63701325, 28864.13446524, 28527.30777714, 28008.25329826,
       29385.09935579, 27515.65646898, 29385.70535471, 20820.5152298 ,
       29616.43753367, 23503.02366358, 26903.63283498, 28336.12994564,
       30943.36743991, 30399.45380919, 27299.66880466, 28092.23031433,
       29778.28630112, 21947.42347288, 30244.27344631, 29892.29625004,
       22168.12323771, 29472.76994361, 31473.35147542, 26678.40927254,
       28928.77392429, 26589.93466207, 30501.32238234, 23368.34278763])

array([20606.97389893, 22111.94543817, 20799.07810538, 24534.34312171,
       23063.04433725, 19835.22133061, 24804.99282839, 24698.00377154,
       23570.41992701, 22503.91601162, 25084.48429451, 22753.56392254,
       21536.87029222, 26733.22854225, 24354.13564396, 19067.36138564,
       23633.6088968 , 19342.61988069, 21535.70621425, 23361.62510231,
       23489.52878204, 25530.94187194, 23681.56372354, 24448.44074977,
       27122.1606065 , 18450.27431994, 25229.2843939 , 25282.41242746,
       23026.50563123, 26633.48620523, 27128.28429844, 22442.36688249,
       26103.20652102, 19061.81368403, 21128.46308503, 23241.71171372])

PandasIndex(Index([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35],
      dtype='int64', name='station'))

az.ess(idataC, var_names=gev_var_names, method="quantile",prob=0.95)

<xarray.Dataset>
Dimensions:  (station: 36)
Coordinates:
  * station  (station) int64 0 1 2 3 4 5 6 7 8 9 ... 27 28 29 30 31 32 33 34 35
Data variables:
    xi_i     (station) float64 2.48e+04 2.332e+04 ... 1.655e+04 2.071e+04
    mu_i     (station) float64 3.045e+04 3.24e+04 ... 3.218e+04 3.121e+04
    sigma_i  (station) float64 2.047e+04 2.172e+04 ... 1.408e+04 2.731e+04

array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35])

array([24799.30094244, 23320.59558907, 16141.76519133, 20480.28927247,
       24439.28752413, 22193.24078948, 19068.89756386, 16569.01117071,
       17000.8366081 , 33234.00672191, 23238.61837911, 21780.05626508,
       19493.28193926, 14164.5318046 , 17863.73216485, 16250.11941655,
       19706.09710587, 16197.18039489, 15909.11736064, 14637.98312047,
       19510.72232118, 25175.2677673 , 23286.81607183, 18219.26055575,
       18760.95113759, 19451.08282319, 16622.21216653, 28748.99892999,
       30477.2249332 , 17174.74201504, 21513.59073288, 16472.40386519,
       20100.32241561, 22823.45857767, 16546.55654458, 20710.21789363])

array([30453.80485017, 32397.86397235, 32364.31833287, 30680.55984654,
       32383.47738081, 33536.3533055 , 32221.64252795, 31985.65847391,
       31961.0871252 , 31425.99647703, 33644.24912307, 33752.2614337 ,
       29929.11368771, 32361.03173707, 32643.14461784, 31082.07111719,
       32993.95828846, 31736.81602949, 29588.80547114, 32501.06580441,
       29794.1762125 , 32130.35207769, 33976.78750243, 31868.68347893,
       30345.55414303, 30193.6978859 , 32287.63035434, 33435.74516027,
       32546.34687422, 31299.32425619, 32942.97977736, 32395.56046113,
       33106.68996801, 33566.36772103, 32178.5537167 , 31210.63832064])

array([20470.95881694, 21720.18942624, 13096.82359087, 24644.85802898,
       15820.08343273, 23060.46820398, 14088.50178724, 20198.17417372,
       14191.30836963, 12225.39318522, 28134.80824198, 25224.56086141,
       20331.34433507, 18302.65873143, 20841.05222503, 29099.57758285,
       22638.9353154 , 25989.46097912, 22881.09507492, 23462.01396068,
       10079.44436839, 21895.05120113, 23545.79207293, 19337.63275879,
       22560.51009512, 25304.91670211, 24124.3816463 , 21992.6298073 ,
       25280.6563978 , 23167.85586097, 13248.53796006, 26537.9583752 ,
       23653.58160645, 23498.09314736, 14075.73540281, 27308.11380027])

PandasIndex(Index([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35],
      dtype='int64', name='station'))

az.ess(idataH, var_names=gev_var_names, method="quantile",prob=0.95)

<xarray.Dataset>
Dimensions:  (station: 36)
Coordinates:
  * station  (station) int64 0 1 2 3 4 5 6 7 8 9 ... 27 28 29 30 31 32 33 34 35
Data variables:
    xi_i     (station) float64 3.802e+04 3.697e+04 ... 3.674e+04 3.858e+04
    mu_i     (station) float64 4.032e+04 3.799e+04 ... 3.996e+04 3.807e+04
    sigma_i  (station) float64 3.666e+04 3.248e+04 ... 3.107e+04 3.409e+04

array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35])

array([38015.07288081, 36973.32930876, 37094.42966464, 36904.95951884,
       37554.22604623, 38446.82795664, 36215.70102189, 37719.2157197 ,
       36987.08429833, 38350.61922446, 36904.91209602, 36359.15571572,
       39033.2578736 , 37877.70568155, 38945.08891829, 37674.65229786,
       36318.1347611 , 38200.42178051, 37306.7460213 , 38733.64479067,
       37821.82458154, 37041.06163876, 37430.62058721, 37310.25612935,
       38001.15325925, 38486.60172938, 36983.38668793, 36355.03377299,
       37219.36066694, 37922.00506032, 36406.13529962, 37729.03599484,
       38839.49406868, 35330.82642384, 36736.19452538, 38575.02960829])

array([40317.58750925, 37993.79904021, 37977.52760133, 39391.86616871,
       38919.87934011, 37983.53292244, 36852.77452067, 38976.46757745,
       38784.05662328, 39374.75164108, 38646.28067607, 39508.68194168,
       38709.56768194, 39691.09643179, 38492.12543811, 39335.50745969,
       37587.91331531, 38004.96394048, 39754.3513237 , 38076.92256333,
       38663.45944128, 37480.6894743 , 40669.5416685 , 39556.81682858,
       38430.55673031, 38600.49422343, 38236.39648265, 38459.46160192,
       39888.0657197 , 39094.72207321, 36777.72625595, 39403.5605605 ,
       37162.64001926, 38224.82207704, 39961.6467733 , 38073.92809785])

array([36658.78450725, 32480.8721865 , 29533.97967048, 35298.71184214,
       31059.43247194, 30387.75615206, 29971.67637031, 30712.04497356,
       31540.01384897, 33087.76022176, 34001.17772779, 35897.51958414,
       32510.07804197, 31696.35245001, 30736.4635252 , 36861.30821431,
       28529.61126528, 38154.79106026, 32153.28777313, 34716.45843925,
       32174.51950817, 29483.52798487, 35676.39989149, 34092.34494708,
       30644.75478952, 38933.64717475, 28379.48130594, 31892.34381114,
       37883.38659898, 33454.29488058, 33755.23904662, 34448.4437697 ,
       30097.57561145, 32278.06104781, 31071.95920132, 34092.35831035])

PandasIndex(Index([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35],
      dtype='int64', name='station'))

az.ess(idataL, var_names=gev_var_names, method="quantile",prob=0.95)

<xarray.Dataset>
Dimensions:  (station: 36)
Coordinates:
  * station  (station) int64 0 1 2 3 4 5 6 7 8 9 ... 27 28 29 30 31 32 33 34 35
Data variables:
    xi_i     (station) float64 3.737e+04 3.702e+04 ... 3.815e+04 3.709e+04
    mu_i     (station) float64 3.89e+04 3.671e+04 ... 3.713e+04 3.739e+04
    sigma_i  (station) float64 3.531e+04 3.323e+04 ... 3.307e+04 3.333e+04

array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35])

array([37371.43083328, 37017.162647  , 36600.64364305, 40122.9012317 ,
       37875.33637889, 38014.25683922, 38028.09238814, 36585.64938585,
       38025.85341274, 38561.37283877, 38805.69755827, 40035.2901855 ,
       38613.28588625, 38450.30013286, 37295.30165656, 37729.06320523,
       36911.65518357, 37414.61497391, 38731.22843655, 38146.18185386,
       38529.74931212, 36263.65706287, 37979.09732733, 38420.73525445,
       37697.40972474, 37508.97104934, 35482.30886164, 36723.56585917,
       38494.76629371, 37005.39559365, 38645.08221195, 38109.99987672,
       37875.93947254, 38790.15944333, 38151.62862088, 37094.72897345])

array([38895.34611687, 36707.55006183, 38206.27119939, 37694.30048725,
       35696.62199334, 35431.71323983, 38034.11096687, 37124.28192391,
       37901.677328  , 39819.80951363, 36664.84212807, 36049.19299685,
       35906.89082879, 39378.95472822, 38458.30715057, 39048.60895127,
       37402.99271323, 37817.76587874, 38967.29736458, 38963.08685766,
       37308.78286058, 38216.01145875, 36682.44399837, 39060.83156487,
       38715.77435517, 37434.50186838, 38396.72373859, 38629.82739382,
       38607.69631251, 37909.67919759, 36520.88210489, 35527.99566531,
       38111.95810745, 37310.36783279, 37127.88664455, 37394.94753361])

array([35309.56879529, 33227.50105967, 32353.13794735, 36042.55096756,
       33107.10447904, 34626.69138252, 33424.42378804, 32093.39146114,
       34066.26814993, 36193.54271556, 34305.95893412, 36759.45522571,
       34437.03030982, 33898.19251   , 32704.08406445, 35706.65774884,
       31992.49577314, 38192.51427718, 34510.14681876, 34887.92097667,
       33883.94169858, 32520.118654  , 34202.84445405, 34105.66792225,
       34138.2374863 , 35629.99368418, 31661.19657703, 33030.79457431,
       34129.69834843, 31101.14958897, 33331.40086572, 35104.82360255,
       34461.54043513, 34013.80204018, 33065.58026177, 33332.73378708])

PandasIndex(Index([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35],
      dtype='int64', name='station'))

az.ess(idataS, var_names=gev_var_names, method="quantile",prob=0.95)

<xarray.Dataset>
Dimensions:  (station: 36)
Coordinates:
  * station  (station) int64 0 1 2 3 4 5 6 7 8 9 ... 27 28 29 30 31 32 33 34 35
Data variables:
    xi_i     (station) float64 2.612e+04 2.346e+04 ... 1.86e+04 1.985e+04
    mu_i     (station) float64 3.189e+04 3.297e+04 ... 3.604e+04 3.136e+04
    sigma_i  (station) float64 1.645e+04 1.604e+04 ... 2.065e+04 2.75e+04

array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35])

array([26121.05195338, 23464.30251948, 17229.6037861 , 20107.19558825,
       25766.20450998, 22523.69476238, 18031.69770189, 16179.55930429,
       19373.37507734, 32023.69241881, 24041.34796057, 25153.98460648,
       18244.57116441, 15437.96269846, 16598.64233199, 15685.43632863,
       22083.40315332, 20347.68451872, 20709.03561859, 17718.9602273 ,
       18490.70965963, 26159.71322294, 22644.64739617, 20921.67481842,
       24712.25042284, 19479.14333802, 19972.98995726, 28763.84810962,
       29630.53648956, 19202.23786993, 24916.21118697, 19922.32794174,
       25058.69610382, 23890.44810322, 18599.94853746, 19854.06075745])

array([31892.21155256, 32971.38005317, 34340.15982315, 32678.53317069,
       34917.97851532, 32663.40689711, 33411.68375854, 34041.19518271,
       33500.88206972, 27563.80062633, 33233.66790525, 32733.58336504,
       33007.62792877, 32737.98355644, 32940.64007237, 31566.79269062,
       31881.74300956, 32848.75958911, 32537.28642696, 31865.05880126,
       32397.6434489 , 33162.34507025, 30608.16829148, 33412.66924509,
       32862.14689583, 30141.24298215, 33807.62641418, 34785.43655084,
       31628.36867718, 32104.48828698, 34249.8744589 , 32978.19477952,
       33523.74828084, 32496.79058329, 36037.2547413 , 31357.69803094])

array([16453.47904115, 16040.75552876, 19604.37944752, 22224.57027341,
       17439.82530475, 13544.64654776, 19669.27595007, 18257.58056224,
       18527.40184773, 21881.13376173, 29878.31231654, 16760.68206577,
       14932.94404313, 21398.86123755, 19497.54215982, 34545.99128251,
       15815.65319339, 17568.2737006 , 15581.31904534, 17428.26129881,
       20390.87151538, 18692.78308377, 18222.3243335 , 17823.06236873,
       21112.20698474, 14853.58586378, 19886.81232615, 19526.59394446,
       22811.90515729, 21020.52323758, 22285.81420024, 31030.26268506,
       22223.21049265, 28880.33950917, 20650.26319663, 27496.25010182])

PandasIndex(Index([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35],
      dtype='int64', name='station'))

az.ess(idataC, var_names=gev_var_names, method="tail",prob=[0.05,0.95])

<xarray.Dataset>
Dimensions:  (station: 36)
Coordinates:
  * station  (station) int64 0 1 2 3 4 5 6 7 8 9 ... 27 28 29 30 31 32 33 34 35
Data variables:
    xi_i     (station) float64 2.48e+04 2.109e+04 ... 1.655e+04 2.071e+04
    mu_i     (station) float64 3.045e+04 3.19e+04 ... 3.218e+04 3.121e+04
    sigma_i  (station) float64 2.047e+04 2.172e+04 ... 1.408e+04 2.731e+04

array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35])

array([24799.30094244, 21092.91868618, 16141.76519133, 20241.00514722,
       21759.71591478, 22193.24078948, 19068.89756386, 16569.01117071,
       17000.8366081 , 17815.13723765, 19647.68851673, 21780.05626508,
       19493.28193926, 14164.5318046 , 17863.73216485, 16250.11941655,
       19706.09710587, 16197.18039489, 15909.11736064, 14637.98312047,
       19510.72232118, 21484.89384089, 23286.81607183, 18219.26055575,
       18760.95113759, 19451.08282319, 16622.21216653, 20177.72280025,
       19379.40027516, 17174.74201504, 20200.64689656, 16472.40386519,
       20100.32241561, 20467.19616383, 16546.55654458, 20710.21789363])

array([30453.80485017, 31899.23563843, 32174.25733671, 30680.55984654,
       30927.25681106, 31913.71190921, 31325.27437698, 31543.35262286,
       31610.60168839, 28825.93075944, 31424.56746446, 33662.73109742,
       29929.11368771, 31477.61113167, 32643.14461784, 30605.64871685,
       29472.71248685, 30192.2413916 , 29588.80547114, 32501.06580441,
       29794.1762125 , 31158.07593289, 33053.68912129, 31868.68347893,
       30210.1589984 , 30193.6978859 , 32245.10950928, 32183.22084865,
       31589.69770895, 31299.32425619, 32527.38096512, 32395.56046113,
       31687.38016649, 31637.89353154, 32178.5537167 , 31210.63832064])

array([20470.95881694, 21720.18942624, 13096.82359087, 24644.85802898,
       15820.08343273, 23060.46820398, 14088.50178724, 20198.17417372,
       14191.30836963, 12225.39318522, 28134.80824198, 25224.56086141,
       20331.34433507, 18302.65873143, 20841.05222503, 29099.57758285,
       22638.9353154 , 25989.46097912, 22881.09507492, 23462.01396068,
       10079.44436839, 21895.05120113, 23545.79207293, 19337.63275879,
       22560.51009512, 25304.91670211, 24124.3816463 , 21992.6298073 ,
       25280.6563978 , 23167.85586097, 13248.53796006, 26537.9583752 ,
       23653.58160645, 23498.09314736, 14075.73540281, 27308.11380027])

PandasIndex(Index([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35],
      dtype='int64', name='station'))

az.ess(idataH, var_names=gev_var_names, method="tail",prob=[0.05,0.95])

<xarray.Dataset>
Dimensions:  (station: 36)
Coordinates:
  * station  (station) int64 0 1 2 3 4 5 6 7 8 9 ... 27 28 29 30 31 32 33 34 35
Data variables:
    xi_i     (station) float64 3.646e+04 3.667e+04 ... 3.674e+04 3.692e+04
    mu_i     (station) float64 3.878e+04 3.799e+04 ... 3.718e+04 3.807e+04
    sigma_i  (station) float64 3.621e+04 3.248e+04 ... 3.107e+04 3.409e+04

array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35])

array([36455.01829253, 36665.71750592, 36391.39778609, 36616.91935493,
       35419.62364904, 37955.84194157, 36215.70102189, 37719.2157197 ,
       36717.99613794, 38350.61922446, 36904.91209602, 36359.15571572,
       38295.24518043, 37877.70568155, 38945.08891829, 37674.65229786,
       34133.69292406, 37294.48512817, 37306.7460213 , 38733.64479067,
       37792.13402973, 35541.651017  , 37430.62058721, 37310.25612935,
       38001.15325925, 37196.87403398, 36036.25461156, 36355.03377299,
       35232.66062507, 36999.8561165 , 36406.13529962, 37729.03599484,
       34527.86848095, 35330.82642384, 36736.19452538, 36915.60518877])

array([38781.84706656, 37993.79904021, 37289.68978635, 39391.86616871,
       36755.18393328, 36620.7048481 , 36852.77452067, 37985.02793606,
       37874.07146425, 39374.75164108, 38220.42859706, 37931.58881169,
       38116.66748349, 39287.79218616, 37594.09799205, 39335.50745969,
       37587.91331531, 38004.96394048, 39102.80651046, 38076.92256333,
       37004.29900612, 37480.6894743 , 36138.52505612, 37423.39317931,
       38430.55673031, 38570.81330199, 38236.39648265, 38283.47539945,
       38654.35880651, 36020.6741093 , 36777.72625595, 39403.5605605 ,
       37162.64001926, 38224.82207704, 37182.95134027, 38073.92809785])

array([36214.55987326, 32480.8721865 , 29533.97967048, 35298.71184214,
       31059.43247194, 30387.75615206, 29971.67637031, 30712.04497356,
       31540.01384897, 33087.76022176, 34001.17772779, 35897.51958414,
       32510.07804197, 31696.35245001, 30736.4635252 , 36464.21003228,
       28529.61126528, 38154.79106026, 32153.28777313, 34716.45843925,
       32174.51950817, 29483.52798487, 35076.00210501, 34092.34494708,
       30644.75478952, 38933.64717475, 28379.48130594, 31892.34381114,
       37883.38659898, 33454.29488058, 33755.23904662, 34448.4437697 ,
       30097.57561145, 32278.06104781, 31071.95920132, 34092.35831035])

PandasIndex(Index([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35],
      dtype='int64', name='station'))

az.ess(idataL, var_names=gev_var_names, method="tail",prob=[0.05,0.95])

<xarray.Dataset>
Dimensions:  (station: 36)
Coordinates:
  * station  (station) int64 0 1 2 3 4 5 6 7 8 9 ... 27 28 29 30 31 32 33 34 35
Data variables:
    xi_i     (station) float64 3.737e+04 3.609e+04 ... 3.716e+04 3.709e+04
    mu_i     (station) float64 3.776e+04 3.671e+04 ... 3.64e+04 3.631e+04
    sigma_i  (station) float64 3.506e+04 3.323e+04 ... 3.307e+04 3.333e+04

array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35])

array([37371.43083328, 36091.81632629, 36600.64364305, 38522.10542072,
       35687.97938986, 38014.25683922, 36448.99186991, 36585.64938585,
       38025.85341274, 38078.23935319, 38805.69755827, 39066.87924206,
       37752.178567  , 38450.30013286, 37295.30165656, 37729.06320523,
       35805.5488406 , 37414.61497391, 38731.22843655, 38146.18185386,
       37881.06986809, 36263.65706287, 37979.09732733, 38420.73525445,
       37697.40972474, 37508.97104934, 35482.30886164, 36723.56585917,
       36062.13765043, 37005.39559365, 38645.08221195, 38109.99987672,
       36221.39661202, 37995.78262843, 37158.79259589, 37094.72897345])

array([37758.40224835, 36707.55006183, 35913.83872914, 37405.70119188,
       35696.62199334, 35431.71323983, 36428.78298845, 35887.03419627,
       36829.07175566, 38354.54596095, 36664.84212807, 36049.19299685,
       35024.04163117, 37961.60740354, 36960.42208849, 38998.03430627,
       37402.99271323, 37817.76587874, 38967.29736458, 36969.35355698,
       36555.57712743, 38216.01145875, 36587.35148274, 37358.60840012,
       37048.51096773, 36618.43674894, 36586.0302687 , 36851.70156142,
       38153.82164679, 35533.55778107, 36520.88210489, 35527.99566531,
       37556.71272661, 37310.36783279, 36404.42467801, 36305.25722436])

array([35061.70356562, 33227.50105967, 32353.13794735, 36042.55096756,
       33107.10447904, 34626.69138252, 33424.42378804, 30142.39816886,
       34066.26814993, 36193.54271556, 34305.95893412, 36635.42521581,
       34437.03030982, 33898.19251   , 32161.74313106, 35706.65774884,
       31275.40055373, 38192.51427718, 34510.14681876, 34887.92097667,
       33883.94169858, 32520.118654  , 34202.84445405, 34105.66792225,
       33867.01797791, 35629.99368418, 31661.19657703, 33030.79457431,
       34129.69834843, 31101.14958897, 33331.40086572, 35104.82360255,
       34461.54043513, 33677.94213744, 33065.58026177, 33332.73378708])

PandasIndex(Index([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35],
      dtype='int64', name='station'))

az.ess(idataS, var_names=gev_var_names, method="tail",prob=[0.05,0.95])

<xarray.Dataset>
Dimensions:  (station: 36)
Coordinates:
  * station  (station) int64 0 1 2 3 4 5 6 7 8 9 ... 27 28 29 30 31 32 33 34 35
Data variables:
    xi_i     (station) float64 2.612e+04 2.346e+04 ... 1.86e+04 1.985e+04
    mu_i     (station) float64 3.189e+04 3.297e+04 ... 3.413e+04 3.136e+04
    sigma_i  (station) float64 1.645e+04 1.604e+04 ... 2.065e+04 2.75e+04

array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35])

array([26121.05195338, 23464.30251948, 17229.6037861 , 20107.19558825,
       21106.8684779 , 22523.69476238, 18031.69770189, 16179.55930429,
       19373.37507734, 23437.17454539, 22506.3810442 , 23219.07715481,
       18244.57116441, 15437.96269846, 16598.64233199, 15685.43632863,
       22083.40315332, 20347.68451872, 20709.03561859, 17718.9602273 ,
       18490.70965963, 19749.10175504, 22644.64739617, 20921.67481842,
       23498.05761398, 19479.14333802, 19972.98995726, 21344.45088538,
       19540.60336361, 19202.23786993, 24916.21118697, 19922.32794174,
       25058.69610382, 23890.44810322, 18599.94853746, 19854.06075745])

array([31892.21155256, 32971.38005317, 29239.74790394, 32050.13023206,
       31996.49964038, 32435.09145176, 33411.68375854, 32846.45737019,
       33500.88206972, 27563.80062633, 33101.09835421, 32733.58336504,
       33007.62792877, 32737.98355644, 32230.63905677, 30201.91381242,
       28732.9842816 , 32647.7289973 , 32537.28642696, 31865.05880126,
       32394.59011948, 32017.42250847, 30608.16829148, 33145.60424416,
       31939.63163666, 30141.24298215, 33807.62641418, 34383.96695474,
       31138.21181465, 32104.48828698, 34249.8744589 , 32978.19477952,
       33181.05218702, 32496.79058329, 34128.55796533, 31357.69803094])

array([16453.47904115, 16040.75552876, 19604.37944752, 22224.57027341,
       17439.82530475, 13544.64654776, 19669.27595007, 18257.58056224,
       18527.40184773, 21881.13376173, 23797.73975052, 16760.68206577,
       14932.94404313, 21398.86123755, 19497.54215982, 20673.91917996,
       15815.65319339, 17568.2737006 , 15581.31904534, 17428.26129881,
       20390.87151538, 18692.78308377, 18222.3243335 , 17823.06236873,
       21112.20698474, 14853.58586378, 19886.81232615, 19526.59394446,
       22811.90515729, 21020.52323758, 22285.81420024, 23968.61852614,
       22223.21049265, 15952.14244593, 20650.26319663, 27496.25010182])

PandasIndex(Index([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35],
      dtype='int64', name='station'))

#model_comp = {"Common": idataC, "Separate":idataS }
model_comp_int = {"Hilbert": idataH, "Latent":idataL }

with model_H:
    model_H.add_coord("test_station",test_longer_sites)
    _XYtest = pm.Data("XYtest",XY_test_longer_array, dims=["test_station", "feature"])
    # model_H.add_coord("test_station",LOO_sites) 
    # _XYtest = pm.Data("XYtest",XY_Loo, dims=["test_station", "feature"])
   
    #Hilbert space GP conditional 
    c_f_mu = hsgp_mu.conditional("c_f_mu",Xnew=_XYtest) #Hilbert
    c_f_sigma = hsgp_sigma.conditional("c_f_sigma", Xnew=_XYtest) #Hilbert
    c_f_xi = hsgp_xi.conditional("c_f_xi", Xnew=_XYtest) #Hilbert
    
    #GEV parameters conditioned on the traning data
    c_mu_i = pm.Deterministic("c_mu_i", at.abs(c_f_mu)+eps_mu, dims="test_station")
    c_sigma_i = pm.Deterministic("c_sigma_i", at.abs(c_f_sigma)+eps_sigma, dims="test_station")
    c_xi_i = pm.Deterministic("c_xi_i",  c_f_xi+eps_xi, dims="test_station")


    c_Ys = pmx.GenExtreme("c_Ys", mu=c_mu_i, sigma=c_sigma_i, xi=c_xi_i, scipy=False, dims=("years", "test_station"))
    #Return level
    c_Zs_5 = pm.Deterministic("c_Zs_5", c_mu_i - c_sigma_i / c_xi_i *  (1 -(-np.log(1-p5))**(-c_xi_i)), dims="test_station")
    c_Zs_10 = pm.Deterministic("c_Zs_10", c_mu_i - c_sigma_i / c_xi_i *  (1 -(-np.log(1-p10))**(-c_xi_i)), dims="test_station")
    c_Zs_20 = pm.Deterministic("c_Zs_20", c_mu_i - c_sigma_i / c_xi_i *  (1 -(-np.log(1-p20))**(-c_xi_i)), dims="test_station")
    c_Zs_30 = pm.Deterministic("c_Zs_30", c_mu_i - c_sigma_i / c_xi_i *  (1 -(-np.log(1-p30))**(-c_xi_i)), dims="test_station")
    c_Zs_50 = pm.Deterministic("c_Zs_50", c_mu_i - c_sigma_i / c_xi_i *  (1 -(-np.log(1-p50))**(-c_xi_i)), dims="test_station")
    c_Zs_100 = pm.Deterministic("c_Zs_100", c_mu_i - c_sigma_i / c_xi_i * (1 -(-np.log(1-p100))**(-c_xi_i)), dims="test_station")
    c_Zs_200 = pm.Deterministic("c_Zs_200", c_mu_i - c_sigma_i / c_xi_i *  (1 -(-np.log(1-p200))**(-c_xi_i)), dims="test_station")
    c_Zs_500 = pm.Deterministic("c_Zs_500", c_mu_i - c_sigma_i / c_xi_i * (1 -(-np.log(1-p500))**(-c_xi_i)), dims="test_station")
    c_Zs_1000 = pm.Deterministic("c_Zs_1000", c_mu_i - c_sigma_i / c_xi_i * (1 -(-np.log(1-p1000))**(-c_xi_i)), dims="test_station")
    c_Zs_5000 = pm.Deterministic("c_Zs_5000", c_mu_i - c_sigma_i / c_xi_i * (1 -(-np.log(1-p5000))**(-c_xi_i)), dims="test_station")
    c_Zs_10000 = pm.Deterministic("c_Zs_10000", c_mu_i - c_sigma_i / c_xi_i * (1 -(-np.log(1-p10000))**(-c_xi_i)), dims="test_station")
    c_Zs_50000 = pm.Deterministic("c_Zs_50000", c_mu_i - c_sigma_i / c_xi_i * (1 -(-np.log(1-p50000))**(-c_xi_i)), dims="test_station")
    c_Zs_100000 = pm.Deterministic("c_Zs_100000", c_mu_i - c_sigma_i / c_xi_i * (1 -(-np.log(1-p100000))**(-c_xi_i)), dims="test_station")

with model_L:
    model_L.add_coord("test_station",test_longer_sites) 
    _XYtest = pm.Data("XYtest",XY_test_longer_array, dims=["test_station", "feature"]) 
    #model_L.add_coord("test_station",LOO_sites) 
    #_XYtest = pm.Data("XYtest",XY_Loo, dims=["test_station", "feature"])
    #GP conditional 
    c_f_mu = gp_mu.conditional("c_f_mu",Xnew=_XYtest) #Hilbert
    c_f_sigma = gp_sigma.conditional("c_f_sigma", Xnew=_XYtest) #Hilbert
    c_f_xi = gp_xi.conditional("c_f_xi", Xnew=_XYtest) #Hilbert

    #GEV parameters conditioned on the traning data
    c_mu_i = pm.Deterministic("c_mu_i", at.abs(c_f_mu)+eps_mu, dims="test_station")
    c_sigma_i = pm.Deterministic("c_sigma_i", at.abs(c_f_sigma)+eps_sigma, dims="test_station")
    c_xi_i = pm.Deterministic("c_xi_i",  c_f_xi+eps_xi, dims="test_station")


    c_Ys = pmx.GenExtreme("c_Ys", mu=c_mu_i, sigma=c_sigma_i, xi=c_xi_i, scipy=False, dims=("years", "test_station"))
    #Return level
    c_Zs_5 = pm.Deterministic("c_Zs_5", c_mu_i - c_sigma_i / c_xi_i *  (1 -(-np.log(1-p5))**(-c_xi_i)), dims="test_station")
    c_Zs_10 = pm.Deterministic("c_Zs_10", c_mu_i - c_sigma_i / c_xi_i *  (1 -(-np.log(1-p10))**(-c_xi_i)), dims="test_station")
    c_Zs_20 = pm.Deterministic("c_Zs_20", c_mu_i - c_sigma_i / c_xi_i *  (1 -(-np.log(1-p20))**(-c_xi_i)), dims="test_station")
    c_Zs_30 = pm.Deterministic("c_Zs_30", c_mu_i - c_sigma_i / c_xi_i *  (1 -(-np.log(1-p30))**(-c_xi_i)), dims="test_station")
    c_Zs_50 = pm.Deterministic("c_Zs_50", c_mu_i - c_sigma_i / c_xi_i *  (1 -(-np.log(1-p50))**(-c_xi_i)), dims="test_station")
    c_Zs_100 = pm.Deterministic("c_Zs_100", c_mu_i - c_sigma_i / c_xi_i * (1 -(-np.log(1-p100))**(-c_xi_i)), dims="test_station")
    c_Zs_200 = pm.Deterministic("c_Zs_200", c_mu_i - c_sigma_i / c_xi_i *  (1 -(-np.log(1-p200))**(-c_xi_i)), dims="test_station")
    c_Zs_500 = pm.Deterministic("c_Zs_500", c_mu_i - c_sigma_i / c_xi_i * (1 -(-np.log(1-p500))**(-c_xi_i)), dims="test_station")
    c_Zs_1000 = pm.Deterministic("c_Zs_1000", c_mu_i - c_sigma_i / c_xi_i * (1 -(-np.log(1-p1000))**(-c_xi_i)), dims="test_station")
    c_Zs_5000 = pm.Deterministic("c_Zs_5000", c_mu_i - c_sigma_i / c_xi_i * (1 -(-np.log(1-p5000))**(-c_xi_i)), dims="test_station")
    c_Zs_10000 = pm.Deterministic("c_Zs_10000", c_mu_i - c_sigma_i / c_xi_i * (1 -(-np.log(1-p10000))**(-c_xi_i)), dims="test_station")
    c_Zs_50000 = pm.Deterministic("c_Zs_50000", c_mu_i - c_sigma_i / c_xi_i * (1 -(-np.log(1-p50000))**(-c_xi_i)), dims="test_station")
    c_Zs_100000 = pm.Deterministic("c_Zs_100000", c_mu_i - c_sigma_i / c_xi_i * (1 -(-np.log(1-p100000))**(-c_xi_i)), dims="test_station")

with model_H:
    pm.set_data({"sites": test_longer_sites})
    #pm.set_data({"sites": LOO_sites}) 
    out_of_sample_H = pm.sample_posterior_predictive(traceH, 
                                                      var_names=["c_f_mu", "c_f_sigma", "c_f_xi", "c_xi_i","c_mu_i","c_sigma_i", "c_Ys","c_Zs_5", "c_Zs_10", 
                                                              "c_Zs_20","c_Zs_30",
                                                              "c_Zs_50","c_Zs_100","c_Zs_200", "c_Zs_500", "c_Zs_1000","c_Zs_5000", "c_Zs_10000", "c_Zs_50000", "c_Zs_100000"]
                                                   , predictions=True, random_seed=rng)

Sampling: [c_Ys]

with model_L:
    pm.set_data({"sites":test_longer_sites}) 
    #pm.set_data({"sites": LOO_sites}) 
    out_of_sample_L = pm.sample_posterior_predictive(traceL, 
                                                     var_names=["c_f_mu", "c_f_sigma", "c_f_xi", "c_xi_i","c_mu_i","c_sigma_i", "c_Ys","c_Zs_5", "c_Zs_10", 
                                                              "c_Zs_20","c_Zs_30",
                                                              "c_Zs_50","c_Zs_100","c_Zs_200", "c_Zs_500", "c_Zs_1000","c_Zs_5000", "c_Zs_10000", "c_Zs_50000", "c_Zs_100000"]
                                                   , predictions=True, random_seed=rng)

Sampling: [c_Ys, c_f_mu, c_f_sigma, c_f_xi]

all_idata = [idataC, idataH,idataL, idataS]
models_names = ["Common", "Hilbert", "Latent", "Separate"]

models_intrp_in_out_sample_dict = {"Hilbert": [idataH, out_of_sample_H], "Latent": [idataL, out_of_sample_L]}
models_out_sample_dict = {"Hilbert": out_of_sample_H, "Latent": out_of_sample_L}

az.plot_density(all_idata, var_names=["xi_i"], hdi_prob=.95,  figsize=(22,25), data_labels=["Common","Hilbert", "Latent","Separate"]);

az.plot_density(all_idata, var_names=["mu_i"], hdi_prob=.95,  figsize=(22,25), data_labels=["Common", "Hilbert", "Latent","Separate"]);

az.plot_density(all_idata, var_names=["sigma_i"], hdi_prob=.95,  figsize=(22,25), data_labels=["Common","Hilbert", "Latent","Separate"]);

BHM_boxplot_models_comparison(models_names, all_idata, train_sites, "Zs_50",  region_BM, train, n_draws=10000)

BHM_boxplot_models_comparison(models_names, all_idata, train_sites, "Zs_100", region_BM, train, n_draws=10000)

BHM_boxplot_models_comparison(models_names, all_idata, train_sites, "Zs_200", region_BM, train, n_draws=10000)

BHM_boxplot_models_comparison(models_names, all_idata, train_sites, "Zs_500", region_BM, train, n_draws=10000)

BHM_boxplot_models_comparison(models_names, all_idata, train_sites, "Zs_1000", region_BM, train, n_draws=10000)

models_intrp_names = ["Hilbert", "Latent"]
models_intrp_idata = [idataH,idataL]

models_no_intrp_names = ["Common", "Separate"]
models_no_intrp_idata = [idataC,idataS]

BHM_violinplot_models_comparison_old(models_names=models_intrp_names, models_idata=models_intrp_idata, train_sites=train_sites, var_name="xi_i", showmeans=False, showextrema=False, showmedians=True,points=40000)

#BHM_violinplot_models_comparison(models_names=models_intrp_names, models_idata=models_intrp_idata, train_sites=train_sites, var_name="xi_i", train=train, region_BM=region_BM,  Baseline_region_BM= Baseline_region_BM,showmeans=False, showextrema=False, showmedians=True,points=40000, n_draws=10000)

BHM_violinplot_models_comparison_old(models_names=models_no_intrp_names, models_idata=models_no_intrp_idata, train_sites=train_sites, var_name="xi_i", showmeans=False, showextrema=False, showmedians=True,points=40000)

#BHM_violinplot_models_comparison(models_names=models_no_intrp_names, models_idata=models_no_intrp_idata, train_sites=train_sites, var_name="xi_i", train=train, region_BM=region_BM,  Baseline_region_BM= Baseline_region_BM,showmeans=False, showextrema=False, showmedians=True,points=40000, n_draws=10000)

BHM_violinplot_models_comparison_old(models_names=models_intrp_names, models_idata=models_intrp_idata, train_sites=train_sites, var_name="mu_i", showmeans=False, showextrema=False, showmedians=True,points=40000)

#BHM_violinplot_models_comparison(models_names=models_intrp_names, models_idata=models_intrp_idata, train_sites=train_sites, var_name="mu_i", train=train, region_BM=region_BM,  Baseline_region_BM= Baseline_region_BM,showmeans=False, showextrema=False, showmedians=True,points=40000, n_draws=10000)

BHM_violinplot_models_comparison_old(models_names=models_no_intrp_names, models_idata=models_no_intrp_idata, train_sites=train_sites, var_name="mu_i", showmeans=False, showextrema=False, showmedians=True,points=40000)

#BHM_violinplot_models_comparison(models_names=models_no_intrp_names, models_idata=models_no_intrp_idata, train_sites=train_sites, var_name="mu_i", train=train, region_BM=region_BM,  Baseline_region_BM= Baseline_region_BM,showmeans=False, showextrema=False, showmedians=True,points=40000, n_draws=10000)

BHM_violinplot_models_comparison_old(models_names=models_intrp_names, models_idata=models_intrp_idata, train_sites=train_sites, var_name="sigma_i", showmeans=False, showextrema=False, showmedians=True,points=40000)

#BHM_violinplot_models_comparison(models_names=models_intrp_names, models_idata=models_intrp_idata, train_sites=train_sites, var_name="sigma_i", train=train, region_BM=region_BM,  Baseline_region_BM= Baseline_region_BM,showmeans=False, showextrema=False, showmedians=True,points=40000, n_draws=10000)

BHM_violinplot_models_comparison_old(models_names=models_no_intrp_names, models_idata=models_no_intrp_idata, train_sites=train_sites, var_name="sigma_i", showmeans=False, showextrema=False, showmedians=True,points=40000)

#BHM_violinplot_models_comparison(models_names=models_no_intrp_names, models_idata=models_no_intrp_idata, train_sites=train_sites, var_name="sigma_i", train=train, region_BM=region_BM,  Baseline_region_BM= Baseline_region_BM,showmeans=False, showextrema=False, showmedians=True,points=40000, n_draws=10000)

BHM_geographical_spread_GEV_param_posterior_plot(models_intrp_idata[0] ,train, Baseline_region_BM, compare_stations_in_top_4_dict, models_intrp_names[0])

BHM_geographical_spread_GEV_param_posterior_plot(models_intrp_idata[1] ,train, Baseline_region_BM, compare_stations_in_top_4_dict, models_intrp_names[1])

BHM_geographical_spread_GEV_param_posterior_plot(models_intrp_idata[0] ,train, Baseline_region_BM, geographical_spread_list[0], models_intrp_names[0])

BHM_geographical_spread_GEV_param_posterior_plot(models_intrp_idata[1] ,train, Baseline_region_BM, geographical_spread_list[0], models_intrp_names[1])

BHM_geographical_spread_GEV_param_posterior_plot(models_intrp_idata[0] ,train, Baseline_region_BM, geographical_spread_list[1], models_intrp_names[0])

BHM_geographical_spread_GEV_param_posterior_plot(models_intrp_idata[1] ,train, Baseline_region_BM, geographical_spread_list[1], models_intrp_names[1])

BHM_geographical_spread_GEV_param_posterior_plot(models_intrp_idata[0] ,train, Baseline_region_BM, geographical_spread_list[2], models_intrp_names[0])

BHM_geographical_spread_GEV_param_posterior_plot(models_intrp_idata[1] ,train, Baseline_region_BM, geographical_spread_list[2], models_intrp_names[1])

BHM_geographical_spread_GEV_param_posterior_plot(models_intrp_idata[0] ,train, Baseline_region_BM, geographical_spread_list[3], models_intrp_names[0])

BHM_geographical_spread_GEV_param_posterior_plot(models_intrp_idata[1] ,train, Baseline_region_BM, geographical_spread_list[3], models_intrp_names[1])

BHM_geographical_spread_GEV_param_posterior_plot(models_intrp_idata[0] ,train, Baseline_region_BM, geographical_spread_list[4], models_intrp_names[0])

BHM_geographical_spread_GEV_param_posterior_plot(models_intrp_idata[1] ,train, Baseline_region_BM, geographical_spread_list[4], models_intrp_names[1])

BHM_geographical_spread_GEV_param_posterior_plot(models_intrp_idata[0] ,train, Baseline_region_BM, geographical_spread_list[5], models_intrp_names[0])

BHM_geographical_spread_GEV_param_posterior_plot(models_intrp_idata[1] ,train, Baseline_region_BM, geographical_spread_list[5], models_intrp_names[1])

compare_stations_in_top_4_dict

{'station': [28, 15, 17, 31], 'location': 'Top 4 (fewest NaNs) stations'}

BHM_supp_GEV_KDE_plot(models_intrp_idata[1], "post", geographical_spread_list[0], train, region_BM,  models_intrp_names[1], region_sites, Baseline_region_BM);

BHM_supp_GEV_KDE_plot(models_intrp_idata[0], "post", geographical_spread_list[0], train, region_BM,  models_intrp_names[0], region_sites, Baseline_region_BM);

BHM_supp_GEV_KDE_plot(models_intrp_idata[1], "post", geographical_spread_list[0], train, region_BM,  models_intrp_names[1], region_sites, Baseline_region_BM);

BHM_supp_GEV_KDE_plot(models_intrp_idata[0], "post", geographical_spread_list[1], train, region_BM,  models_intrp_names[0], region_sites, Baseline_region_BM);

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.

BHM_supp_GEV_KDE_plot(models_intrp_idata[1], "post", geographical_spread_list[1], train, region_BM,  models_intrp_names[1], region_sites, Baseline_region_BM);

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.

BHM_supp_GEV_KDE_plot(models_intrp_idata[0], "post", geographical_spread_list[2], train, region_BM,  models_intrp_names[0], region_sites, Baseline_region_BM);

BHM_supp_GEV_KDE_plot(models_intrp_idata[1], "post", geographical_spread_list[2], train, region_BM,  models_intrp_names[1], region_sites, Baseline_region_BM);

BHM_supp_GEV_KDE_plot(models_intrp_idata[0], "post", geographical_spread_list[3], train, region_BM,  models_intrp_names[0], region_sites, Baseline_region_BM);

BHM_supp_GEV_KDE_plot(models_intrp_idata[1], "post", geographical_spread_list[3], train, region_BM,  models_intrp_names[1], region_sites, Baseline_region_BM);

BHM_supp_GEV_KDE_plot(models_intrp_idata[0], "post", geographical_spread_list[4], train, region_BM,  models_intrp_names[0], region_sites, Baseline_region_BM);

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.

BHM_supp_GEV_KDE_plot(models_intrp_idata[1], "post", geographical_spread_list[4], train, region_BM,  models_intrp_names[1], region_sites, Baseline_region_BM);

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.

BHM_supp_GEV_KDE_plot(models_intrp_idata[0], "post", geographical_spread_list[5], train, region_BM,  models_intrp_names[0], region_sites, Baseline_region_BM);

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.

BHM_supp_GEV_KDE_plot(models_intrp_idata[1], "post", geographical_spread_list[5], train, region_BM,  models_intrp_names[1], region_sites, Baseline_region_BM);

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.

BHM_hdpi_and_MLM_ci_return_lvl_multiple_sites_plot(models_intrp_idata[0], compare_stations_in_top_4_dict, train=train, region_BM=region_base_BM,model_name=models_intrp_names[0], region_sites=region_sites, n_draws=10000,  figsize=(22,22), show=True)

BHM_hdpi_and_MLM_ci_return_lvl_multiple_sites_plot(models_intrp_idata[1], compare_stations_in_top_4_dict, train, region_base_BM,models_intrp_names[1], region_sites, n_draws=10000,  figsize=(22,22), show=True)

BHM_hdpi_and_MLM_ci_return_lvl_multiple_sites_plot(models_intrp_idata[0], geographical_spread_list[0], train, region_base_BM,models_intrp_names[0], region_sites, n_draws=10000 ,figsize=(22,22), show=True)

BHM_hdpi_and_MLM_ci_return_lvl_multiple_sites_plot(models_intrp_idata[1],geographical_spread_list[0], train, region_base_BM,models_intrp_names[1], region_sites, n_draws=10000 ,figsize=(22,22), show=True)

BHM_hdpi_and_MLM_ci_return_lvl_multiple_sites_plot(models_intrp_idata[0], geographical_spread_list[1], train, region_base_BM,models_intrp_names[0], region_sites, n_draws=10000 ,figsize=(22,22), show=True)

BHM_hdpi_and_MLM_ci_return_lvl_multiple_sites_plot(models_intrp_idata[1], geographical_spread_list[1], train, region_base_BM,models_intrp_names[1], region_sites, n_draws=10000 ,figsize=(22,22), show=True)

BHM_hdpi_and_MLM_ci_return_lvl_multiple_sites_plot(models_intrp_idata[0], geographical_spread_list[2], train, region_base_BM,models_intrp_names[0], region_sites, n_draws=10000 ,figsize=(22,22), show=True)

BHM_hdpi_and_MLM_ci_return_lvl_multiple_sites_plot(models_intrp_idata[1], geographical_spread_list[2], train, region_base_BM,models_intrp_names[1], region_sites, n_draws=10000 ,figsize=(22,22), show=True)

BHM_hdpi_and_MLM_ci_return_lvl_multiple_sites_plot(models_intrp_idata[0], geographical_spread_list[3], train, region_base_BM,models_intrp_names[0], region_sites, n_draws=10000  ,figsize=(22,22), show=True)

BHM_hdpi_and_MLM_ci_return_lvl_multiple_sites_plot(models_intrp_idata[1], geographical_spread_list[3], train, region_base_BM,models_intrp_names[1], region_sites, n_draws=10000  ,figsize=(22,22), show=True)

BHM_hdpi_and_MLM_ci_return_lvl_multiple_sites_plot(models_intrp_idata[0], geographical_spread_list[4], train, region_base_BM,models_intrp_names[0], region_sites, n_draws=10000  ,figsize=(22,22), show=True)

BHM_hdpi_and_MLM_ci_return_lvl_multiple_sites_plot(models_intrp_idata[1], geographical_spread_list[4], train, region_base_BM,models_intrp_names[1], region_sites, n_draws=10000  ,figsize=(22,22), show=True)

BHM_hdpi_and_MLM_ci_return_lvl_multiple_sites_plot(models_intrp_idata[0], geographical_spread_list[5], train, region_base_BM,models_intrp_names[0], region_sites, n_draws=10000  ,figsize=(22,22), show=True)

BHM_hdpi_and_MLM_ci_return_lvl_multiple_sites_plot(models_intrp_idata[1], geographical_spread_list[5], train, region_base_BM,models_intrp_names[1], region_sites, n_draws=10000  ,figsize=(22,22), show=True)

BHM_subset_pair_plot(models_intrp_idata[0],geographical_spread_list[0], var_names="xi_i", model_name=models_intrp_names[0], figsize=(10,10))

BHM_subset_pair_plot(models_intrp_idata[1],geographical_spread_list[0], var_names="xi_i", model_name=models_intrp_names[1], color="Oranges", figsize=(15,15))

BHM_subset_pair_plot(models_intrp_idata[0],geographical_spread_list[1], var_names="xi_i", model_name=models_intrp_names[0], figsize=(10,10))

BHM_subset_pair_plot(models_intrp_idata[1],geographical_spread_list[1], var_names="xi_i", model_name=models_intrp_names[1], color="Oranges", figsize=(15,15))

BHM_subset_pair_plot(models_intrp_idata[0],geographical_spread_list[2], var_names="xi_i", model_name=models_intrp_names[0], figsize=(15,15))

BHM_subset_pair_plot(models_intrp_idata[1],geographical_spread_list[2], var_names="xi_i", model_name=models_intrp_names[1], color="Oranges", figsize=(15,15))

BHM_subset_pair_plot(models_intrp_idata[0],geographical_spread_list[3], var_names="xi_i", model_name=models_intrp_names[0], figsize=(10,10))

BHM_subset_pair_plot(models_intrp_idata[1],geographical_spread_list[3], var_names="xi_i", model_name=models_intrp_names[1], color="Oranges", figsize=(10,10))

BHM_subset_pair_plot(models_intrp_idata[0],geographical_spread_list[4], var_names="xi_i", model_name=models_intrp_names[0], figsize=(15,15))

BHM_subset_pair_plot(models_intrp_idata[1],geographical_spread_list[4], var_names="xi_i", model_name=models_intrp_names[1], color="Oranges", figsize=(15,15))

BHM_subset_pair_plot(models_intrp_idata[0],geographical_spread_list[5], var_names="xi_i", model_name=models_intrp_names[0], figsize=(15,15))

BHM_subset_pair_plot(models_intrp_idata[1],geographical_spread_list[5], var_names="xi_i", model_name=models_intrp_names[1], color="Oranges", figsize=(15,15))

BHM_subset_pair_plot(models_intrp_idata[0],geographical_spread_list[0], var_names="mu_i", model_name=models_intrp_names[0], figsize=(15,15))

BHM_subset_pair_plot(models_intrp_idata[1],geographical_spread_list[0], var_names="mu_i", model_name=models_intrp_names[1], color="Oranges", figsize=(15,15))

BHM_subset_pair_plot(models_intrp_idata[0],geographical_spread_list[1], var_names="mu_i", model_name=models_intrp_names[0], figsize=(15,15))

BHM_subset_pair_plot(models_intrp_idata[1],geographical_spread_list[1], var_names="mu_i", model_name=models_intrp_names[1], color="Oranges", figsize=(15,15))

BHM_subset_pair_plot(models_intrp_idata[0],geographical_spread_list[2], var_names="mu_i", model_name=models_intrp_names[0], figsize=(15,15))

BHM_subset_pair_plot(models_intrp_idata[1],geographical_spread_list[2], var_names="mu_i", model_name=models_intrp_names[1], color="Oranges", figsize=(15,15))

BHM_subset_pair_plot(models_intrp_idata[0],geographical_spread_list[3], var_names="mu_i", model_name=models_intrp_names[0], figsize=(10,10))

BHM_subset_pair_plot(models_intrp_idata[1],geographical_spread_list[3], var_names="mu_i", model_name=models_intrp_names[1], color="Oranges", figsize=(10,10))

BHM_subset_pair_plot(models_intrp_idata[0],geographical_spread_list[4], var_names="mu_i", model_name=models_intrp_names[0], figsize=(15,15))

BHM_subset_pair_plot(models_intrp_idata[1],geographical_spread_list[4], var_names="mu_i", model_name=models_intrp_names[1], color="Oranges", figsize=(15,15))

BHM_subset_pair_plot(models_intrp_idata[0],geographical_spread_list[5], var_names="mu_i", model_name=models_intrp_names[0], figsize=(10,10))

BHM_subset_pair_plot(models_intrp_idata[1],geographical_spread_list[5], var_names="mu_i", model_name=models_intrp_names[1], color="Oranges", figsize=(10,10))

BHM_subset_pair_plot(models_intrp_idata[0],geographical_spread_list[0], var_names="sigma_i", model_name=models_intrp_names[0], figsize=(10,10))

BHM_subset_pair_plot(models_intrp_idata[1],geographical_spread_list[0], var_names="sigma_i", model_name=models_intrp_names[1], color="Oranges", figsize=(15,15))

BHM_subset_pair_plot(models_intrp_idata[0],geographical_spread_list[1], var_names="sigma_i", model_name=models_intrp_names[0], figsize=(10,10))

BHM_subset_pair_plot(models_intrp_idata[1],geographical_spread_list[1], var_names="sigma_i", model_name=models_intrp_names[1], color="Oranges", figsize=(10,10))

BHM_subset_pair_plot(models_intrp_idata[0],geographical_spread_list[2], var_names="sigma_i", model_name=models_intrp_names[0], figsize=(15,15))

BHM_subset_pair_plot(models_intrp_idata[1],geographical_spread_list[2], var_names="sigma_i", model_name=models_intrp_names[1], color="Oranges", figsize=(15,15))

BHM_subset_pair_plot(models_intrp_idata[0],geographical_spread_list[3], var_names="sigma_i", model_name=models_intrp_names[0], figsize=(15,15))

BHM_subset_pair_plot(models_intrp_idata[1],geographical_spread_list[3], var_names="sigma_i", model_name=models_intrp_names[1], color="Oranges", figsize=(10,10))

BHM_subset_pair_plot(models_intrp_idata[0],geographical_spread_list[4], var_names="sigma_i", model_name=models_intrp_names[0], figsize=(15,15))

BHM_subset_pair_plot(models_intrp_idata[1],geographical_spread_list[4], var_names="sigma_i", model_name=models_intrp_names[1], color="Oranges", figsize=(15,15))

BHM_subset_pair_plot(models_intrp_idata[0],geographical_spread_list[5], var_names="sigma_i", model_name=models_intrp_names[0], figsize=(15,15))

BHM_subset_pair_plot(models_intrp_idata[1],geographical_spread_list[5], var_names="sigma_i", model_name=models_intrp_names[1], color="Oranges", figsize=(15,15))

axes = az.plot_forest([idataH, idataL],model_names= ["Hilbert",  "Latent"],
                            kind='forestplot',
                           var_names=["xi_i"],
                           #filter_vars="regex",
                             combined=True,
                      hdi_prob=.95,
                           figsize=(14, 14));
axes[0].set_title("Estimated shape param of the GEV");

axes = az.plot_forest([idataH, idataL],model_names= ["Hilbert",  "Latent"],
                            kind='forestplot',
                           var_names=["mu_i"],
                           #filter_vars="regex",
                             combined=True,
                      hdi_prob=.95,
                           figsize=(14, 14));
axes[0].set_title("Estimated location param of the GEV");

axes = az.plot_forest([idataH, idataL],model_names= ["Hilbert",  "Latent"],
                            kind='forestplot',
                           var_names=["sigma_i"],
                           #filter_vars="regex",
                             combined=True,
                      hdi_prob=.95,
                           figsize=(14, 14));
axes[0].set_title("Estimated scale param of the GEV");

#BHM_hpdi_and_MLM_ci_return_lvl_Loo_multiple_models_plot(models_out_sample_dict, test_sites, region_BM, Loo_column, region_sites)

#BHM_GEV_param_LOO_posterior_plot(models_out_sample_dict, Loo_column,region_BM, region_sites)

#LOO_pair_plot(out_of_sample_H, "Hilbert", region_sites, Loo_column)

#LOO_pair_plot(out_of_sample_L, "Latent", region_sites, Loo_column, colors="orange")

BHM_comparative_mean_and_hpdi_var_name_surface_plots(models_intrp_in_out_sample_dict, XY_train, XY_test_longer_array, X, Y, var_name="xi_i",figsize=(36,36))

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

BHM_comparative_mean_and_hpdi_var_name_surface_plots(models_intrp_in_out_sample_dict, XY_train, XY_test_longer_array, X, Y, var_name="xi_i",figsize=(36,36))

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

BHM_comparative_mean_and_hpdi_var_name_surface_plots(models_intrp_in_out_sample_dict, XY_train, XY_test_longer_array, X, Y, var_name="xi_i",figsize=(36,36))

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

BHM_comparative_mean_and_hpdi_var_name_surface_plots(models_intrp_in_out_sample_dict, XY_train, XY_test_longer_array, X, Y, var_name="mu_i",figsize=(36,36))

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

BHM_comparative_mean_and_hpdi_var_name_surface_plots(models_intrp_in_out_sample_dict, XY_train, XY_test_longer_array, X, Y, var_name="sigma_i",figsize=(36,36))

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

BHM_comparative_mean_and_hpdi_var_name_surface_plots(models_intrp_in_out_sample_dict, XY_train, XY_test_longer_array, X, Y, var_name="Zs_100", return_period=100,figsize=(36,36))

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

BHM_comparative_mean_and_hpdi_var_name_surface_plots(models_intrp_in_out_sample_dict, XY_train, XY_test_longer_array, X, Y, var_name="Zs_500", return_period=500,figsize=(36,36))

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

BHM_comparative_mean_and_hpdi_var_name_surface_plots(models_intrp_in_out_sample_dict, XY_train, XY_test_longer_array, X, Y, var_name="Zs_1000", return_period=1000,figsize=(36,36))

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

BHM_comparative_mean_and_hpdi_var_name_surface_plots(models_intrp_in_out_sample_dict, XY_train, XY_test_longer_array, X, Y, var_name="Zs_5000", return_period=5000,figsize=(36,36))

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

BHM_comparative_mean_and_hpdi_var_name_surface_plots(models_intrp_in_out_sample_dict, XY_train, XY_test_longer_array, X, Y, var_name="Zs_10000", return_period=10000,figsize=(36,36))

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

BHM_comparative_mean_and_hpdi_var_name_surface_plots(models_intrp_in_out_sample_dict, XY_train, XY_test_longer_array, X, Y, var_name="Zs_50000", return_period=50000, figsize=(36,36))

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

BHM_comparative_mean_and_hpdi_var_name_surface_plots(models_intrp_in_out_sample_dict, XY_train, XY_test_longer_array, X, Y, var_name="Zs_100000", return_period=100000, figsize=(36,36))

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "ok" (-> marker='o'). The keyword argument will take precedence.

/home/sm_seuba/DataAnalysisPlt.py:802: UserWarning:

color is redundantly defined by the 'color' keyword argument and the fmt string "ok" (-> color='k'). The keyword argument will take precedence.

	rank	elpd_loo	p_loo	elpd_diff	weight	se	dse	warning	scale
Hilbert	0	319.036085	39.267076	0.000000	0.513596	29.169579	0.000000	False	log
Latent	1	318.173279	35.403501	0.862807	0.323106	29.237279	2.249843	False	log
Common	2	310.592203	78.147082	8.443882	0.073388	29.973057	8.218994	True	log
Separate	3	309.916475	80.546255	9.119610	0.089909	30.821881	8.892290	True	log

	rank	elpd_loo	p_loo	elpd_diff	weight	se	dse	warning	scale
Hilbert	0	319.036085	39.267076	0.000000	3.174852e-03	30.049570	0.000000	False	log
Latent	1	318.173279	35.403501	0.862807	5.853058e-01	30.198486	2.249843	False	log
Common	2	310.592203	78.147082	8.443882	1.570696e-15	30.690284	8.218994	True	log
Separate	3	309.916475	80.546255	9.119610	4.115194e-01	31.552372	8.892290	True	log

	rank	elpd_loo	p_loo	elpd_diff	weight	se	dse	warning	scale
Hilbert	0	319.036085	39.267076	0.000000	0.667422	30.049570	0.000000	False	log
Latent	1	318.173279	35.403501	0.862807	0.332578	30.198486	2.249843	False	log

	rank	elpd_loo	p_loo	elpd_diff	weight	se	dse	warning	scale
Hilbert	0	319.036085	39.267076	0.000000	0.667422	30.049570	0.000000	False	log
Latent	1	318.173279	35.403501	0.862807	0.332578	30.198486	2.249843	False	log