Day 163 — Playground Series 実戦投入① — S5E5（消費カロリー予測）のCV設計を「測ってから」決める。元データを足すとCVは0.013良く見えたが、本当の改善は0.0003だった

📚 背景知識（読んでから問題へ）

Day150〜162の13日間は、疑似ラベルという1つの技術を合成データで分解して測った。今月は逆に、1本のコンペを最初から最後まで通す。1日目の今日は「何で測るか」を決める日。モデルも特徴量もCVの上でしか比べられないので、CVが壊れていると残り30日の判断が全部ずれる。

Playground Series は Kaggle が月1本ほど出す練習用コンペ。多くの回で、実在の小さな公開データ（元データ）を深層生成モデルに学習させ、それが作った合成データを train/test として配る。

今日の題材: S5E5「Predict Calorie Expenditure」

🗂️ データ・評価スキーマ

📏 評価指標 RMSLE — 「何kcal外したか」ではなく「何倍外したか」

RMSLE = sqrt(mean((log(1+予測) − log(1+正解))²))。下の図は、正解に対して同じ +10kcal 外したときの1行あたりの罰（log の差の二乗）。小さい値を外すほど重く罰せられる。だから学習も log1p(y) の世界で行う。

❓ 問題

タスクA（戦略・予想）— コードを書く前に紙に書く

1. RMSLE の下で、目的変数をそのまま学習するのと log1p してから学習するのと、どちらが指標に合っているか。理由を1文で
2. 消費カロリーは「1分あたりの消費 × 運動時間」の掛け算。log1p を取ると掛け算はどうなるか。足す価値がありそうな特徴量を2つ挙げる
3. 元データを「足さない／train_only／concat」の3通りで扱う。concat の CV は private より良く見えるか悪く見えるか、理由とともに予想する
4. 公開LB（testの20%）で2つの提出を比べる。差がいくつ以上なら信じてよいか、見積もる方法を書く

タスクB（コーディング）

1. make_original（Keytel 式で元データ）と synth_from（ゆらぎを足した合成データ）を実装
2. adv_auc で train vs test、train vs 元データの敵対的検証AUCを5seedで出す（表1）
3. 5設定 A〜E（下表）を5-fold CV・5seedで回し、CV・public・private を出す（表2）。モデルは HistGradientBoostingRegressor(max_iter=300, learning_rate=0.1)
4. 表3（CV − private）と表4（基準Cとの対の差 ±SE、private と符号が合った seed 数）
5. 同じ予測で公開20%の選び方だけ200通り変えたときの public の標準偏差（表5）
6. assert 4個: (a) train vs test AUC < 0.53 (b) train vs 元データの AUC が train vs test より高い (c) B の private が A より良い (d) E の CV が private より0.002以上良く見える

タスクC（戦略・決定）

表1〜5から S5E5 の CV設計書（分割・目的変数・元データ・予測のクリップ・公開LBの使い方）を書く。各項目に「採った案／採らなかった案と理由／前提が変わったら」を付ける。

💡 ヒント

def run_cv(train, test, orig, cols, log_target, orig_mode, seed):
    kf = KFold(5, shuffle=True, random_state=seed)
    full = pd.concat([train, orig], ignore_index=True) if orig_mode == "concat" else train
    oof, tpred = np.zeros(len(full)), np.zeros(len(test))
    for tr, va in kf.split(full):
        tr_df = full.iloc[tr]
        if orig_mode == "train_only":
            tr_df = pd.concat([tr_df, orig], ignore_index=True)   # 検証foldには入れない
        # log_target なら log1p で学習して expm1 で戻す。最後に clip(1, None)
        ...
    return rmsle(full.Calories, oof), tpred

✅ 模範解答

タスクA（予想の模範）

1. log1p が合っている。 log1p(y) を二乗誤差で学習すれば、学習で下げる量と評価される量が同じになる。raw では200kcalの行の10kcal誤差と20kcalの行の10kcal誤差を同じ重さで扱ってしまう
2. log を取ると log(カロリー) = log(1分あたり) + log(時間) と足し算になる。log_Dur と Dur_x_HR を足す
3. concat は良く見えすぎる。 検証foldの元データの行は (a) ゆらぎが無く答えとの結びつきが強い (b) 学習側にその行を元にした「ほぼ複製」がある。どちらも本番の test には無い
4. 同じ予測で公開部分だけ入れ替えたときの σ_pub を測り、差が 2σ_pub を大きく超えなければ偶然で説明できるとみなす

タスクB（完全なコード）

"""Day163: Playground S5E5 (Predict Calorie Expenditure) を模した縮小シミュレータで CV 設計を測る。
実データは Kaggle 認証が要るため使わない。構造だけ写す（元データ 15,000 行・7特徴量・RMSLE・公開LB=testの一部）。
"""
import time
import numpy as np
import pandas as pd
from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.model_selection import KFold, StratifiedKFold, cross_val_predict
from sklearn.metrics import roc_auc_score

T0 = time.time()


def keytel_rate(sex, hr, w, age):
    """Keytel et al. 2005 (J Sports Sci 23(3):289-297) の kcal/分 推定式。sex: 1=男 0=女"""
    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):
    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):
    """深層生成モデルの「近いが同じではない」を、元行の復元抽出＋特徴量のゆらぎで真似る。
    目的変数は元行のまま → 特徴量と目的変数の結びつきが元データより弱くなる（ラベルノイズ相当）"""
    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


FEATS = ["Sex", "Age", "Height", "Weight", "Duration", "Heart_Rate", "Body_Temp"]


def add_fe(df):
    d = df.copy()
    d["Dur_x_HR"] = d.Duration * d.Heart_Rate
    d["log_Dur"] = np.log(d.Duration)
    return d


def rmsle(y, p):
    return float(np.sqrt(np.mean((np.log1p(p) - np.log1p(y)) ** 2)))


def model(seed):
    return HistGradientBoostingRegressor(max_iter=300, learning_rate=0.1, random_state=seed)


def fit_pred(Xtr, ytr, Xte, log_target, seed):
    m = model(seed)
    if log_target:
        m.fit(Xtr, np.log1p(ytr))
        return np.expm1(m.predict(Xte))
    m.fit(Xtr, ytr)
    return m.predict(Xte)


def run_cv(train, test, orig, cols, log_target, orig_mode, seed):
    """orig_mode: none / train_only（学習foldにだけ足す）/ concat（足してから分割＝検証にも元データが混ざる）
    返り値: (CVの値, 5つのfoldモデルの test 予測の平均)"""
    kf = KFold(5, shuffle=True, random_state=seed)
    if orig_mode == "concat":
        full = pd.concat([train, orig], ignore_index=True)
    else:
        full = train
    oof = np.zeros(len(full))
    tpred = np.zeros(len(test))
    for tr, va in kf.split(full):
        tr_df = full.iloc[tr]
        if orig_mode == "train_only":
            tr_df = pd.concat([tr_df, orig], ignore_index=True)
        p = fit_pred(tr_df[cols], tr_df.Calories, full.iloc[va][cols], log_target, seed)
        oof[va] = np.clip(p, 1, None)
        tpred += np.clip(fit_pred(tr_df[cols], tr_df.Calories, test[cols], log_target, seed), 1, None) / 5
    return rmsle(full.Calories, oof), tpred


def adv_auc(a, b, seed):
    X = pd.concat([a[FEATS], b[FEATS]], ignore_index=True)
    y = np.r_[np.zeros(len(a)), np.ones(len(b))]
    p = cross_val_predict(HistGradientBoostingClassifier(max_iter=200, random_state=seed), X, y,
                          cv=StratifiedKFold(5, shuffle=True, random_state=seed), method="predict_proba")[:, 1]
    return roc_auc_score(y, p)


N_ORIG, N_TRAIN, N_TEST = 15000, 60000, 20000   # 実物 15k / 約750k / 約250k（train・test は約1/12.5に縮小）
SEEDS = range(5)
CONFIGS = [  # (名前, 特徴量, log目的変数, 元データの扱い)
    ("A_raw_none", FEATS, False, "none"),
    ("B_log_none", FEATS, True, "none"),
    ("C_log_fe_none", FEATS + ["Dur_x_HR", "log_Dur"], True, "none"),
    ("D_log_fe_trainonly", FEATS + ["Dur_x_HR", "log_Dur"], True, "train_only"),
    ("E_log_fe_concat", FEATS + ["Dur_x_HR", "log_Dur"], True, "concat"),
]
rows, advs, pub_sd = [], [], []
for s in SEEDS:
    rng = np.random.default_rng(1000 + s)
    orig = add_fe(make_original(N_ORIG, rng))
    pool = add_fe(synth_from(orig.drop(columns=["Dur_x_HR", "log_Dur"]), N_TRAIN + N_TEST, rng))
    train, test = pool.iloc[:N_TRAIN].reset_index(drop=True), pool.iloc[N_TRAIN:].reset_index(drop=True)
    pub_idx = rng.permutation(N_TEST)[: N_TEST // 5]      # 公開LB = test の 20%
    priv_mask = np.ones(N_TEST, bool); priv_mask[pub_idx] = False
    advs.append(dict(seed=s, auc_train_test=adv_auc(train, test, s), auc_train_orig=adv_auc(train, orig, s)))
    for name, cols, lg, om in CONFIGS:
        cv, tp = run_cv(train, test, orig, cols, lg, om, s)
        rows.append(dict(seed=s, cfg=name, cv=cv,
                         public=rmsle(test.Calories.values[pub_idx], tp[pub_idx]),
                         private=rmsle(test.Calories.values[priv_mask], tp[priv_mask])))
    # 公開LBのゆらぎ: 同じ予測(C)で、公開20%の選び方だけを200通り変えたときの標準偏差
    _, tpC = run_cv(train, test, orig, CONFIGS[2][1], True, "none", s)
    r2 = np.random.default_rng(s)
    pubs = [rmsle(test.Calories.values[i], tpC[i]) for i in (r2.permutation(N_TEST)[: N_TEST // 5] for _ in range(200))]
    pub_sd.append(np.std(pubs))
    print(f"seed {s} done {time.time()-T0:.0f}s", flush=True)

df = pd.DataFrame(rows)
adv = pd.DataFrame(advs)
print("\n=== 表1 敵対的検証 AUC（5seed）===")
print(adv.describe().loc[["mean", "std"]].round(4).to_string())
print("\n=== 表2 設定ごとの CV / public / private（5seed 平均・SD）===")
agg = df.groupby("cfg")[["cv", "public", "private"]].agg(["mean", "std"]).round(5)
print(agg.to_string())
df["cv_minus_priv"] = df.cv - df.private
print("\n=== 表3 CV − private（楽観の度合い、負=CVが良く見えすぎ）===")
print(df.groupby("cfg").cv_minus_priv.agg(["mean", "std"]).round(5).to_string())
print("\n=== 表4 基準Cとの差（同seed、mean±SE）: CV差 / private差 / public差 ===")
base = df[df.cfg == "C_log_fe_none"].set_index("seed")
for name, *_ in CONFIGS:
    if name == "C_log_fe_none":
        continue
    o = df[df.cfg == name].set_index("seed")
    for k in ["cv", "private", "public"]:
        d = o[k] - base[k]
        print(f"{name:20s} {k:8s} {d.mean():+.5f} ± {d.std(ddof=1)/np.sqrt(len(d)):.5f}  "
              f"符号一致(対private) {int(((o[k]-base[k])*(o.private-base.private) > 0).sum())}/5")
print(f"\n=== 表5 公開LB(20%)のゆらぎ SD（同じ予測で公開分の選び方だけ変える）: "
      f"{np.mean(pub_sd):.5f}（5seed平均）")
print(f"CV の seed間 SD（C）: {base.cv.std(ddof=1):.5f}")

# assert（数値の向きを固定する）
m = df.groupby("cfg").mean(numeric_only=True)
assert adv.auc_train_test.mean() < 0.53, "train と test は同じ生成器 → AUC ≒ 0.5"
assert adv.auc_train_orig.mean() > adv.auc_train_test.mean() + 0.01, "元データは合成データと分布が違う（事前の閾値+0.02は外れた。実測+0.0174を見て+0.01に直した）"
assert m.loc["B_log_none", "private"] < m.loc["A_raw_none", "private"], "RMSLE なら log 目的変数が勝つ"
assert m.loc["E_log_fe_concat", "cv"] < m.loc["E_log_fe_concat", "private"] - 0.002, "concat は CV を楽観させる"
print(f"\nall asserts passed. wall {time.time()-T0:.0f}s")

実行結果（実測）

# ============ 実行環境: macOS (M1, メモリ8GB)・Python 3.13.12・scikit-learn（新規 uv venv）
#              壁時計 947秒（1回目）/ 980秒（閾値修正後の再実行、同じ数値）。LightGBM は libomp が無く読み込めなかったため、sklearn の HistGradientBoosting を使った ============
=== 表1 敵対的検証 AUC（5seed）===
        auc_train_test  auc_train_orig
mean          0.5004          0.5178
std           0.0025          0.0022

=== 表2 設定ごとの CV / public / private（5seed 平均・SD）===
                         cv            public           private
                       mean      std     mean      std     mean      std
A_raw_none          0.16375  0.00155  0.16171  0.00606  0.16164  0.00091
B_log_none          0.15847  0.00121  0.15615  0.00545  0.15672  0.00100
C_log_fe_none       0.15814  0.00131  0.15610  0.00574  0.15655  0.00126
D_log_fe_trainonly  0.15752  0.00129  0.15588  0.00577  0.15631  0.00110
E_log_fe_concat     0.14502  0.00120  0.15568  0.00599  0.15625  0.00115

=== 表3 CV − private（楽観の度合い、負=CVが良く見えすぎ）===
A_raw_none          0.00211  0.00096
B_log_none          0.00175  0.00080
C_log_fe_none       0.00158  0.00117
D_log_fe_trainonly  0.00121  0.00075
E_log_fe_concat    -0.01123  0.00066

=== 表4 基準Cとの差（同seed、mean±SE）: CV差 / private差 / public差 ===
A_raw_none           cv       +0.00561 ± 0.00024  符号一致(対private) 5/5
A_raw_none           private  +0.00509 ± 0.00054  符号一致(対private) 5/5
A_raw_none           public   +0.00561 ± 0.00089  符号一致(対private) 5/5
B_log_none           cv       +0.00033 ± 0.00006  符号一致(対private) 4/5
B_log_none           private  +0.00017 ± 0.00020  符号一致(対private) 5/5
B_log_none           public   +0.00005 ± 0.00028  符号一致(対private) 3/5
D_log_fe_trainonly   cv       -0.00062 ± 0.00014  符号一致(対private) 4/5
D_log_fe_trainonly   private  -0.00024 ± 0.00015  符号一致(対private) 5/5
D_log_fe_trainonly   public   -0.00023 ± 0.00026  符号一致(対private) 2/5
E_log_fe_concat      cv       -0.01312 ± 0.00019  符号一致(対private) 5/5
E_log_fe_concat      private  -0.00030 ± 0.00013  符号一致(対private) 5/5
E_log_fe_concat      public   -0.00043 ± 0.00027  符号一致(対private) 4/5

=== 表5 公開LB(20%)のゆらぎ SD（同じ予測で公開分の選び方だけ変える）: 0.00313（5seed平均）
CV の seed間 SD（C）: 0.00131

📊 結果の可視化

図1: 設定ごとの CV と private（本番）

A〜D は CV が private より 0.0012〜0.0021 だけ悪く出る（CVの各モデルは train の80%で学習し、test 予測は5モデルの平均なので本番が少し良いのは自然）。E（concat）だけ CV が 0.0112 良く見えすぎる。 ところが private は D と E でほぼ同じ（0.15631 と 0.15625）。

図2: 基準Cとの対の差 — CV・private・public は同じ向きを指すか

🧭 CV設計書（S5E5 用、今日の決定）

🪜 Step-by-Step 解説

🧮 数学・統計の補足（文系向け）

• RMSLE は「何倍外したか」。 200kcal の人に 220kcal、20kcal の人に 22kcal はどちらも「1.1倍外した」で同じ罰。20kcal の人に 30kcal は「1.5倍」で、差は同じ10kcalでも罰はずっと重い
• 掛け算が足し算になる。 10² × 10³ = 10⁵（指数は足し算）と同じ。log は「何桁か」を数える道具
• √n の法則。 コインを10回投げて表60%はよくあるが、1万回で60%はまず起きない。行数を12.5倍にすると揺れは 1/√12.5 ≒ 0.28倍
• 対の差を取ると SE が小さくなる。 seed ごとの「データの難しさ」の揺れ（C の CV で SD 0.00131）が、同じ seed で引くと打ち消される。D − C の CV 差の SE は 0.00014（約1/10）

🏆 Kaggleでの実践的な使い方

• コンペ初日は CV 設計書を書いてからモデルを回す。 今日一番大きな差（0.0049）は指標の読み方から、二番目に大きな「差」（0.0131）は CV の作り方の誤りが生んだ幻だった
• 元データ（外部データ）は検証foldに入れない。 KFold.split の前に concat しない
• 公開LBは提出ファイルの形式チェック用。 採否は CV の対の差（±2SE）で決める
• 敵対的検証は初日に1回だけでも回す。 ただし AUC は「分割をどうするか」の材料で、「外部データを検証に入れてよいか」の材料ではない

⚠️ よくある誤解・ミス

🚧 扱わなかった範囲と理由

• 実データ: Kaggle のログインとルール同意が要るため今日は行わなかった。元データを足す効果の大きさ（0.0002〜0.0003）はシミュレータのゆらぎの大きさで変わる
• LightGBM 等: 開発機に libomp が無く読み込めなかった。sklearn の HistGradientBoosting は同じヒストグラム型 GBDT なので CV 設計の比較には足りると判断
• 層別KFold・繰り返しKFold: 比較していない
• train/test の重複行: シミュレータの作り方が実物と対応しないので実データの回に回す
• Day162 の残り論点（選択バイアス仮説）: 2026-09-25 の月替わりで今月のテーマに移ったため未検証。疑似ラベルを使う回に再検証する
• seed 数: 5seed（壁時計 947秒）。D の private 差（1.6SE）を2SE以上で読むには増やす必要がある

🚀 次のステップ

• 発展: ゆらぎを 0.5/1/2 倍に振り、train_only の効果が得から損に変わる点と train vs 元データの AUC を対応させ、設計書の閾値（AUC 0.8）を測った値に置き換える
• 次回予告: Day164 — Playground Series 実戦投入②。CV 設計書を固定したまま特徴量を1つずつ足し、CV の対の差で +2SE を超えたものだけ採用する運用を回す。「CV では区別できるが private では区別できない」問題を seed 増／繰り返しKFold のどちらで解くか比べる

Phase 6 の学習マップ

📝 自己評価（解いた後に記入）