Day 164 — Playground Series 実戦投入② — 特徴量を1つずつ足して採否を決める。1回の5分割では本当に効く8件のうち4件を見逃し、3回繰り返すと見逃し1件・偽採用1件になった。昨日採用した log_Dur は差がちょうど0だった

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

昨日の宿題。 Day163 で特徴量 Dur_x_HR・log_Dur を足すと、CV では +5.5SE で区別できたのに private（16,000行）では 0.9SE で区別できなかった。今日はこれを「特徴量をどう採否判定するか」という運用の問題として扱う。

今日の問い: 「CV の対の差が −2SE を下回ったら採用」というルールは、SE の取り方で採用・見逃し・偽採用がどれだけ変わるか。正解（真の効果）は train とは別に作った test 6万行で測る。

🗂️ データ・候補スキーマ

📏 評価の仕組み — 対の差と「−2SE」の判定線

指標は Day163 と同じ RMSLE（何倍外したか）。今日評価するのは採否ルールの方で、判定は「対の差の平均が −2SE を下回るか」。下の図は同じ候補（真の改善 −0.0005 前後）を、SE が違う2つのルールで判定したときの模式図。SE が大きいと、同じ改善でも判定線をまたげず見逃す。

点の位置 −4.8 は Dur_x_HR の R2 平均差、SE は表2の中央値（×10⁻⁴）。

❓ 問題

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

1. 候補7つそれぞれについて「真に効く／効かない／悪化する」を予想する
2. log_Dur = log(Duration) は GBDT に対してどんな効果を持つか。理由とともに
3. 陰性対照（noise・HR_shuffled）を足すと、真の効果は 0 になるか、悪化するか
4. R1（1回の5分割・fold差SE）と R2（3回繰り返し・15fold差の素朴なSE）では、どちらが見逃しが多く、どちらが偽採用が多いか

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

1. Day163 の make_original・synth_from を再利用し、add_candidates(df, rng) で候補7列を足す
2. run(train, test, orig, cols, split_seed) を CV 設計書どおりに実装（ランダム5分割・log1p・元データは学習foldにだけ・下限1クリップ）。1回の fit で検証fold と test を両方予測する（Day163 は2回 fit していた）
3. 4ワールド × 3分割seed × (基準 + 候補7) を回し、R1・R2・R3 の対の差と SE、真の効果とその SE を出す
4. 表1（候補ごとの平均と採用回数）、表2（SE の大きさ）、表3（ルールごとの正しい採用・偽採用・見逃し）、表4（CV差と真の効果の順位相関）
5. assert 2つ: (a) Keytel_kcal は4ワールド中3以上で真に効く (b) 15fold差の SE は1回の5分割の SE より小さい（28件の70%超）

タスクC（分析・決定）

表1〜4から (1) 採否ルールを1つ決め、(2) 今日の候補の採否を決めて CV 設計書に「決定6: 特徴量の採否」を書き足す。「採った案／採らなかった案と理由／前提が変わったら」を付ける。

💡 ヒント

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(max_iter=300, learning_rate=0.1, random_state=split_seed)
    m.fit(tr_df[cols], np.log1p(tr_df.Calories))
    pv, pt = [np.clip(np.expm1(m.predict(X)), 1, None) for X in (train.iloc[va][cols], test[cols])]
    ...
fd = fold_scores_cand - fold_scores_base          # 対の差（5個）
adopt = fd.mean() < -2 * fd.std(ddof=1) / np.sqrt(5)

✅ 模範解答

タスクA（予想の模範と実測の照合）

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

"""Day164: Playground S5E5 縮小シミュレータで「特徴量を1つずつ足し、CVの対の差で採否を決める」運用を測る。
実コンペでは train は1つしかない（Day163 の「seed＝データを作り直す」は使えない）。
使える揺れの源は (1) fold 間の差 (2) 分割 seed を変えた繰り返し KFold の2つ。
採否ルール3種の「偽採用（効かない特徴量を採る）」「見逃し」を、大きな test（真の効果）と突き合わせる。
"""
import time
import numpy as np
import pandas as pd
from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.model_selection import KFold

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):  # Day163 と同じ
    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 と同じ
    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


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


def add_candidates(d, rng):
    d = d.copy()
    d["Dur_x_HR"] = d.Duration * d.Heart_Rate
    d["log_Dur"] = np.log(d.Duration)
    d["BMI"] = d.Weight / (d.Height / 100) ** 2
    d["HR_x_W"] = d.Heart_Rate * d.Weight
    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
    d["noise"] = rng.normal(0, 1, len(d))                  # 陰性対照1: 目的変数と無関係
    d["HR_shuffled"] = rng.permutation(d.Heart_Rate.values)  # 陰性対照2: 分布は本物、中身は無関係
    return d


CANDS = ["Dur_x_HR", "log_Dur", "BMI", "HR_x_W", "Keytel_kcal", "noise", "HR_shuffled"]
NEG = {"noise", "HR_shuffled"}


def rmsle_rows(y, p):  # 行ごとの二乗対数誤差（対の差の SE を行で出すため）
    return (np.log1p(p) - np.log1p(y)) ** 2


def fit_pred(Xtr, ytr, Xs, seed):
    m = HistGradientBoostingRegressor(max_iter=300, learning_rate=0.1, random_state=seed)
    m.fit(Xtr, np.log1p(ytr))
    return [np.clip(np.expm1(m.predict(X)), 1, None) for X in Xs]


def run(train, test, orig, cols, split_seed):
    """CV 設計書（Day163）どおり: ランダム5分割・log1p・元データは学習foldにだけ・下限1クリップ。
    返り値: fold ごとの RMSLE（5個）, 全体 CV, test 予測（5モデル平均）"""
    kf = KFold(5, shuffle=True, random_state=split_seed)
    oof, tp, fold_scores = np.zeros(len(train)), np.zeros(len(test)), []
    for tr, va in kf.split(train):
        tr_df = pd.concat([train.iloc[tr], orig], ignore_index=True)
        pv, pt = fit_pred(tr_df[cols], tr_df.Calories, [train.iloc[va][cols], test[cols]], split_seed)
        oof[va] = pv
        tp += pt / 5
        fold_scores.append(np.sqrt(rmsle_rows(train.Calories.values[va], pv).mean()))
    return np.array(fold_scores), np.sqrt(rmsle_rows(train.Calories.values, oof).mean()), tp


N_ORIG, N_TRAIN, N_TEST = 15000, 30000, 60000
WORLDS, REPEATS = range(4), range(3)
rec = []
for w in WORLDS:
    rng = np.random.default_rng(2000 + w)
    orig = make_original(N_ORIG, rng)
    pool = synth_from(orig, N_TRAIN + N_TEST, rng)
    orig, pool = add_candidates(orig, rng), add_candidates(pool, rng)
    train, test = pool.iloc[:N_TRAIN].reset_index(drop=True), pool.iloc[N_TRAIN:].reset_index(drop=True)
    res = {}
    for name in ["base"] + CANDS:
        cols = BASE if name == "base" else BASE + [name]
        res[name] = [run(train, test, orig, cols, r) for r in REPEATS]
    yt = test.Calories.values
    for c in CANDS:
        # (1) 1回の5分割: fold ごとの対の差
        fd = res[c][0][0] - res["base"][0][0]
        # (2) 3回繰り返しの15個の fold 差（重なりを無視した素朴な SE）
        fd15 = np.concatenate([res[c][r][0] - res["base"][r][0] for r in REPEATS])
        # (3) 繰り返しごとの全体 CV 差（3個）
        cvd = np.array([res[c][r][1] - res["base"][r][1] for r in REPEATS])
        # 真の効果: test 6万行で、3繰り返し×5モデルの平均予測どうしの差（行で対にした SE）
        pc = np.mean([res[c][r][2] for r in REPEATS], axis=0)
        pb = np.mean([res["base"][r][2] for r in REPEATS], axis=0)
        ec, eb = rmsle_rows(yt, pc), rmsle_rows(yt, pb)
        true_d = np.sqrt(ec.mean()) - np.sqrt(eb.mean())
        # 行ごとの差の SE を、平方根の前の平均の差から近似（デルタ法）
        true_se = (ec - eb).std(ddof=1) / np.sqrt(len(yt)) / (2 * np.sqrt(eb.mean()))
        rec.append(dict(world=w, cand=c,
                        d1=fd.mean(), se1=fd.std(ddof=1) / np.sqrt(5),
                        d15=fd15.mean(), se15=fd15.std(ddof=1) / np.sqrt(15),
                        d3=cvd.mean(), se3=cvd.std(ddof=1) / np.sqrt(3),
                        true_d=true_d, true_se=true_se))
    print(f"world {w} done {time.time()-T0:.0f}s", flush=True)

df = pd.DataFrame(rec)
df["truth_help"] = df.true_d < -2 * df.true_se          # 真に効く
df["truth_hurt"] = df.true_d > 2 * df.true_se           # 真に悪化
for k in ["1", "15", "3"]:
    df[f"adopt{k}"] = df[f"d{k}"] < -2 * df[f"se{k}"]
pd.set_option("display.width", 200)
print("\n=== 表1 候補ごとの 対の差（4ワールド平均）と 採用回数/4 ===")
g = df.groupby("cand", sort=False)
t1 = pd.DataFrame({
    "true_d": g.true_d.mean(), "true_se(平均)": g.true_se.mean(), "真に効く": g.truth_help.sum(),
    "d1": g.d1.mean(), "採用R1": g.adopt1.sum(),
    "d15": g.d15.mean(), "採用R2": g.adopt15.sum(),
    "d3": g.d3.mean(), "採用R3": g.adopt3.sum()})
print(t1.round(5).to_string())
print("\n=== 表2 SE の大きさ（4ワールド×7候補の中央値）===")
print(df[["se1", "se15", "se3", "true_se"]].median().round(6).to_string())
print("\n=== 表3 ルールごとの判定 対 真の効果（28 = 4ワールド×7候補）===")
for k, nm in [("1", "R1 1回5分割・fold差SE"), ("15", "R2 3回繰返し・15fold差SE(素朴)"), ("3", "R3 3回繰返し・CV差SE(n=3)")]:
    a = df[f"adopt{k}"]
    print(f"{nm:32s} 採用 {a.sum():2d} | 正しい採用 {(a & df.truth_help).sum():2d} | "
          f"偽採用(真に効かない) {(a & ~df.truth_help).sum():2d} | うち陰性対照 {(a & df.cand.isin(NEG)).sum():2d} | "
          f"見逃し {(~a & df.truth_help).sum():2d}")
print("\n=== 表4 CV差と真の効果の順位相関（候補×ワールド 28点、Spearman）===")
for k in ["1", "15", "3"]:
    print(f"d{k:3s} vs true_d: {df[f'd{k}'].rank().corr(df.true_d.rank()):.3f}")
print(f"\nwall {time.time()-T0:.0f}s")
df.to_csv("day164_records.csv", index=False)

# assert（数値の向きを固定する）
assert df[df.cand == "Keytel_kcal"].truth_help.sum() >= 3, "式そのものを渡す特徴量は真に効く"
assert (df.se15 < df.se1).mean() > 0.7, "15個の fold 差を独立扱いすると SE は小さく見える"
print("all asserts passed")

実行結果（実測）

# ============ 実行環境: macOS (M1, メモリ8GB)・Python 3.13.12・scikit-learn 1.9.1（新規 uv venv）
#              壁時計 1,191秒（ワールド0 はスワップ中で 546秒、1〜3 は各 180〜270秒）============
=== 表1 候補ごとの 対の差（4ワールド平均）と 採用回数/4 ===
              true_d  true_se(平均)  真に効く       d1  採用R1      d15  採用R2       d3  採用R3
cand                                                                                
Dur_x_HR    -0.00051      0.00010     4 -0.00063     2 -0.00048     4 -0.00048     4
log_Dur      0.00000      0.00000     0  0.00000     0  0.00000     0  0.00000     0
BMI          0.00014      0.00005     1  0.00014     0  0.00015     0  0.00015     0
HR_x_W       0.00006      0.00006     0  0.00016     0  0.00025     0  0.00025     0
Keytel_kcal -0.00044      0.00017     3 -0.00110     2 -0.00100     4 -0.00099     4
noise        0.00054      0.00005     0  0.00060     0  0.00055     0  0.00055     0
HR_shuffled  0.00048      0.00005     0  0.00050     0  0.00051     0  0.00051     0

=== 表2 SE の大きさ（4ワールド×7候補の中央値）===
se1        0.000200
se15       0.000127
se3        0.000083
true_se    0.000057

=== 表3 ルールごとの判定 対 真の効果（28 = 4ワールド×7候補）===
R1 1回5分割・fold差SE                 採用  4 | 正しい採用  4 | 偽採用(真に効かない)  0 | うち陰性対照  0 | 見逃し  4
R2 3回繰返し・15fold差SE(素朴)           採用  8 | 正しい採用  7 | 偽採用(真に効かない)  1 | うち陰性対照  0 | 見逃し  1
R3 3回繰返し・CV差SE(n=3)              採用  8 | 正しい採用  7 | 偽採用(真に効かない)  1 | うち陰性対照  0 | 見逃し  1

=== 表4 CV差と真の効果の順位相関（候補×ワールド 28点、Spearman）===
d1   vs true_d: 0.941
d15  vs true_d: 0.940
d3   vs true_d: 0.939

wall 1191s
all asserts passed

28件の内訳（真の効果で分類）: 真に効く 8件（Dur_x_HR 4・Keytel_kcal 3・BMI 1）、真に悪化 11件（noise 4・HR_shuffled 4・BMI 2・HR_x_W 1）、どちらとも言えない 9件（log_Dur 4 を含む）。

📊 結果の可視化

図1 候補ごとの CV の対の差（R2）と真の効果（4ワールド平均、×10⁻⁴）

図2 ルールごとの判定（28件のうち真に効く8件 + 誤判定）

図3 SE の大きさ（28件の中央値、×10⁻⁴）

🧭 決定 — 採否ルールと CV 設計書「決定6」

(1) 採否ルール: R2 を使い、R3 を照合に使う

(2) CV 設計書 決定6: 特徴量の採否（S5E5 用）

🪜 Step-by-Step 解説

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

• 5個の差の SE が揺れる理由: SE は「差のばらつき ÷ √個数」。5個だとばらつき自体がたまたま大きくも小さくも出る。5人に聞いたアンケートのばらつきを信じにくいのと同じ。今日の R1 の SE は 0〜0.00061 まで散らばった
• 15個を独立とみなすと楽観的になる理由: 3回の5分割は同じ3万行を並べ替えただけで、同じ行が検証側に3回出る。同じ人に3回聞いた答えを「3人分」と数えるようなもの。補正方法（Nadeau & Bengio 2003 の補正付き t 検定）はあるが今日は実装していない
• 単調変換: 大小の順番を変えない変換。1, 2, 3 の log は 0, 0.69, 1.10 で順番は同じ。木は「上か下か」しか見ないので同じ木になる。線形モデルや NN には間隔も効くので log に意味がある
• デルタ法: RMSLE は平均の平方根。平方根のグラフの傾き 1 ÷ (2×√値) を掛けて、平均の差の SE を RMSLE の差の SE に直す。小さな変化なら直線で近似してよい、という微分の考え方

🏆 Kaggleでの実践的な使い方

• 特徴量の採否は「3回繰り返し5分割の fold 差 −2SE」で判定し、判定表を残す。 1回の5分割は本当に効く特徴量の半分を見逃した。GBDT 1本が数分で済む Playground 級なら3倍の時間を払う価値がある
• 陰性対照（乱数列）を1本、同じ手順に通す。 乱数を採用したらルールか CV が壊れている
• GBDT に単調変換した列を足さない。 効くのは線形モデル・NN を混ぜるときと、目的変数側の変換（Day163 の log1p で −0.0049）
• ドメインの式は採用候補の最上位に置くが、CV の改善幅は割り引く。 今日は真の効果の約2.3倍に出た

⚠️ よくある誤解・ミス

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

• 実データ（S5E5）: Kaggle のログインとルール同意が必要なため未使用。Keytel_kcal がよく効くのはシミュレータの目的変数を Keytel 式で作ったからでもある。実物の目的変数の作り方は未確認
• モデルの seed 平均: HistGradientBoosting は1万行を超えると early stopping が自動で有効になり、判定用の10%を random_state で選ぶ。今日は random_state を分割 seed と同じにしたので、分割の揺れとモデルの揺れを分けていない
• Nadeau & Bengio の補正付き t 検定: R2 の楽観を補正する標準的な方法。今日は素朴な比較にとどめた
• 特徴量の組み合わせ（2つ同時・前向き選択）: 1つずつだけ。Dur_x_HR と Keytel_kcal の効果が足し算になるかは未測定
• ワールド数: 4ワールド（28件）。偽採用率を「1/20」の精度で言うには少ない。言えるのは件数そのもの（偽採用 R1 0・R2 1・R3 1）まで
• Day162 の選択バイアス仮説: 疑似ラベルを使う回に再検証（Day163 から持ち越し）

🚀 次のステップ

• 発展: R2 に Nadeau & Bengio の補正を入れ、偽採用と見逃しがどう動くかを同じ28件で測る。陰性対照を「害の無い列」に替えて偽採用率そのものを測る
• 次回予告: Day165 — Playground Series 実戦投入③。Dur_x_HR・Keytel_kcal を入れた状態で GBDT のハイパーパラメータを Optuna で探し、探索回数を増やすほど CV が楽観する「探索の過学習」を R2 と test 6万行で測る

Phase 6 の学習マップ

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