§6.1 子集選擇法（Subset Selection）

📚 理論核心

🔢 三種子集選擇法

1. 最佳子集選擇（Best Subset Selection）

暴力法。對 \(p\) 個預測變數，嘗試所有 \(2^p\) 種可能組合，對每種組合擬合模型，選出最佳者。

2. 向前逐步選擇（Forward Stepwise Selection）

從空模型開始，每次加入「最有幫助」的一個變數，逐步建構。總共只需擬合約 \(1 + \sum_{k=0}^{p-1} (p-k) \approx 1 + p(p+1)/2\) 個模型。

3. 向後逐步選擇（Backward Stepwise Selection）

從全模型開始，每次踢出「最沒貢獻」的一個變數。要求 \(n > p\)（樣本數大於變數數）。

📏 模型選擇指標

選出 \(\mathcal{M}_0, \mathcal{M}_1, \ldots, \mathcal{M}_p\) 後，用什麼標準挑最終模型？

🐍 實作程式碼

# §6.1 子集選擇 — 完整 Python 實作
# 可在 Google Colab 直接執行

import numpy as np
import pandas as pd
import itertools
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import cross_val_score

# ============================================
# 1. 讀取資料（ISLP 課本 Credit 資料集）
# ============================================
url = "https://www.statlearning.com/s/Credit.csv"
credit = pd.read_csv(url)
credit = pd.get_dummies(credit, columns=['Student', 'Gender', 'Married', 'Ethnicity'], drop_first=True)
credit = credit.drop(columns=['Unnamed: 0'])

y = credit['Balance'].values
X_all = credit.drop(columns=['Balance'])
feature_names = X_all.columns.tolist()
X_all = X_all.values
n, p = X_all.shape

print(f"樣本數 n={n}, 預測變數 p={p}")
print(f"若做最佳子集 = 2^{p} = {2**p:,} 個模型 → 不可行，改用向前逐步")

# ============================================
# 2. 向前逐步選擇（Forward Stepwise）
# ============================================
selected = []           # 已選入的變數 index
remaining = list(range(p))  # 尚未選入的變數 index
models_log = []         # 記錄每個 k 的最佳模型

for k in range(p):
    best_rss = float('inf')
    best_var = None
    best_model = None
    
    for j in remaining:
        trial_vars = selected + [j]
        X_trial = X_all[:, trial_vars]
        model = LinearRegression().fit(X_trial, y)
        y_pred = model.predict(X_trial)
        rss = np.sum((y - y_pred) ** 2)
        
        if rss < best_rss:
            best_rss = rss
            best_var = j
            best_model = model
    
    selected.append(best_var)
    remaining.remove(best_var)
    models_log.append({
        'k': k + 1,
        'vars': [feature_names[i] for i in selected],
        'RSS': best_rss,
        'coef_count': len(selected) + 1  # +1 for intercept
    })

# ============================================
# 3. 計算模型選擇指標
# ============================================
sigma2_hat = models_log[-1]['RSS'] / (n - p - 1)  # 用全模型估計 σ²

for m in models_log:
    d = m['coef_count']  # 參數數（含截距）
    m['Cp'] = (m['RSS'] + 2 * d * sigma2_hat) / n
    m['AIC'] = n * np.log(m['RSS'] / n) + 2 * d
    m['BIC'] = n * np.log(m['RSS'] / n) + np.log(n) * d
    m['Adj_R2'] = 1 - (m['RSS'] / (n - d)) / (np.sum((y - np.mean(y))**2) / (n - 1))

# ============================================
# 4. 輸出結果
# ============================================
results_df = pd.DataFrame(models_log)
results_df = results_df.round({'RSS': 0, 'Cp': 1, 'AIC': 0, 'BIC': 0, 'Adj_R2': 4})

print("\n📊 向前逐步選擇結果：")
print(results_df[['k', 'Cp', 'AIC', 'BIC', 'Adj_R2']].to_string(index=False))

# 找出各指標的最優模型
best_cp_k = results_df.loc[results_df['Cp'].idxmin(), 'k']
best_bic_k = results_df.loc[results_df['BIC'].idxmin(), 'k']
best_adjr2_k = results_df.loc[results_df['Adj_R2'].idxmax(), 'k']

print(f"\n✅ C_p 最優：{best_cp_k} 個變數")
print(f"✅ BIC 最優：{best_bic_k} 個變數（最簡約）")
print(f"✅ Adj-R² 最優：{best_adjr2_k} 個變數")

# 顯示 BIC 最優模型的變數
best_bic_vars = results_df.loc[results_df['BIC'].idxmin(), 'vars']
print(f"\n🏆 BIC 選出的關鍵變數：{best_bic_vars}")

# ============================================
# 5. 交叉驗證確認
# ============================================
print("\n🔬 5-fold CV 確認 BIC 最優模型：")
X_best = X_all[:, [feature_names.index(v) for v in best_bic_vars]]
model_best = LinearRegression()
cv_scores = cross_val_score(model_best, X_best, y, cv=5, scoring='neg_mean_squared_error')
rmse_cv = np.sqrt(-cv_scores.mean())
print(f"5-Fold CV RMSE = {rmse_cv:.1f}")

# 對比全模型的 CV
model_full = LinearRegression()
cv_scores_full = cross_val_score(model_full, X_all, y, cv=5, scoring='neg_mean_squared_error')
rmse_cv_full = np.sqrt(-cv_scores_full.mean())
print(f"全模型 (p={p}) 5-Fold CV RMSE = {rmse_cv_full:.1f}")
print(f"➡️ 變數從 {p} 降到 {best_bic_k}，CV RMSE 變化: {rmse_cv:.1f} vs {rmse_cv_full:.1f}")

🌍 應用場景

⚖️ 優缺點分析

🆚 與其他方法的比較

🧠 自我內化