6.1 子集選擇

為什麼需要選擇變數？

想像你要搬家——你有一整間房子的家具，但新家只有一半大。你會怎麼做？當然是挑最重要的帶走。統計學習中的「子集選擇」(subset selection) 就是這個道理：從 \(p\) 個可能的預測變數中，挑出真正有用的那一組。

為什麼不全部留著？因為多餘的變數不但不會幫忙，還會增加模型複雜度、放大變異數、導致過擬合。好的模型不是變數最多的模型，而是在測試資料上表現最好的模型。

6.1.1 最佳子集選擇 (Best Subset Selection)

概念：暴力解法，全部試過

最佳子集選擇的邏輯非常直覺：把所有可能的變數組合都試一遍，挑最好的。就像去 buffet 每一道菜都吃一口，再決定最喜歡的組合。

具體來說，對於 \(p\) 個預測變數，我們要擬合所有包含恰好 \(k\) 個變數的模型（\(k = 0, 1, \dots, p\)），總共有 \(2^p\) 個可能的模型。演算法如下：

計算瓶頸：\(2^p\) 的詛咒

最佳子集選擇的問題是計算量隨 \(p\) 指數成長。當 \(p = 10\)，約需考慮 1,024 個模型；當 \(p = 20\)，超過一百萬個；當 \(p = 40\)，超過一兆個——即使用最快的電腦也跑不動。

# 最佳子集選擇的計算複雜度說明
import numpy as np

# 示範 2^p 的爆炸性成長
for p in [5, 10, 15, 20, 30, 40]:
    n_models = 2**p
    print(f"p = {p:2d}: 需要考慮 {n_models:>15,} 個模型")
# p =  5: 需要考慮              32 個模型
# p = 10: 需要考慮           1,024 個模型
# p = 15: 需要考慮          32,768 個模型
# p = 20: 需要考慮       1,048,576 個模型
# p = 30: 需要考慮   1,073,741,824 個模型
# p = 40: 需要考慮 1,099,511,627,776 個模型

用 Credit 資料集示範

try:
    from google.colab import drive
    drive.mount('/content/drive')
    DATA_PATH = '/content/drive/MyDrive/ISLP_data/'
except ImportError:
    DATA_PATH = '/tmp/'

import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from itertools import combinations
from ISLP import load_data
from sklearn.linear_model import LinearRegression

# 載入 Credit 資料集
credit = load_data('Credit')
# 只取數值型變數 + 將 region 轉為 dummy（仿課本 Figure 6.1）
X_full = pd.get_dummies(credit.drop('Balance', axis=1), drop_first=True)
y = credit['Balance']
p = X_full.shape[1]

# 只示範前 100 個組合（實際 2^11 = 2048 太多）
np.random.seed(1)
rss_list = []
r2_list = []
n_combos = min(100, 2**p)
sampled = [np.random.choice(p, size=k+1, replace=False)
           for k in range(p) for _ in range(min(30, n_combos//p))]

for cols in sampled:
    cols = list(set(cols))  # 去重
    if len(cols) == 0:
        continue
    model = LinearRegression().fit(X_full.iloc[:, cols], y)
    y_pred = model.predict(X_full.iloc[:, cols])
    rss = np.sum((y - y_pred)**2)
    tss = np.sum((y - np.mean(y))**2)
    r2 = 1 - rss / tss
    rss_list.append((len(cols), rss))
    r2_list.append((len(cols), r2))

# 找出每個 k 的最佳 RSS（紅色前沿）
rss_df = pd.DataFrame(rss_list, columns=['k', 'RSS'])
r2_df = pd.DataFrame(r2_list, columns=['k', 'R2'])
best_rss = rss_df.groupby('k')['RSS'].min().reset_index()
best_r2 = r2_df.groupby('k')['R2'].max().reset_index()

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4.5))
ax1.scatter(rss_df['k'], rss_df['RSS']/1e6, alpha=0.3, s=15, color='#58a6ff')
ax1.plot(best_rss['k'], best_rss['RSS']/1e6, 'r-o', lw=2, markersize=5, label='Best (red frontier)')
ax1.set_xlabel('Number of Predictors'); ax1.set_ylabel('RSS (×10⁶)')
ax1.set_title('RSS vs Number of Predictors')
ax1.legend(fontsize=8)

ax2.scatter(r2_df['k'], r2_df['R2'], alpha=0.3, s=15, color='#58a6ff')
ax2.plot(best_r2['k'], best_r2['R2'], 'r-o', lw=2, markersize=5, label='Best (red frontier)')
ax2.set_xlabel('Number of Predictors'); ax2.set_ylabel('R²')
ax2.set_title('R² vs Number of Predictors')
ax2.legend(fontsize=8)
plt.tight_layout(); plt.show()
print('Done: Best subset selection demo (like Figure 6.1)')

6.1.2 逐步選擇 (Stepwise Selection)

當 \(p\) 很大時，\(2^p\) 個模型根本跑不動。逐步選擇提供了一條「聰明捷徑」——不全部試，而是在模型空間中進行引導式搜尋。就像逛大賣場：你不會每條走道都走過，而是按照清單有方向地找。

前向逐步選擇 (Forward Stepwise Selection)

從零變數開始，每次加入一個「最能改善擬合」的變數，直到所有變數都被納入。

前向逐步選擇只需擬合 \(1 + p(p+1)/2\) 個模型，遠少於 \(2^p\)。當 \(p=20\) 時是 211 個 vs 1,048,576 個。

後向逐步選擇 (Backward Stepwise Selection)

從全模型（包含所有 \(p\) 個變數）開始，每次移除一個「最不有用」的變數。

用 sklearn 實作前向與後向選擇

try:
    from google.colab import drive
    drive.mount('/content/drive')
    DATA_PATH = '/content/drive/MyDrive/ISLP_data/'
except ImportError:
    DATA_PATH = '/tmp/'

import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from ISLP import load_data
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import cross_val_score

# 載入 Credit 資料
credit = load_data('Credit')
X = pd.get_dummies(credit.drop('Balance', axis=1), drop_first=True)
y = credit['Balance']
feat_names = list(X.columns)

def forward_stepwise(X, y, max_features=None):
    """前向逐步選擇：回傳每步選入的變數與對應的模型"""
    n, p = X.shape
    if max_features is None:
        max_features = min(p, n - 1)
    selected = []
    remaining = list(range(p))
    models = []
    scores = []
    for k in range(max_features):
        best_rss = float('inf')
        best_feat = None
        best_model = None
        for j in remaining:
            cols = selected + [j]
            model = LinearRegression().fit(X.iloc[:, cols], y)
            rss = np.sum((y - model.predict(X.iloc[:, cols]))**2)
            if rss < best_rss:
                best_rss = rss
                best_feat = j
                best_model = model
        selected.append(best_feat)
        remaining.remove(best_feat)
        models.append(best_model)
        scores.append(best_rss)
    return selected, models, scores

# 執行前向選擇
sel_idx, models, rss_vals = forward_stepwise(X, y, max_features=10)
print("變數選擇順序（前向）：")
for i, idx in enumerate(sel_idx):
    print(f"  Step {i+1}: {feat_names[idx]:<20s} RSS={rss_vals[i]/1e6:.2f}M")

# 視覺化：RSS 下降曲線
fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(range(1, len(rss_vals)+1), np.array(rss_vals)/1e6, 'b-o', lw=2, markersize=6)
ax.set_xlabel('Number of Predictors'); ax.set_ylabel('RSS (×10⁶)')
ax.set_title('Forward Stepwise: RSS vs Model Size (Credit Data)')
ax.grid(True, alpha=0.3)
plt.tight_layout(); plt.show()

try:
    from google.colab import drive
    drive.mount('/content/drive')
    DATA_PATH = '/content/drive/MyDrive/ISLP_data/'
except ImportError:
    DATA_PATH = '/tmp/'

import numpy as np
import pandas as pd
from ISLP import load_data
from sklearn.linear_model import LinearRegression

# 後向逐步選擇
credit = load_data('Credit')
X = pd.get_dummies(credit.drop('Balance', axis=1), drop_first=True)
y = credit['Balance']
feat_names = list(X.columns)

def backward_stepwise(X, y):
    """後向逐步選擇"""
    n, p = X.shape
    selected = list(range(p))
    removal_order = []
    rss_vals = []
    for k in range(p, 1, -1):
        worst_rss = float('inf')
        worst_feat = None
        for j in selected:
            cols = [c for c in selected if c != j]
            model = LinearRegression().fit(X.iloc[:, cols], y)
            rss = np.sum((y - model.predict(X.iloc[:, cols]))**2)
            if rss < worst_rss:
                worst_rss = rss
                worst_feat = j
        selected.remove(worst_feat)
        removal_order.append(worst_feat)
        rss_vals.append(worst_rss)
    # 最後剩一個變數
    model = LinearRegression().fit(X.iloc[:, selected], y)
    rss_vals.append(np.sum((y - model.predict(X.iloc[:, selected]))**2))
    return removal_order, rss_vals

removal, rss_bw = backward_stepwise(X, y)
print("變數移除順序（後向）：")
for i, idx in enumerate(removal):
    print(f"  Step {i+1}: 移除 {feat_names[idx]:<20s}")
print(f"  最後保留: {[feat_names[s] for s in [c for c in range(X.shape[1]) if c not in removal]]}")

混合方法 (Hybrid Approaches)

混合方法結合了前向與後向的優點：像前向選擇一樣逐步加入變數，但每加入一個新變數後，也檢查是否可以移除已選入但不再有用的變數。這種「邊加邊減」的策略更接近最佳子集選擇，同時保留了計算效率。

6.1.3 選擇最佳模型

子集選擇方法會產出一系列的模型 \(\mathcal{M}_0, \mathcal{M}_1, \dots, \mathcal{M}_p\)，現在的問題是：哪一個最好？因為訓練誤差（RSS / \(R^2\)）只會隨變數增加而改善，不能用來選。我們需要能估計測試誤差的指標。

\(C_p\)、AIC、BIC 與 Adjusted \(R^2\)

這些指標都在訓練誤差的基礎上加入懲罰項，懲罰模型複雜度（變數太多）。

四種指標的比較實作

try:
    from google.colab import drive
    drive.mount('/content/drive')
    DATA_PATH = '/content/drive/MyDrive/ISLP_data/'
except ImportError:
    DATA_PATH = '/tmp/'

import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from ISLP import load_data
from sklearn.linear_model import LinearRegression
from itertools import combinations
from math import comb

credit = load_data('Credit')
X = pd.get_dummies(credit.drop('Balance', axis=1), drop_first=True)
y = credit['Balance']
n, p = X.shape

# 用全模型估計 sigma^2
full_model = LinearRegression().fit(X, y)
sigma2_hat = np.sum((y - full_model.predict(X))**2) / (n - p - 1)
TSS = np.sum((y - np.mean(y))**2)

# 對每個 k 找隨機子集計算 Cp, AIC, BIC, AdjR2
np.random.seed(42)
results = []
for k in range(1, p+1):
    best_cp = float('inf')
    best_bic = float('inf')
    best_adjr2 = -float('inf')
    best_r2 = -float('inf')
    # 隨機取若干組合
    for _ in range(min(200, int(comb(p, k)) if k <= p//2 else 200)):
        cols = list(np.random.choice(p, size=k, replace=False))
        model = LinearRegression().fit(X.iloc[:, cols], y)
        rss = np.sum((y - model.predict(X.iloc[:, cols]))**2)
        r2 = 1 - rss / TSS
        cp = (rss + 2 * k * sigma2_hat) / n
        bic = (rss + np.log(n) * k * sigma2_hat) / n
        adjr2 = 1 - (rss/(n-k-1)) / (TSS/(n-1))
        if cp < best_cp: best_cp = cp
        if bic < best_bic: best_bic = bic
        if adjr2 > best_adjr2: best_adjr2 = adjr2
        if r2 > best_r2: best_r2 = r2
    results.append({'k': k, 'R2': best_r2, 'Cp': best_cp, 'BIC': best_bic, 'AdjR2': best_adjr2})

rdf = pd.DataFrame(results)

fig, axes = plt.subplots(1, 3, figsize=(14, 4))
axes[0].plot(rdf['k'], rdf['Cp'], 'b-o', lw=2, markersize=5)
axes[0].set_xlabel('Number of Predictors'); axes[0].set_ylabel('Cp')
axes[0].set_title('Cp (lower is better)')
axes[0].axvline(rdf.loc[rdf['Cp'].idxmin(), 'k'], color='red', ls='--', alpha=0.7)
axes[0].grid(True, alpha=0.3)

axes[1].plot(rdf['k'], rdf['BIC'], 'orange', marker='s', lw=2, markersize=5)
axes[1].set_xlabel('Number of Predictors'); axes[1].set_ylabel('BIC')
axes[1].set_title('BIC (lower is better)')
axes[1].axvline(rdf.loc[rdf['BIC'].idxmin(), 'k'], color='red', ls='--', alpha=0.7)
axes[1].grid(True, alpha=0.3)

axes[2].plot(rdf['k'], rdf['AdjR2'], 'g-^', lw=2, markersize=5)
axes[2].set_xlabel('Number of Predictors'); axes[2].set_ylabel('Adjusted R²')
axes[2].set_title('Adjusted R² (higher is better)')
axes[2].axvline(rdf.loc[rdf['AdjR2'].idxmax(), 'k'], color='red', ls='--', alpha=0.7)
axes[2].grid(True, alpha=0.3)

plt.tight_layout(); plt.show()
print(f"Cp 建議: {rdf.loc[rdf['Cp'].idxmin(), 'k']} 個變數")
print(f"BIC 建議: {rdf.loc[rdf['BIC'].idxmin(), 'k']} 個變數（傾向更小模型）")
print(f"AdjR² 建議: {rdf.loc[rdf['AdjR2'].idxmax(), 'k']} 個變數")
print(f"R² (僅供參考，不能用來選!): {rdf.loc[rdf['R2'].idxmax(), 'k']} 個變數 (=全模型)")

驗證集與交叉驗證

除了上述間接估計測試誤差的指標，我們也可以直接用驗證集或交叉驗證來估計測試誤差（詳見第五章）。這是最直接的方法：在不同資料分割上評估模型效能，選出平均驗證誤差最小的模型大小。

try:
    from google.colab import drive
    drive.mount('/content/drive')
    DATA_PATH = '/content/drive/MyDrive/ISLP_data/'
except ImportError:
    DATA_PATH = '/tmp/'

import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from ISLP import load_data
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import cross_val_score

credit = load_data('Credit')
X = pd.get_dummies(credit.drop('Balance', axis=1), drop_first=True)
y = credit['Balance']
n, p = X.shape

# 用前向逐步選擇的順序，對每個 k 做交叉驗證
def forward_order(X, y):
    selected = []
    remaining = list(range(X.shape[1]))
    for _ in range(min(X.shape[1], X.shape[0]-1)):
        best_rss = float('inf')
        best_feat = None
        for j in remaining:
            cols = selected + [j]
            model = LinearRegression().fit(X.iloc[:, cols], y)
            rss = np.sum((y - model.predict(X.iloc[:, cols]))**2)
            if rss < best_rss:
                best_rss = rss
                best_feat = j
        selected.append(best_feat)
        remaining.remove(best_feat)
    return selected

order = forward_order(X, y)
cv_errors = []
for k in range(1, len(order)+1):
    cols = order[:k]
    scores = cross_val_score(LinearRegression(), X.iloc[:, cols], y,
                             cv=5, scoring='neg_mean_squared_error')
    cv_errors.append(-scores.mean())

best_k = np.argmin(cv_errors) + 1
fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(range(1, len(cv_errors)+1), cv_errors, 'b-o', lw=2, markersize=6)
ax.axvline(best_k, color='red', ls='--', alpha=0.7,
           label=f'Best k={best_k} (min CV MSE)')
ax.set_xlabel('Number of Predictors'); ax.set_ylabel('5-Fold CV MSE')
ax.set_title('Cross-Validation for Model Size Selection (Credit Data)')
ax.legend(); ax.grid(True, alpha=0.3)
plt.tight_layout(); plt.show()
print(f"交叉驗證建議的最佳模型大小: k = {best_k}")

方法總比較

模型選擇準則比較

子集選擇的優缺點

應用場景

🔗 對 Hermes Agent 系統設計的啟發