import os
os.environ["KMP_DUPLICATE_LIB_OK"]="TRUE"
import numpy as np, matplotlib.pyplot as plt, warnings
warnings.filterwarnings("ignore")
from sklearn.model_selection import cross_val_score
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
import statsmodels.api as sm
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; GREY="#a0aec0"; PURP="#6b46c1"
rng=np.random.default_rng(0); n=150
x=np.sort(rng.uniform(-3,3,n)); y=0.5*x**3-2*x+rng.normal(0,3,n)
from sklearn.model_selection import KFold
degs=range(1,11); r2=[]; cvmse=[]; cvblock=[]; aic=[]; bic=[]
for d in degs:
    X=PolynomialFeatures(d).fit_transform(x[:,None]); m=sm.OLS(y,X).fit()
    r2.append(m.rsquared); aic.append(m.aic); bic.append(m.bic)
    sc=lambda cvs: -cross_val_score(LinearRegression(),X,y,cv=cvs,scoring="neg_mean_squared_error").mean()
    cvmse.append(sc(KFold(5,shuffle=True,random_state=0)))    # random folds: predicting a new draw
    cvblock.append(sc(KFold(5)))                              # contiguous folds on sorted x: extrapolation
fig,ax=plt.subplots(1,2,figsize=(13,4.2))
xx=np.linspace(-3,3,200)
ax[0].scatter(x,y,s=10,color=GREY,alpha=.6)
for d,c in [(1,ORANGE),(3,GREEN),(10,RED)]:
    b=np.polyfit(x,y,d); ax[0].plot(xx,np.polyval(b,xx),color=c,lw=2,label=f"degree {d}")
ax[0].plot(xx,0.5*xx**3-2*xx,"k--",lw=1,label="truth (deg 3)"); ax[0].set_title("Polynomial fits"); ax[0].legend(fontsize=8); ax[0].set_ylim(-20,20)
ax[1].plot(list(degs),r2,"o-",color=GREY,lw=2,label="in-sample R² (rises -> trap)")
a2=ax[1].twinx()
a2.plot(list(degs),cvmse,"s-",color=BLUE,lw=2,label="5-fold CV MSE (shuffled folds)")
a2.plot(list(degs),cvblock,"^--",color=ORANGE,lw=1.6,label="contiguous folds on sorted x")
a2.set_yscale("log"); a2.set_ylabel("CV MSE (log scale)",color=BLUE); a2.legend(fontsize=7,loc="upper left")
ax[1].axvline(3,color=GREEN,ls="--"); ax[1].set_xlabel("polynomial degree"); ax[1].set_ylabel("in-sample R²",color=GREY)
ax[1].set_title(f"CV picks degree {list(degs)[int(np.argmin(cvmse))]} (truth=3); R² just keeps rising")
plt.tight_layout(); plt.show()
print(f"in-sample R² rises monotonically to {r2[-1]:.2f} -- the trap. Cross-validation is U-shaped and bottoms out at")
print(f"degree {list(degs)[int(np.argmin(cvmse))]} with shuffled folds, and at degree {list(degs)[int(np.argmin(cvblock))]} with contiguous ones. Same answer here, but not the same measurement:")
print(f"   degree  {'shuffled':>12} {'contiguous':>14}")
for i,dd in enumerate(degs):
    print(f"   {dd:>6}  {cvmse[i]:>12.1f} {cvblock[i]:>14,.0f}")
print(f"The contiguous version has to extrapolate, and a degree-10 polynomial outside its training range is a disaster:")
print(f"its CV MSE reaches {max(cvblock):,.0f} against {max(cvmse):.0f} for shuffled folds. That is why the axis is logarithmic --")
print("on a linear scale the curve leaves the plot entirely and the 'U' is really an explosion. Both curves select the")
print("true degree here, but only the shuffled one is estimating the quantity the section claims to estimate.")

in-sample R² rises monotonically to 0.43 -- the trap. Cross-validation is U-shaped and bottoms out at
degree 3 with shuffled folds, and at degree 3 with contiguous ones. Same answer here, but not the same measurement:
   degree      shuffled     contiguous
        1          14.8             18
        2          15.4             50
        3          10.8             11
        4          11.0             12
        5          11.5             39
        6          11.6          1,146
        7          11.7          2,412
        8          11.6         52,821
        9          12.0          6,102
       10          12.0      8,488,046
The contiguous version has to extrapolate, and a degree-10 polynomial outside its training range is a disaster:
its CV MSE reaches 8,488,046 against 15 for shuffled folds. That is why the axis is logarithmic --
on a linear scale the curve leaves the plot entirely and the 'U' is really an explosion. Both curves select the
true degree here, but only the shuffled one is estimating the quantity the section claims to estimate.

fig,ax=plt.subplots(figsize=(7.5,4.2))
ax.plot(list(degs),aic,"o-",color=BLUE,lw=2,label=f"AIC (min @ deg {list(degs)[int(np.argmin(aic))]})")
ax.plot(list(degs),bic,"s-",color=RED,lw=2,label=f"BIC (min @ deg {list(degs)[int(np.argmin(bic))]})")
ax.axvline(3,color=GREEN,ls="--",label="truth = 3"); ax.set_xlabel("polynomial degree"); ax.set_ylabel("information criterion (lower=better)"); ax.set_title("AIC & BIC both select the true complexity"); ax.legend(fontsize=8)
plt.tight_layout(); plt.show()
print(f"AIC selects degree {list(degs)[int(np.argmin(aic))]}, BIC degree {list(degs)[int(np.argmin(bic))]} (BIC's heavier ln(n) penalty makes it the more parsimonious).")
print("Both agree with cross-validation -- and need only one fit each, no held-out folds. The catch: they require a likelihood.")

AIC selects degree 3, BIC degree 3 (BIC's heavier ln(n) penalty makes it the more parsimonious).
Both agree with cross-validation -- and need only one fit each, no held-out folds. The catch: they require a likelihood.

import pymc as pm, arviz as az
bdeg=range(1,9); elpd=[]; se=[]; idatas={}
for d in bdeg:
    X=PolynomialFeatures(d).fit_transform(x[:,None]); Xs=X.copy(); Xs[:,1:]=(X[:,1:]-X[:,1:].mean(0))/X[:,1:].std(0)
    with pm.Model():
        b=pm.Normal("b",0,10,shape=X.shape[1]); s=pm.HalfNormal("s",5); pm.Normal("y",Xs@b,s,observed=y)
        idata=pm.sample(600,tune=600,chains=2,cores=1,progressbar=False,random_seed=0,idata_kwargs={"log_likelihood":True})
    lo=az.loo(idata); elpd.append(float(lo.elpd)); se.append(float(lo.se)); idatas[f"deg {d}"]=idata
elpd=np.array(elpd); se=np.array(se); best=list(bdeg)[int(np.argmax(elpd))]
cmp=az.compare(idatas)                                          # paired comparison: elpd_diff and its own SE
fig,ax=plt.subplots(1,2,figsize=(13,4.2))
ax[0].errorbar(list(bdeg),elpd,yerr=se,fmt="o-",color=PURP,lw=2,capsize=3)
ax[0].axvline(3,color=GREEN,ls="--",label="truth = 3"); ax[0].legend(fontsize=8)
ax[0].set_xlabel("polynomial degree"); ax[0].set_ylabel("LOO ELPD"); ax[0].set_title("Error bars = SE(ELPD): everything overlaps")
_o=[int(k.split()[1]) for k in cmp.index]; _d=cmp["elpd_diff"].values; _ds=cmp["dse"].values
ax[1].errorbar(_o,_d,yerr=_ds,fmt="o",color=PURP,capsize=3); ax[1].axhline(0,color=GREY,lw=.8)
ax[1].axvline(3,color=GREEN,ls="--"); ax[1].set_xlabel("polynomial degree")
ax[1].set_ylabel("ELPD difference from the best"); ax[1].set_title("Error bars = SE(difference): the honest comparison")
plt.tight_layout(); plt.show()

print("Sampling and PSIS diagnostics -- an ELPD from a bad posterior is not worth comparing:")
print(f"   {'deg':>4} {'ELPD':>9} {'SE':>7} {'max Pareto-k':>14} {'k>0.7':>7} {'max R-hat':>11} {'divergences':>12}")
for d,idata,e,s_ in zip(bdeg,idatas.values(),elpd,se):
    _lo=az.loo(idata,pointwise=True); _k=np.asarray(_lo.pareto_k); _su=az.summary(idata,var_names=["b","s"])
    print(f"   {d:>4} {e:>9.1f} {s_:>7.1f} {float(_k.max()):>14.2f} {int((_k>0.7).sum()):>7d} "
          f"{float(_su['r_hat'].max()):>11.3f} {int(idata.sample_stats['diverging'].sum()):>12d}")
print("   All Pareto-k well below 0.7, no divergences, R-hat at 1.00-1.01: the LOO estimates are trustworthy.")

print(f"\nLOO ELPD is maximized at degree {best} (truth=3), agreeing with CV, AIC and BIC on this sample.")
print("But read the two panels together. The left one, plotting SE(ELPD), suggests nothing can be distinguished at all.")
print("The right one uses the paired standard error and tells a sharper and more useful story:")
print(cmp[["rank","elpd_diff","dse","weight"]].to_string())
print("\nDegrees 1 and 2 are decisively rejected -- gaps of 20 and 30 ELPD against a paired SE of about 6.5. Among degrees")
print("3 through 8 nothing is separable: the gaps are 0 to 2 with paired SEs of 0.7 to 1.2, and arviz flags them all as")
print("|elpd_diff| < 4. LOO is not really saying 'the answer is 3'. It is saying 'at least 3, and beyond that this sample")
print("cannot tell you' -- which is the correct answer for a criterion that targets PREDICTION, because a degree-5 fit")
print("predicts a cubic almost exactly as well as a cubic does.")
print("\nLOO IS cross-validation done analytically over the posterior: az.loo estimates leave-one-out predictive accuracy")
print("directly, which is why it behaves like AIC and unlike BIC -- the subject of the next section.")

g++ not available, if using conda: `conda install gxx`

Initializing NUTS using jitter+adapt_diag...

Sequential sampling (2 chains in 1 job)

NUTS: [b, s]

Sampling 2 chains for 600 tune and 600 draw iterations (1_200 + 1_200 draws total) took 0 seconds.

We recommend running at least 4 chains for robust computation of convergence diagnostics

Initializing NUTS using jitter+adapt_diag...

Sequential sampling (2 chains in 1 job)

NUTS: [b, s]

Sampling 2 chains for 600 tune and 600 draw iterations (1_200 + 1_200 draws total) took 0 seconds.

from sklearn.model_selection import KFold as _KF
def draw(seed,n=150):
    r=np.random.default_rng(seed); xx=np.sort(r.uniform(-3,3,n))
    return xx, 0.5*xx**3-2*xx+r.normal(0,3,n)
def pick(xx,yy,seed=0,dmax=10):
    _r2,_cv,_a,_b=[],[],[],[]
    for d in range(1,dmax+1):
        Xd=PolynomialFeatures(d).fit_transform(xx[:,None]); mm=sm.OLS(yy,Xd).fit()
        _a.append(mm.aic); _b.append(mm.bic)
        _cv.append(-cross_val_score(LinearRegression(),Xd,yy,cv=_KF(5,shuffle=True,random_state=seed),
                                    scoring="neg_mean_squared_error").mean())
    return int(np.argmin(_cv))+1, int(np.argmin(_a))+1, int(np.argmin(_b))+1

NS=200; sel=np.array([pick(*draw(s),seed=s) for s in range(NS)])
names=["cross-validation","AIC","BIC"]
print(f"{NS} fresh datasets from the same cubic process, n=150:")
print(f"   {'criterion':18s} {'P(picks 3)':>11} {'P(too complex)':>15} {'P(too simple)':>14}")
for j,nm in enumerate(names):
    v=sel[:,j]; print(f"   {nm:18s} {np.mean(v==3):>11.2f} {np.mean(v>3):>15.2f} {np.mean(v<3):>14.2f}")
print(f"   all three agree on 3 in the SAME draw: {np.mean((sel==3).all(1)):.2f}")
print(f"\nSeed 0 -- the draw shown in sections 1-3 -- is one of the {np.mean((sel==3).all(1)):.0%} where they agree. It is a representative")
print("illustration, not a demonstration that the criteria are interchangeable.")

# the Bayesian criterion, on fewer draws because each one costs a full posterior
LS=10; bpick=[]
for s in range(LS):
    xs,ys_=draw(s); ee=[]
    for d in range(1,9):
        Xd=PolynomialFeatures(d).fit_transform(xs[:,None]); Xsd=Xd.copy()
        Xsd[:,1:]=(Xd[:,1:]-Xd[:,1:].mean(0))/Xd[:,1:].std(0)
        with pm.Model():
            bb=pm.Normal("b",0,10,shape=Xd.shape[1]); ss=pm.HalfNormal("s",5); pm.Normal("y",Xsd@bb,ss,observed=ys_)
            _id=pm.sample(600,tune=600,chains=2,cores=1,progressbar=False,random_seed=s,
                          idata_kwargs={"log_likelihood":True})
        ee.append(float(az.loo(_id).elpd))
    bpick.append(int(np.argmax(ee))+1)
print(f"\nPSIS-LOO over {LS} draws: picks {bpick}  ->  P(picks 3) = {np.mean(np.array(bpick)==3):.2f}")
print("LOO tracks AIC and CV, as theory says it should: all three are estimates of out-of-sample predictive accuracy.")

fig,ax=plt.subplots(1,2,figsize=(13,4.2))
w=0.27
for j,(nm,c) in enumerate(zip(names,[BLUE,ORANGE,RED])):
    vals,cnts=np.unique(sel[:,j],return_counts=True)
    ax[0].bar(vals+(j-1)*w,cnts/NS,width=w,color=c,label=nm)
ax[0].axvline(3,color=GREEN,ls="--"); ax[0].set_xlabel("degree selected"); ax[0].set_ylabel("frequency")
ax[0].set_title(f"Which degree does each criterion choose? ({NS} draws)"); ax[0].legend(fontsize=8)

NN=[150,400,1000,2500]; res={nm:[] for nm in names}
for n_ in NN:
    ss_=np.array([pick(*draw(s,n_),seed=s) for s in range(80)])
    for j,nm in enumerate(names): res[nm].append(np.mean(ss_[:,j]==3))
for (nm,c) in zip(names,[BLUE,ORANGE,RED]):
    ax[1].plot(NN,res[nm],"o-",color=c,lw=2,label=nm)
ax[1].set_xscale("log"); ax[1].set_xlabel("sample size n"); ax[1].set_ylabel("P(selects the true degree)")
ax[1].set_ylim(0,1.05); ax[1].set_title("More data does not rescue AIC or CV"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.show()
print(f"\n   {'n':>6} " + " ".join(f"{nm:>18}" for nm in names))
for i,n_ in enumerate(NN):
    print(f"   {n_:>6} " + " ".join(f"{res[nm][i]:>18.2f}" for nm in names))
print("\nBIC is already at the ceiling by n=150 and stays there. AIC and cross-validation sit near 0.7 at EVERY sample")
print("size tried, including n=2500 -- more data does not help them. That is what inconsistency means: the probability")
print("of over-selecting does not go to zero as the sample grows, because the penalty per parameter does not grow with n.")
print("\nThis is not AIC failing. It is AIC succeeding at a different task: among models that predict about equally well,")
print("nothing in a predictive criterion prefers the smaller one, and the extra parameters cost almost no accuracy. The")
print("cost of AIC's mistake is small in the currency AIC cares about, and total in the currency BIC cares about.")
print("Choose the criterion by the question -- 'which model will predict best?' or 'which model generated this?' -- and")
print("stop expecting one number to answer both.")

Initializing NUTS using jitter+adapt_diag...

200 fresh datasets from the same cubic process, n=150:
   criterion           P(picks 3)  P(too complex)  P(too simple)
   cross-validation          0.68            0.33           0.00
   AIC                       0.72            0.28           0.00
   BIC                       0.97            0.03           0.00
   all three agree on 3 in the SAME draw: 0.60

Seed 0 -- the draw shown in sections 1-3 -- is one of the 60% where they agree. It is a representative
illustration, not a demonstration that the criteria are interchangeable.

Sequential sampling (2 chains in 1 job)

NUTS: [b, s]

Sampling 2 chains for 600 tune and 600 draw iterations (1_200 + 1_200 draws total) took 0 seconds.

We recommend running at least 4 chains for robust computation of convergence diagnostics

Initializing NUTS using jitter+adapt_diag...

Sequential sampling (2 chains in 1 job)

NUTS: [b, s]

Sampling 2 chains for 600 tune and 600 draw iterations (1_200 + 1_200 draws total) took 0 seconds.

criterion	how	cost	needs
cross-validation	refit on held-out folds	high (K refits)	nothing (model-agnostic)
AIC	$-2\ell+2k$	one fit	a likelihood
BIC	$-2\ell+k\ln n$ (parsimonious)	one fit	a likelihood
PSIS-LOO / WAIC	posterior predictive, LOO via importance sampling	one posterior	a Bayesian fit

Model Evaluation I — Model selection: CV vs AIC/BIC vs Bayesian LOO¶

Four routes to "which model?", and why they agree (Akaike, Schwarz; Watanabe; Vehtari-Gelman-Gabry)¶

1. The overfitting trap, and cross-validation¶

2. Information criteria — AIC and BIC¶

3. Bayesian model comparison — PSIS-LOO¶

4. Does the agreement survive a new sample?¶

5. Summary — four routes, two questions¶