Day 165 — Playground Series 実戦投入③ — Optuna でハイパーパラメータを60回探す。探索に使った CV は改善を最大1.7倍に見せ、探索回数を増やすと本当は悪い設定を選んだワールドもあった。探索に使っていない分割で測り直すと、採否6件がすべて真の効果と一致した

"""Day165: Playground S5E5 縮小シミュレータで「Optuna の探索回数を増やすほど CV が楽観する（探索の過学習）」を測る。
探索は分割seed 0 の5分割（探索用CV）だけを見て行う。選ばれた設定を
 (a) 探索に使っていない分割seed 1〜3 の R2（Day164 の採否ルール）と
 (b) test 6万行（真の効果）
で測り直し、既定値との差が探索用CVでどれだけ大きく見えていたかを数える。
使い方: OMP_NUM_THREADS=2 python day165.py <world>
"""
import sys, time, json
import numpy as np, pandas as pd, optuna
from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.model_selection import KFold

T0 = time.time()
W = int(sys.argv[1])
N_ORIG, N_TRAIN, N_TEST, N_TRIALS = 15000, 30000, 60000, 60
CHECK = [1, 5, 10, 20, 40, 60]


def keytel_rate(sex, hr, w, age):
    m = (-55.0969 + 0.6309 * hr + 0.1988 * w + 0.2017 * age) / 4.184
    f = (-20.4022 + 0.4472 * hr - 0.1263 * w + 0.074 * age) / 4.184
    return np.where(sex == 1, m, f)


def make_original(n, rng):  # Day163/164 と同じ
    sex = rng.integers(0, 2, n); age = rng.integers(20, 80, n)
    h = np.where(sex == 1, rng.normal(175, 7, n), rng.normal(162, 6.5, n)).round()
    w = (rng.normal(24.5, 2.5, n) * (h / 100) ** 2).round()
    dur = rng.integers(1, 31, n).astype(float)
    hr = np.clip(80 + 1.2 * dur + rng.normal(0, 5, n), 67, 128).round()
    temp = (37 + 0.35 * dur ** 0.5 + rng.normal(0, 0.2, n)).round(1)
    rate = np.clip(keytel_rate(sex, hr, w, age), 0.3, None)
    cal = np.clip((rate * dur * np.exp(rng.normal(0, 0.05, n))).round(), 1, None)
    return pd.DataFrame(dict(Sex=sex, Age=age, Height=h, Weight=w, Duration=dur,
                             Heart_Rate=hr, Body_Temp=temp, Calories=cal))


def synth_from(orig, n, rng):  # Day163/164 と同じ
    s = orig.sample(n, replace=True, random_state=int(rng.integers(1 << 30))).reset_index(drop=True)
    s["Age"] = np.clip((s.Age + rng.normal(0, 1.5, n)).round(), 20, 79)
    s["Height"] = (s.Height + rng.normal(0, 2, n)).round()
    s["Weight"] = (s.Weight + rng.normal(0, 2, n)).round()
    s["Duration"] = np.clip((s.Duration + rng.normal(0, 1.0, n)).round(), 1, 30)
    s["Heart_Rate"] = (s.Heart_Rate + rng.normal(0, 2, n)).round()
    s["Body_Temp"] = (s.Body_Temp + rng.normal(0, 0.1, n)).round(1)
    s["Calories"] = np.clip((s.Calories * np.exp(rng.normal(0, 0.03, n))).round(), 1, None)
    return s


def add_feats(d):  # 決定6（Day164）で採用した2列だけ
    d = d.copy()
    d["Dur_x_HR"] = d.Duration * d.Heart_Rate
    d["Keytel_kcal"] = np.clip(keytel_rate(d.Sex.values, d.Heart_Rate.values, d.Weight.values, d.Age.values),
                               0.3, None) * d.Duration
    return d


COLS = ["Sex", "Age", "Height", "Weight", "Duration", "Heart_Rate", "Body_Temp", "Dur_x_HR", "Keytel_kcal"]
DEFAULT = dict(learning_rate=0.1, max_iter=300, max_leaf_nodes=31, min_samples_leaf=20, l2_regularization=0.0)


def sle(y, p):
    return (np.log1p(p) - np.log1p(y)) ** 2


def run(params, split_seed):
    """CV 設計書どおり。early_stopping=False にして、モデル側の乱数を無くす（同じ設定・同じ分割なら同じ予測）"""
    oof, tp, fs = np.zeros(len(train)), np.zeros(len(test)), []
    for tr, va in KFold(5, shuffle=True, random_state=split_seed).split(train):
        tr_df = pd.concat([train.iloc[tr], orig], ignore_index=True)
        m = HistGradientBoostingRegressor(early_stopping=False, random_state=0, **params)
        m.fit(tr_df[COLS], np.log1p(tr_df.Calories))
        pv = np.clip(np.expm1(m.predict(train.iloc[va][COLS])), 1, None)
        tp += np.clip(np.expm1(m.predict(test[COLS])), 1, None) / 5
        oof[va] = pv
        fs.append(np.sqrt(sle(ytr[va], pv).mean()))
    return np.array(fs), np.sqrt(sle(ytr, oof).mean()), tp


rng = np.random.default_rng(2000 + W)        # Day164 のワールド W と同じデータ
orig = make_original(N_ORIG, rng)
pool = synth_from(orig, N_TRAIN + N_TEST, rng)
orig, pool = add_feats(orig), add_feats(pool)
train, test = pool.iloc[:N_TRAIN].reset_index(drop=True), pool.iloc[N_TRAIN:].reset_index(drop=True)
ytr, yt = train.Calories.values, test.Calories.values

trials, preds = [], []


def objective(t):
    p = dict(learning_rate=t.suggest_float("learning_rate", 0.02, 0.3, log=True),
             max_iter=t.suggest_int("max_iter", 100, 600, step=50),
             max_leaf_nodes=t.suggest_int("max_leaf_nodes", 7, 63, log=True),
             min_samples_leaf=t.suggest_int("min_samples_leaf", 5, 200, log=True),
             l2_regularization=t.suggest_float("l2_regularization", 1e-4, 10, log=True))
    fs, cv, tp = run(p, 0)                   # 探索が見るのは分割seed 0 の CV だけ
    trials.append(dict(trial=t.number, cv_search=cv, test=np.sqrt(sle(yt, tp).mean()), **p,
                       folds=json.dumps(list(fs))))
    preds.append(tp)
    return cv


optuna.logging.set_verbosity(optuna.logging.WARNING)
study = optuna.create_study(direction="minimize", sampler=optuna.samplers.TPESampler(seed=W))
d0 = dict(DEFAULT); d0["l2_regularization"] = 1e-4   # 探索空間の下限（0 は log 空間に置けない）
study.enqueue_trial(d0)                      # trial 0 = 既定値（探索0回の状態）
study.optimize(objective, n_trials=N_TRIALS)
tr_df = pd.DataFrame(trials)
print(f"world {W} search done {time.time()-T0:.0f}s", flush=True)

# 探索回数 N ごとの「その時点のベスト」
best_at = {n: int(tr_df.cv_search.iloc[:n].idxmin()) for n in CHECK}
fresh = {}
for i in sorted(set(best_at.values())):
    p = {k: tr_df.loc[i, k] for k in DEFAULT}
    p["max_iter"], p["max_leaf_nodes"], p["min_samples_leaf"] = int(p["max_iter"]), int(p["max_leaf_nodes"]), int(p["min_samples_leaf"])
    fresh[i] = [run(p, s) for s in (1, 2, 3)]   # 探索に使っていない分割 seed
print(f"world {W} re-eval done {time.time()-T0:.0f}s", flush=True)

rows = []
f0 = fresh[0]
for n in CHECK:
    i = best_at[n]
    fs0 = np.array(json.loads(tr_df.folds[0])); fsi = np.array(json.loads(tr_df.folds[i]))
    dS = tr_df.cv_search[i] - tr_df.cv_search[0]
    seS = (fsi - fs0).std(ddof=1) / np.sqrt(5) if i else 0.0
    d15 = np.concatenate([fresh[i][r][0] - f0[r][0] for r in range(3)])
    dF, seF = d15.mean(), (d15.std(ddof=1) / np.sqrt(15) if i else 0.0)
    pi = np.mean([fresh[i][r][2] for r in range(3)], axis=0)   # 真の効果: 3繰返し×5モデル平均（Day164 と同じ）
    pb = np.mean([f0[r][2] for r in range(3)], axis=0)
    ei, eb = sle(yt, pi), sle(yt, pb)
    dT = np.sqrt(ei.mean()) - np.sqrt(eb.mean())
    seT = (ei - eb).std(ddof=1) / np.sqrt(len(yt)) / (2 * np.sqrt(eb.mean())) if i else 0.0
    rows.append(dict(world=W, N=n, best_trial=i, dS=dS, seS=seS, dF=dF, seF=seF, dT=dT, seT=seT,
                     **{k: tr_df.loc[i, k] for k in DEFAULT}))
out = pd.DataFrame(rows)
tr_df.to_csv(f"trials_w{W}.csv", index=False)
out.to_csv(f"check_w{W}.csv", index=False)
pd.set_option("display.width", 250)
print(out.round(5).to_string())
print(f"world {W} wall {time.time()-T0:.0f}s")

import glob, numpy as np, pandas as pd
pd.set_option("display.width", 250)
ck = pd.concat([pd.read_csv(f) for f in sorted(glob.glob("check_w*.csv"))])
tr = {int(f[-5]): pd.read_csv(f) for f in sorted(glob.glob("trials_w*.csv"))}
ck["opt_S_T"] = ck.dS - ck.dT          # 探索用CVが見せた改善 − 真の改善（負 = 楽観の逆）
ck["opt_S_F"] = ck.dS - ck.dF
print("=== 表1 探索回数Nごとのベスト設定: 既定値との差（×1e-4）===")
t = ck.copy()
for c in ["dS", "seS", "dF", "seF", "dT", "seT", "opt_S_T", "opt_S_F"]:
    t[c] = (t[c] * 1e4).round(2)
print(t[["world", "N", "best_trial", "dS", "seS", "dF", "seF", "dT", "seT", "opt_S_T", "opt_S_F",
         "learning_rate", "max_iter", "max_leaf_nodes", "min_samples_leaf", "l2_regularization"]].to_string(index=False))
print("\n=== 表2 全60試行: 探索用CV と test の関係 ===")
for w, d in tr.items():
    rho = d.cv_search.rank().corr(d.test.rank())
    top = d.nsmallest(10, "cv_search")
    rho10 = top.cv_search.rank().corr(top.test.rank())
    best_cv, best_test = d.cv_search.idxmin(), d.test.idxmin()
    print(f"world {w}: Spearman 全60 {rho:.3f} | 上位10 {rho10:.3f} | CV最良 trial {best_cv} (test順位 {int(d.test.rank()[best_cv])}/60) | "
          f"test最良 trial {best_test} (CV順位 {int(d.cv_search.rank()[best_test])}/60) | "
          f"CV幅 {1e4*(d.cv_search.max()-d.cv_search.min()):.1f}e-4 test幅 {1e4*(d.test.max()-d.test.min()):.1f}e-4 | "
          f"CV最良から+0.0002以内 {(d.cv_search <= d.cv_search.min()+2e-4).sum()}試行")

print("\n=== 表3 採否判定（既定値以外が選ばれた点のみ）: −2SE ルール 対 真の効果 ===")
nd = ck[ck.best_trial != 0].copy()
nd["truth"] = nd.dT < -2 * nd.seT
nd["adopt_search"] = nd.dS < -2 * nd.seS
nd["adopt_fresh"] = nd.dF < -2 * nd.seF
nd["z_S"], nd["z_F"], nd["z_T"] = nd.dS / nd.seS, nd.dF / nd.seF, nd.dT / nd.seT
print(nd[["world", "N", "z_S", "z_F", "z_T", "adopt_search", "adopt_fresh", "truth"]].round(2).to_string(index=False))
for k in ["adopt_search", "adopt_fresh"]:
    print(f"{k}: 一致 {(nd[k] == nd.truth).sum()}/{len(nd)} | 見逃し {(~nd[k] & nd.truth).sum()} | 偽採用 {(nd[k] & ~nd.truth).sum()}")
print("\n=== 表4 楽観（真の改善 − 見えた改善、×1e-4。正 = 見かけの方が大きい）===")
nd["楽観_探索CV"] = ((nd.dT - nd.dS) * 1e4).round(2)
nd["楽観_新分割R2"] = ((nd.dT - nd.dF) * 1e4).round(2)
print(nd[["world", "N", "楽観_探索CV", "楽観_新分割R2"]].to_string(index=False))
a = nd[nd.N == 60]
assert ((a.dT - a.dF).abs() < (a.dT - a.dS).abs()).all(), "N=60 では新しい分割の R2 の方が真の効果に近い"
w1 = nd[nd.world == 1].set_index("N")
assert w1.loc[60, "dT"] - w1.loc[60, "dS"] > w1.loc[10, "dT"] - w1.loc[10, "dS"], "探索回数が増えると楽観が増える（world 1）"
assert (nd.adopt_fresh == nd.truth).all(), "新しい分割の R2 で判定すると6件すべて真の効果と一致"
print("all asserts passed")