Day 105 — 外れ値の影響と対策 — Huber Loss・量子回帰

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

🧩 外れ値対策の3つの視点

🎯 問題

中古マンション価格予測（回帰）を模した合成データ。学習データには5%の外れ値（データ入力ミス・超高級物件を想定し価格が3〜8倍に膨らんだ行）を混入させ、holdoutは完全にクリーン（本来知りたい「普通の物件」の集団）とする。

import numpy as np
import pandas as pd

def make_price_df(seed, n, outlier_frac=0.05):
    rng_state = np.random.get_state()
    np.random.seed(seed)
    size_sqm = np.random.normal(70, 20, n).clip(20, 200)
    age_years = np.random.exponential(15, n).clip(0, 80)
    distance_station_km = np.random.exponential(2, n).clip(0.1, 15)
    num_rooms = np.random.randint(1, 6, n)
    renovated = np.random.binomial(1, 0.3, n)

    base_price = (
        3000000
        + 45000 * size_sqm
        - 20000 * age_years
        - 80000 * distance_station_km
        + 500000 * renovated
        + np.random.normal(0, 600000, n)
    )
    price = base_price.clip(500000, None)

    n_outliers = int(n * outlier_frac)
    if n_outliers > 0:
        outlier_idx = np.random.choice(n, n_outliers, replace=False)
        price[outlier_idx] = price[outlier_idx] * np.random.uniform(3, 8, n_outliers)

    df = pd.DataFrame({'SizeSqm': size_sqm, 'AgeYears': age_years,
                        'DistanceStationKm': distance_station_km,
                        'NumRooms': num_rooms, 'Renovated': renovated, 'Price': price})
    np.random.set_state(rng_state)
    return df

FEATURES = ['SizeSqm', 'AgeYears', 'DistanceStationKm', 'NumRooms', 'Renovated']

df_train = make_price_df(seed=42, n=2000, outlier_frac=0.05)   # 学習データ（外れ値5%混入）
df_holdout = make_price_df(seed=999, n=1000, outlier_frac=0.0)  # クリーンなholdout

X, y = df_train[FEATURES].values, df_train['Price'].values
X_hold, y_hold = df_holdout[FEATURES].values, df_holdout['Price'].values

タスク

📐 損失関数の形を比較

横軸は残差 r（実際値−予測値）、縦軸は損失の大きさ。同じ残差でも損失関数によって「罰の増え方」がまったく違うことを図解する。

📊 MAE・係数・カバレッジの比較（実行結果）

実際にコードを実行して得た値（乱数シード固定・環境により多少前後する）。バーはそれぞれのチャート内での相対値。

クリーンholdoutでのMAE比較

係数の歪み: Renovated（真の係数=500,000）

係数の歪み: NumRooms（真の係数=0・価格と無関係）

epsilonスイープ（HuberRegressor、holdout MAE）

outlier_frac スイープ（holdoutは常にクリーン固定）

P10-P90 予測区間の経験的カバレッジ

💡 ヒント

from sklearn.linear_model import LinearRegression, HuberRegressor
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.metrics import mean_absolute_error, mean_squared_error

lr = LinearRegression().fit(X, y)
hub = HuberRegressor(epsilon=___, max_iter=1000).fit(X, y)
gbr_q50 = GradientBoostingRegressor(loss=___, alpha=0.5, n_estimators=200, max_depth=3, random_state=42).fit(X, y)

for name, model in [('LinearRegression', lr), ('Huber', hub), ('Quantile median', gbr_q50)]:
    pred = model.predict(___)  # ← X_hold
    mae = mean_absolute_error(y_hold, pred)
    rmse = mean_squared_error(y_hold, pred) ** 0.5
    print(name, mae, rmse)

# 量子回帰で予測区間
gbr_q10 = GradientBoostingRegressor(loss='quantile', alpha=___, n_estimators=200, max_depth=3, random_state=42).fit(X, y)
gbr_q90 = GradientBoostingRegressor(loss='quantile', alpha=___, n_estimators=200, max_depth=3, random_state=42).fit(X, y)
coverage = np.mean((y_hold >= gbr_q10.predict(X_hold)) & (y_hold <= ___))

✅ 模範解答

import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression, HuberRegressor
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.metrics import mean_absolute_error, mean_squared_error

# ====== データ準備（省略: 問題文のmake_price_dfをそのまま使用） ======
df_train = make_price_df(seed=42, n=2000, outlier_frac=0.05)
df_holdout = make_price_df(seed=999, n=1000, outlier_frac=0.0)
X, y = df_train[FEATURES].values, df_train['Price'].values
X_hold, y_hold = df_holdout[FEATURES].values, df_holdout['Price'].values

# ====== タスク2: 3モデルを学習しholdoutで比較 ======
lr = LinearRegression().fit(X, y)
hub = HuberRegressor(epsilon=1.35, max_iter=1000).fit(X, y)
gbr_q50 = GradientBoostingRegressor(loss='quantile', alpha=0.5, n_estimators=200, max_depth=3, random_state=42).fit(X, y)

for name, model in [('LinearRegression', lr), ('HuberRegressor(eps=1.35)', hub), ('GBR quantile(alpha=0.5)', gbr_q50)]:
    pred = model.predict(X_hold)
    mae = mean_absolute_error(y_hold, pred)
    rmse = mean_squared_error(y_hold, pred) ** 0.5
    print(f"{name:28s}: MAE={mae:,.0f}  RMSE={rmse:,.0f}")

print("\nLinearRegression coef:", np.round(lr.coef_, 1))
print("HuberRegressor coef:  ", np.round(hub.coef_, 1))
print("真の係数:              [45000.0 -20000.0 -80000.0 0.0 500000.0]")

# ====== タスク3: epsilon・outlier_frac のスイープ ======
print("\n--- epsilon スイープ（HuberRegressor） ---")
for eps in [1.1, 1.35, 1.5, 2.0, 5.0]:
    h = HuberRegressor(epsilon=eps, max_iter=1000).fit(X, y)
    pred = h.predict(X_hold)
    print(f"epsilon={eps:4.2f}: MAE={mean_absolute_error(y_hold, pred):,.0f}")

print("\n--- outlier_frac スイープ（holdoutは常にクリーン固定） ---")
for frac in [0.0, 0.05, 0.15, 0.30]:
    df_tr = make_price_df(seed=42, n=2000, outlier_frac=frac)
    Xf, yf = df_tr[FEATURES].values, df_tr['Price'].values
    lr_f = LinearRegression().fit(Xf, yf)
    hub_f = HuberRegressor(epsilon=1.35, max_iter=1000).fit(Xf, yf)
    mae_lr = mean_absolute_error(y_hold, lr_f.predict(X_hold))
    mae_hub = mean_absolute_error(y_hold, hub_f.predict(X_hold))
    print(f"outlier_frac={frac:.2f}: LinearRegression MAE={mae_lr:,.0f}  |  Huber MAE={mae_hub:,.0f}")

# ====== タスク4: 量子回帰で予測区間を作りカバレッジを確認 ======
gbr_q10 = GradientBoostingRegressor(loss='quantile', alpha=0.1, n_estimators=200, max_depth=3, random_state=42).fit(X, y)
gbr_q90 = GradientBoostingRegressor(loss='quantile', alpha=0.9, n_estimators=200, max_depth=3, random_state=42).fit(X, y)

pred_q10, pred_q50, pred_q90 = gbr_q10.predict(X_hold), gbr_q50.predict(X_hold), gbr_q90.predict(X_hold)
coverage = np.mean((y_hold >= pred_q10) & (y_hold <= pred_q90))
print(f"\nP10-P90 経験的カバレッジ: {coverage:.4f}（目標: 約0.80）")
print(f"平均区間幅: {np.mean(pred_q90 - pred_q10):,.0f}")

出力イメージ（実際に実行して確認した値。乱数・環境により多少前後する）:

LinearRegression            : MAE=1,209,810  RMSE=1,393,665
HuberRegressor(eps=1.35)    : MAE=472,350  RMSE=593,536
GBR quantile(alpha=0.5)     : MAE=497,668  RMSE=626,191

LinearRegression coef: [  63067.3  -20091.8  -67398.3  -81853.5 1152766.7]
HuberRegressor coef:   [  45514.3  -17840.6  -78983.8  -20275.2  415242.2]
真の係数:              [45000.0 -20000.0 -80000.0 0.0 500000.0]

--- epsilon スイープ（HuberRegressor） ---
epsilon=1.10: MAE=472,547
epsilon=1.35: MAE=472,350
epsilon=1.50: MAE=471,840
epsilon=2.00: MAE=471,080
epsilon=5.00: MAE=1,031,850

--- outlier_frac スイープ（holdoutは常にクリーン固定） ---
outlier_frac=0.00: LinearRegression MAE=472,530  |  Huber MAE=475,303
outlier_frac=0.05: LinearRegression MAE=1,209,810  |  Huber MAE=472,350
outlier_frac=0.15: LinearRegression MAE=3,885,645  |  Huber MAE=478,855
outlier_frac=0.30: LinearRegression MAE=7,764,981  |  Huber MAE=651,325

P10-P90 経験的カバレッジ: 0.8020（目標: 約0.80）
平均区間幅: 2,012,673

🪜 Step-by-Step 解説

13モデルを外れ値混入データで学習し、クリーンなholdoutで評価する

lr = LinearRegression().fit(X, y)
hub = HuberRegressor(epsilon=1.35, max_iter=1000).fit(X, y)
gbr_q50 = GradientBoostingRegressor(loss='quantile', alpha=0.5, n_estimators=200, max_depth=3, random_state=42).fit(X, y)

2係数を比較し、外れ値が生む「幻の相関」を確認する

print("LinearRegression coef:", lr.coef_)
print("HuberRegressor coef:  ", hub.coef_)

3epsilon・outlier_fracを振って、ロバスト性の限界と代償を確認する

for eps in [1.1, 1.35, 1.5, 2.0, 5.0]:
    h = HuberRegressor(epsilon=eps, max_iter=1000).fit(X, y)
    ...
for frac in [0.0, 0.05, 0.15, 0.30]:
    ...

4量子回帰で予測区間を作り、カバレッジを検証する

pred_q10, pred_q90 = gbr_q10.predict(X_hold), gbr_q90.predict(X_hold)
coverage = np.mean((y_hold >= pred_q10) & (y_hold <= pred_q90))

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

数式を見たら

🏆 Kaggleでの実践的な使い方

よく使われるコンペカテゴリ: ☑ 表形式データ（Tabular） / ☐ 自然言語処理（NLP） / ☐ 画像認識（CV） / ☑ 時系列（Time Series）

• LightGBM/XGBoostのロバスト目的関数: LightGBMはobjective='huber'やobjective='quantile'、objective='fair'を直接指定できる。XGBoostにもreg:pseudohubererrorがある。保険金額・不動産価格・売上予測など外れ値の多い実データでは、前処理でクリッピングする代わりにこれらを試すのが定石の一つ
• M5 Forecasting - Uncertainty（Kaggle過去コンペ）: 点予測ではなく複数パーセンタイルの予測区間を提出させるコンペで、まさに今日学んだ量子回帰（Pinball Loss/WSPLでの評価）が中心的な手法だった
• 保険金請求額予測: 少数の巨大請求（火災・大事故）が学習データに混ざる業界では、二乗誤差ベースのモデルがその少数例に引っ張られやすい。Huber Lossや対数変換（Phase 3で学んだ目的変数変換）と組み合わせて使われる
• 前処理（クリッピング・対数変換）との使い分け: Phase 3 Day 072で学んだ「外れ値をクリップ/対数変換する」方法はデータ側を変える発想、今日のHuber/量子回帰は損失関数（モデル側）を変える発想。どちらか一方でなく併用されることも多い

⚠️ よくある誤解・ミス

🚀 次のステップ

• 発展: 3つの量子回帰モデル（P10・P50・P90）の予測値を全holdout行で比較し、quantile crossing（P10がP50を上回るような矛盾）が実際に起きていないか確認してみる。もし起きていた場合、np.sortで3つの予測値を強制的に単調化する後処理を試す
• 次回予告: Day 106「不均衡データ対策（SMOTE・class_weight・閾値調整）」 — 外れ値と並んで実務・Kaggle頻出の「データの偏り」問題として、少数派クラスへの対処法を学ぶ

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