"""Implied Gaussian width of the IRI objective probabilities vs the realised error of the plume mean, 2011-12..2025-12.
sigma from each tail with 3-97% probability: (+-thr - mu) / Phi^-1(p), thr = 0.45 C before 2013-05 and 0.5 C after (IRI's thresholds), mu = parsed plume 'COMBINED AVG' (approximate:
the two tails disagree by ~0.1 C median, so mu or the model set differ slightly from IRI's). Reads data/obj_pairs_v5.csv
and data/plume.csv. Writes results/sigma_vs_rmse_v5.csv."""
import pandas as pd, numpy as np, os
from scipy.stats import norm
R = os.environ.get('ROOT', '/work/enso-objective')
q = pd.read_csv(f'{R}/data/obj_pairs_v5.csv'); pl = pd.read_csv(f'{R}/data/plume.csv')
w = pl[pl.model == 'COMBINED AVG'].set_index(['issue', 'season']).nino34.rename('mu')
q = q.join(w, on=['issue', 'season']); q = q[(q.issue >= '2011-12') & (q.issue <= '2025-12')].dropna(subset=['mu']).copy()
thr = np.where(q.issue < '2013-05', 0.45, 0.5)
s1 = np.where(q.q_lan.between(.03, .97), (-thr - q.mu) / norm.ppf(q.q_lan.clip(1e-3, 1 - 1e-3)), np.nan)
s2 = np.where(q.q_eln.between(.03, .97), (thr - q.mu) / norm.ppf(1 - q.q_eln.clip(1e-3, 1 - 1e-3)), np.nan)
S = np.column_stack([s1, s2]); q['sig'] = np.where(np.isnan(S).all(axis=1), np.nan, np.nanmean(np.where(np.isnan(S), np.nan, S), axis=1))
q.loc[(q.sig <= 0) | (q.sig > 3), 'sig'] = np.nan; q['e'] = q.mu - q.oni
g = q.groupby('lead').agg(n=('e', 'size'), bias=('e', 'mean'), rmse=('e', lambda x: np.sqrt((x ** 2).mean())), n_sig=('sig', 'count'), sigma_median=('sig', 'median'))
g['ratio'] = g.sigma_median / g.rmse; g.to_csv(f'{R}/results/sigma_vs_rmse_v5.csv'); print(g.round(3).to_string())
print('tail disagreement median |s1-s2|:', round(float(np.nanmedian(np.abs(s1 - s2))), 3))