Day 169 — 2026-10-01 | Playground Series 実戦投入（10/1 のバックフィル分・独立1本）— 「効かない特徴量を R2 はどれだけ通すか」。真の効果がほぼ0の候補3つ×36ワールド＝108件で、naive の −2SE は 11件（10.2%、名目 2.3% の約4.5倍）を採用した。ただし通した 11件の真の効果の合計は −1.69（わずかに得）で、採否を締めると効く列の得（正味 35.9 → 18.3 / 3.9）の方が大きく減る

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

Day168 の宿題。 Day168 では、R2 の SE（15個の fold 差を独立とみなす naive）は小さすぎ（被覆 0.82）、Nadeau & Bengio 補正（NB、2.18倍）は広すぎる（被覆 0.99）と分かった。しかし偽採用率は測れていなかった。入れた陰性対照 noise（真の効果 +8.5）と BMI（+2.1）は真に悪い列で、−2SE に近づきもしないため、「効かない列を誤って採用する頻度」の証拠にならなかった。決定10「採否は naive の −2SE のまま」は、偽採用を測らずに下した判断だった。今日はそこを埋める。

「真の効果がちょうど0の候補」は作れるか。 先に試した（パイロット2ワールド）。既存列の完全な複製 Age_dup = Age を足すと、HistGB は d_cv・sd_fold・d_true がすべてちょうど 0.0（2ワールドとも）。木は値の順序だけで分割するので、同じ列が1本増えても同じ分割が選ばれる。この場合 SE も 0 になり、差 < −2×0 は偽なので採用は絶対に起きない。つまり「ちょうど0」は SE の検査にならない（Day164 の log_Dur と同じ性質）。偽採用の検査には「予測は変わるが、平均の効果が0付近」の列が要る。

今日の候補（パイロットで 2ワールドの真の効果を見てから選んだ）。 既存列の積で、木が自力で作れるが明示するとわずかに分割が変わる列。

データ・モデル・R2 は Day168 と同じ（縮小シミュレータ、元1.5万・train3万・test6万、HistGB(early_stopping=False, random_state=0)、log1p 学習、下限1クリップ、元データは学習 fold だけに追加）。ワールドのシードも 3000+w で、Day168 の24ワールドと同じデータ（0〜23）に 12 ワールド（24〜35）を足した。差は ×10⁻⁴ で書く。数値は仮定の上の実測（実物の S5E5 ではない）。

用語整理:

❓ 問題

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

1. ほぼ0の3候補 108件のうち、naive の −2SE は何件採用するか。NB は何件か
2. 採用された列の真の効果は、平均して「害」と「得」のどちらか
3. 採否の閾値を −3SE に締める／naive に 1.5 倍を掛ける／NB にすると、効く列 Dur_x_HR（36ワールド中34で真に効く）の採用はどう減るか。正味はどう変わるか

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

1. Day168 の one_world(w) を、候補を環境変数 CANDS で切り替えられる形に直す（add に積の列を足す）
2. 36ワールドを 6プロセスに分け（OMP_NUM_THREADS=1）、CSV にまとめる
3. report で表A（真の効果・SE）、表B（z の SD・被覆）、表C（閾値ごとの採否）、表D（d_cv/SE の分布）、表F（正味）を出す

タスクC（決定・戦略）

結果から、決定6・10（採否は naive の −2SE）を見直すか決め、「決定11」を足す。「採った案／採らなかった案と理由／前提が変わったら」を付ける。

💡 ヒント

ヒント1（方向性）

• 偽採用は「採用の数」でなく、「採用した列が本当に悪化させた量」まで見る。件数だけ見ると、害の小さい誤りと大きい誤りが同じ1件になる
• 真の効果がちょうど0の列では SE も 0 になり得る。検査に使うなら、予測が少し変わる列を選ぶ

ヒント2（アプローチ）

• ほぼ0の候補はパイロット（2ワールド）で真の効果を見てから選ぶ。ただし選んだ後の本番は別ワールドを含め、パイロットの2ワールド（0, 1）は本番にも入る点を明記する
• 閾値を比べるときは同じ 144 行に当てはめる。閾値を選ぶ根拠にした行で評価するので、結論は「その行での比較」に留まる

ヒント3（コード骨格）

a = df.d_cv < -2 * df.naive                     # 採用
z = a & df.cand.isin(ZERO)                      # ほぼ0の列を通した件
loss = df.d_true[z].sum()                       # 正なら害、負なら（結果的に）得
gain = -df.d_true[a].sum()                      # 正味

✅ 模範解答

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

注意（「ほぼ0」の意味）。 3候補の真の効果は、ワールドをまたいだ平均が Age_x_Height −0.10（範囲 −2.26〜+2.19）、Temp_x_Age +0.31（−0.89〜+1.60）、Wt_x_Age −0.59（−2.66〜+1.66）、SD は約 0.95。ワールドごとに符号が変わる弱い効果であって、ちょうど0ではない。そのため「偽採用率 10.2%」は「効かない列を通す確率」でなく、「平均効果が小さく符号が揺れる列を通す確率」。真の効果がちょうど0の列（Age_dup）は SE=0 で採用が起きず、測れなかった。

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

day169.py（world・report の2モード）。Day168 のコードとの差は、CANDS 環境変数・add の列・report の中身だけ:

"""Day169: 真の効果が0付近の候補で、R2(15 fold差)の -2SE 採否の偽採用率を測る（Day168 の未達の消化）。
使い方:
  OMP_NUM_THREADS=1 CANDS=Dur_x_HR,Temp_x_Age,Wt_x_Age,Age_x_Height python day169.py world <w0> <w1> <out.csv>
  python day169.py report out1.csv out2.csv ...
"""
import os, sys, time
import numpy as np, pandas as pd
from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.model_selection import KFold

N_ORIG, N_TRAIN, N_TEST = 15000, 30000, 60000
BASE = ["Sex", "Age", "Height", "Weight", "Duration", "Heart_Rate", "Body_Temp"]
CAND_NAMES = os.environ.get("CANDS", "Dur_x_HR,Sex_x_Age,Age_dup,Wt_x_Age").split(",")
CANDS = {c: BASE + [c] for c in CAND_NAMES}
PARAMS = dict(learning_rate=0.1, max_iter=200, max_leaf_nodes=31, min_samples_leaf=20, l2_regularization=1e-4)
SPLITS = (1, 2, 3)                      # R2: 分割seed 3通り × 5分割 = 15 fold
K = 5


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〜168 と同じ
    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〜168 と同じ
    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(d, rng):
    d = d.copy(); d["Dur_x_HR"] = d.Duration * d.Heart_Rate
    d["Sex_x_Age"] = d.Sex * d.Age; d["Age_dup"] = d.Age; d["Wt_x_Age"] = d.Weight * d.Age
    d["Temp_x_Age"] = d.Body_Temp * d.Age; d["Age_x_Height"] = d.Age * d.Height; return d


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


def one_world(w):
    rng = np.random.default_rng(3000 + w)
    orig = make_original(N_ORIG, rng)
    pool = synth_from(orig, N_TRAIN + N_TEST, rng)
    orig, pool = add(orig, rng), add(pool, rng)
    train, test = pool.iloc[:N_TRAIN].reset_index(drop=True), pool.iloc[N_TRAIN:].reset_index(drop=True)
    yt = test.Calories.values

    def run(cols, split_seed):  # 各 fold の検証 RMSLE と、各 fold モデルの test RMSLE
        va_s, te_s = [], []
        for tr, va in KFold(K, shuffle=True, random_state=split_seed).split(train):
            d = pd.concat([train.iloc[tr], orig], ignore_index=True)
            m = HistGradientBoostingRegressor(early_stopping=False, random_state=0, **PARAMS)
            m.fit(d[cols], np.log1p(d.Calories))
            pv = np.clip(np.expm1(m.predict(train.iloc[va][cols])), 1, None)
            va_s.append(rmsle(train.Calories.values[va], pv))
            te_s.append(rmsle(yt, np.clip(np.expm1(m.predict(test[cols])), 1, None)))
        return np.array(va_s), np.array(te_s)

    res = {s: {"base": run(BASE, s), **{c: run(cols, s) for c, cols in CANDS.items()}} for s in SPLITS}
    out = []
    for c in CANDS:
        fd = np.concatenate([res[s][c][0] - res[s]["base"][0] for s in SPLITS])   # 15 個の fold 差
        td = np.concatenate([res[s][c][1] - res[s]["base"][1] for s in SPLITS])   # 15 個の「1モデルの test 差」
        out.append(dict(world=w, cand=c, d_cv=fd.mean(), sd_fold=fd.std(ddof=1), d_true=td.mean()))
    return out


if __name__ == "__main__":
    if sys.argv[1] == "world":
        w0, w1, path = int(sys.argv[2]), int(sys.argv[3]), sys.argv[4]
        rows, T0 = [], time.time()
        for w in range(w0, w1):
            rows += one_world(w); pd.DataFrame(rows).to_csv(path, index=False)
            print(f"world {w} done {time.time()-T0:.0f}s", flush=True)

    if sys.argv[1] == "report":
        df = pd.concat([pd.read_csv(p) for p in sys.argv[2:]], ignore_index=True)
        J, ratio = 15, N_TRAIN / K / (N_TRAIN - N_TRAIN / K)
        df["naive"] = df.sd_fold / np.sqrt(J)
        df["nb"] = df.sd_fold * np.sqrt(1 / J + ratio)
        df["err"] = df.d_cv - df.d_true
        ZERO = [c for c in df.cand.unique() if c != "Dur_x_HR"]
        # 表A〜D（真の効果・z の SD・被覆・採否・d_cv/SE の分布）と、下の表F:
        for lab, mult, th in (("naive -2.0SE", 1.0, 2.0), ("naive -2.5SE", 1.0, 2.5), ("naive -3.0SE", 1.0, 3.0),
                              ("x1.5 -2SE", 1.5, 2.0), ("NB -2SE", None, 2.0)):
            se = df.nb if mult is None else df.naive * mult
            a = df.d_cv < -th * se
            z = a & df.cand.isin(ZERO); p = a & (df.cand == "Dur_x_HR")
            print(f"{lab:13s} 正味 {-df.d_true[a].sum()*1e4:7.2f} | ほぼ0 採用{int(z.sum()):2d}/{(df.cand.isin(ZERO)).sum()} "
                  f"損益 {df.d_true[z].sum()*1e4:6.2f} | Dur_x_HR 採用{int(p.sum()):2d}/36 得 {-df.d_true[p].sum()*1e4:6.2f}")

（表A〜D の print は Day168 の report と同じ形で、day169.py の全文は scratchpad に保存した。ここでは差分の核だけを載せた。）

実行結果（実測）

# ============ 実行環境: macOS (M1, メモリ8GB)・Python 3.13.12・scikit-learn 1.9.1・numpy 2.5.3・pandas 3.0.6（新規 uv venv）
# world モード: 36ワールドを 6プロセス × 6ワールドに分け OMP_NUM_THREADS=1 で並列（壁時計 804秒）。パイロット 2ワールドは単独 171秒/ワールド（6モデル）
worlds: 36  行: 144

=== 表A 候補ごと（×1e-4）===
Age_x_Height  true平均  -0.10 範囲  -2.26~   2.19 true<0 22/36 | cv平均  -0.25 err SD 1.81 | SE naive中央 1.35 NB中央 2.95
Dur_x_HR      true平均  -2.14 範囲  -5.27~   1.19 true<0 34/36 | cv平均  -2.57 err SD 2.42 | SE naive中央 1.38 NB中央 3.01
Temp_x_Age    true平均   0.31 範囲  -0.89~   1.60 true<0 11/36 | cv平均   0.10 err SD 1.51 | SE naive中央 1.18 NB中央 2.57
Wt_x_Age      true平均  -0.59 範囲  -2.66~   1.66 true<0 25/36 | cv平均  -0.79 err SD 2.02 | SE naive中央 1.42 NB中央 3.09

=== 表B z=err/SE の SD と 被覆（1.96SE）===
Age_x_Height  n= 36 SD(z): naive 1.37 NB 0.63 | 被覆 naive 0.86 NB 1.00
Dur_x_HR      n= 36 SD(z): naive 1.88 NB 0.86 | 被覆 naive 0.75 NB 0.97
Temp_x_Age    n= 36 SD(z): naive 1.30 NB 0.60 | 被覆 naive 0.86 NB 1.00
Wt_x_Age      n= 36 SD(z): naive 1.48 NB 0.68 | 被覆 naive 0.89 NB 1.00
ほぼ0の3候補    n=108 SD(z): naive 1.37 NB 0.63 | 被覆 naive 0.87 NB 1.00

=== 表C -2SE 採否（採用 = d_cv < -2SE）===
[naive] Age_x_Height 3/36（害2） Dur_x_HR 16/36（真<0 15・害1） Temp_x_Age 2/36（害2） Wt_x_Age 6/36（真<0 4・害2）  ほぼ0計 11/108 率0.102 害6
[nb]    Age_x_Height 0/36 Dur_x_HR 2/36 Temp_x_Age 0/36 Wt_x_Age 0/36                                            ほぼ0計 0/108
[x1.5]  Age_x_Height 0/36 Dur_x_HR 8/36 Temp_x_Age 0/36 Wt_x_Age 1/36（害1）                                      ほぼ0計 1/108 率0.009

=== 表D ほぼ0の候補 d_cv/SE の分布（名目 P(<-2)=0.023）===
naive  SD 1.26 平均 -0.23 P(<-2) 0.102 P(<-1.5) 0.139 P(>2) 0.037
nb     SD 0.58 平均 -0.10 P(<-2) 0.000 P(<-1.5) 0.000 P(>2) 0.000
真の効果を引いた z (naive) P(<-2) 0.093  真の効果の SD 0.95 平均 -0.13

=== 表F 採否ルール別の正味（×1e-4。採用した列の真の効果の合計の符号反転。大きいほど良い）===
naive -2.0SE  正味 35.90 | ほぼ0 採用11/108 損益 -1.69 | Dur_x_HR 採用16/36 得 34.21 | 不採用の Dur_x_HR の真の得 42.88
naive -2.5SE  正味 21.69 | ほぼ0 採用 6/108 損益 -0.60 | Dur_x_HR 採用11/36 得 21.09 | 不採用の Dur_x_HR の真の得 56.01
naive -3.0SE  正味 18.27 | ほぼ0 採用 1/108 損益  0.98 | Dur_x_HR 採用 8/36 得 19.25 | 不採用の Dur_x_HR の真の得 57.85
x1.5 -2SE     正味 18.27 | ほぼ0 採用 1/108 損益  0.98 | Dur_x_HR 採用 8/36 得 19.25 | 不採用の Dur_x_HR の真の得 57.85
NB -2SE       正味  3.94 | ほぼ0 採用 0/108 損益  0.00 | Dur_x_HR 採用 2/36 得  3.94 | 不採用の Dur_x_HR の真の得 73.16
偽採用率 0.102 Wilson95% [0.058, 0.173]
ワールドをまたいだ相関: d_cv(Wt_x_Age) と d_cv(Age_x_Height) の相関 0.24
CV差の符号が真と逆: Age_x_Height 15/36, Temp_x_Age 16/36, Wt_x_Age 13/36

読み方: 「損益」は採用した列の真の効果の合計（正なら害、負なら得）。「採用しなかった Dur_x_HR の真の得」が、各ルールが捨てた分。

タスクC（決定）

(1) 見立て: R2 の naive −2SE は、ほぼ0の列を約1割通す（名目の約4.5倍）。ただし通した分の害は小さく、締めると効く列の取りこぼしの方が大きいので、採否のルールは変えない。

(2) CV 設計書に書き足す「決定11: 偽採用の見積もり」（S5E5 用）

📊 結果の可視化

偽採用（ほぼ0の列を通した件数）は 11 → 6 → 1 → 0 と減るが、効く列の得を捨てる量の方が大きい。

🪜 Step-by-Step 解説

Step 1: 「ちょうど0」の候補を試し、検査に使えないと確かめる

Age_dup = Age は d_cv・sd_fold・d_true がすべて 0.0（2ワールド）。木は順序だけを見るので、複製列は分割を変えない。SE が 0 なので d < −2SE が起きない。「偽採用 0 件」は SE の良さを示さない。

Step 2: ほぼ0の候補を選ぶ

既存列の積5つを2ワールドで試し、真の効果が小さい3つ（Temp_x_Age・Wt_x_Age・Age_x_Height）を選んだ。Sex_x_Age（−2.35, −1.31）は効きすぎなので外した。選び方の注意: パイロットで「小さい」と見て選んだので、本番36ワールドでの範囲は −2.66〜+2.19 と広がった。「ほぼ0」は2ワールドの印象で、36ワールドでは符号が揺れる弱い効果だった。

Step 3: SE の足りなさを再確認する（表B）

3候補・108件で、naive の z の SD は 1.37、被覆 0.87。Day168 の 3候補（SD 1.67、被覆 0.82）より少し良いが、同じ向き。NB は SD 0.63・被覆 1.00 で広すぎる。

Step 4: 偽採用率（表C・D）

naive −2SE の採用は 11/108（10.2%）。名目 2.3% ならば約 2.5 件のはずで、4.5倍。ただし d_cv/SE の平均は −0.23 と、真の効果の平均（−0.13）に近い。採用の過剰は「SE が小さい」ことと「真の効果が揺れる」ことの両方から来ている（真の効果を引いた z でも P(<−2)=9.3%）。分けてはいない。

Step 5: 採用した列の害を数える（表F）

採用 11件の真の効果は −2.49〜+1.66。害は 6件（+0.01〜+1.66）、得は 5件。合計は −1.69。つまり通した弱い列は、平均では test をわずかに良くした。1件ごとの害の上限は +1.66 で、Dur_x_HR の平均の得（−2.14）の約 0.8 倍。

Step 6: 閾値を締めた場合の損得

Dur_x_HR が真に効くワールドは 34/36。naive −2SE は 16件採用して 34.2 の得を取る。−3SE は 8件で 19.3。捨てた得は 42.9 → 57.9。偽採用を 11件から 1件に減らして得られるのは損益 +2.7（−1.69 → +0.98）だが、捨てるのは 15.0 の得。この比較は 1シミュレータ・効く列1種類で、効く列の割合が違えば変わる。

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

名目 2.3% の意味。 正規分布で「−2 より小さい」ことは約 2.3%（片側）。真の効果が0で、SE が正しければ、−2SE を下回る確率はそこまで下がる。実測 10.2% は、SE が小さいこと（Day168）と、本当に少し効いている／悪いことの重なり。

なぜ「害」と「得」が混ざるか。 試験で「たまたま点が高かった」生徒を選ぶと、次の試験では平均に戻る。CV で「たまたま良く見えた」弱い列を選ぶと、本当の効果は小さく、符号も不確か。ただ、元々少し効く列が混ざっていれば、選んだ集団の平均はわずかに得になる。

Wilson 区間。 108件中11件という少ない件数での割合の幅。[5.8%, 17.3%] なので、「約1割、ただし 6〜17% の幅」と読む。

🏆 Kaggleでの実践的な使い方

• CV で −2SE を超えた特徴量は「本当に効く」と「弱い列がたまたま」の混合。今日の見積もりでは約1割の弱い列が通る。通しても害は +1.7×10⁻⁴ 以内だったが、これは今日のシミュレータの話
• 完全な複製列は木モデルでは何も変えない。特徴量の重複を「検査用の陰性対照」に使うと、SE が0になって検査にならない。木には効かない変換（単調変換）も同じ（Day164 の log_Dur）
• 候補を1つずつ足す運用では、取りこぼしと偽採用の大きさを揃えて比べる。今日は偽採用 1 件分の害（平均 +0.8）より、取りこぼす 1 件分の得（平均 −2.1）の方が大きい。実物（S5E5 本物のデータ）では未確認

⚠️ よくある誤解・ミス

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

• 真の効果がちょうど0で SE が0にならない候補: 作れなかった（複製・単調変換は SE=0）。ランダムに予測を変えるが平均効果0の列（例: 学習ごとに seed を変えた乱数列は noise と同じく悪化する）は試していない。本当の「偽採用率」の下限・上限は未測定
• 世界の実物との比較: 1シミュレータ・HistGB・36ワールド・候補4つ（うち効く列は1種類）。Nadeau & Bengio（2003, Machine Learning 52:239–281）は、Springer の論文ページが認証ページへリダイレクトされ、今日も開けていない（式は記憶）。実物の S5E5 は未使用
• 候補プールの構成を変えた場合: 効く列1：ほぼ0の列3 の割合のみ。効く列が 1割以下のときの損益の逆転点は計算していない
• 行対デルタ法（Day164）・R3 との比較: 今日は fold 差の SE（R2）のみ
• パイロットの2ワールドが本番に重なる: ワールド 0, 1 は選定にも評価にも使った。除外した場合の数字は出していない
• Day168 の 24ワールドとの合算: 候補が違う（BMI・noise は今日は入れていない）ので、合算していない

🚀 次のステップ

• 発展: 学習のたびに予測が少し変わるが平均効果が0の列（例: 木の max_features で列抽出を入れたときの seed 違い）を候補にして、偽採用率の下限を測る。効く列の割合を 1:3 から 1:9 に変えて、決定11 の損益の逆転点を数える
• 次回予告: （この回は独立したバックフィル。state の次回予告「Day 167 — 最終提出の選び方」は変更しない。同予告の持ち越しのうち、Day168 の「偽採用率は未測定」は今日で部分的に消化。なお予告の「Day 167」は実際には Day170 以降の番号になる）

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