7.1 多項式迴歸（Polynomial Regression）

從直線到曲線：為什麼線性不夠用？

到目前為止，整本書的焦點都在 線性模型 上。線性迴歸簡單、解讀容易、推論方便，但它在預測能力上有先天限制——現實世界的關係往往是彎的，不是直的。

想像你要用「年齡」預測「薪資」：20 歲和 40 歲的薪資差距，會跟 40 歲和 60 歲的差距一樣嗎？顯然不會。薪資通常呈現倒 U 型：年輕時期快速成長，中年到達高峰，接近退休時平緩甚至下滑。一條直線根本無法捕捉這種曲線。

用生活化比喻：線性迴歸像用 一把直尺 量身高——你只能畫一條直線穿過資料。多項式迴歸則像拿 一卷軟尺，它可以彎曲來貼合曲線，彎曲的程度由多項式的次數 \( d \) 控制。

實戰：用多項式迴歸預測薪資

我們用 ISLP 內建的 Wage 資料集來示範。這個資料包含 3,000 名男性勞工的薪資和人口統計資訊。我們用「年齡」來預測「薪資」，從 1 次（直線）逐漸增加到 4 次（曲線），觀察擬合的變化。

# 7.1 多項式迴歸：Wage 資料集的示範
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import numpy as np
from ISLP import load_data
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures

# 讀取 Wage 資料
wage = load_data('Wage')
age = wage['age'].values.reshape(-1, 1)
wage_val = wage['wage'].values

# 比較不同次數的多項式擬合
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
degrees = [1, 3, 4]
titles = ['Degree 1 (線性)', 'Degree 3 (三次)', 'Degree 4 (四次)']

age_grid = np.linspace(age.min(), age.max(), 200).reshape(-1, 1)

for ax, d, title in zip(axes, degrees, titles):
    poly = PolynomialFeatures(degree=d, include_bias=False)
    X_poly = poly.fit_transform(age)
    model = LinearRegression().fit(X_poly, wage_val)
    X_grid_poly = poly.transform(age_grid)
    y_pred = model.predict(X_grid_poly)
    ax.scatter(age[::20], wage_val[::20], alpha=0.3, s=4, color='#58a6ff')
    ax.plot(age_grid, y_pred, color='#f85149', linewidth=2, label=f'd={d}')
    ax.set(xlabel='Age', ylabel='Wage ($1,000)', title=title)
    ax.legend()

plt.tight_layout()
plt.show()
print(f"Degree-4 R-squared = {LinearRegression().fit(PolynomialFeatures(degree=4, include_bias=False).fit_transform(age), wage_val).score(PolynomialFeatures(degree=4, include_bias=False).fit_transform(age), wage_val):.3f}")
print(f"Degree-1 R-squared = {LinearRegression().fit(age, wage_val).score(age, wage_val):.3f}")

觀察結果：degree=1（直線）的 \( R^2 \) 相當低——一條直線解釋不了薪資的倒 U 型軌跡。degree=4 的曲線明顯更貼近資料的整體趨勢，\( R^2 \) 顯著提升。但注意，更高次數不代表更好——下一節會討論如何選擇最佳次數。

擬合的變異：信賴區間怎麼算？

看到圖上的虛線了嗎？那兩條虛線代表估計的 95% 信賴區間——是 \( 2 \times \) 點標準誤的範圍。

假設我們在某個年齡 \( x_0 \) 上計算了擬合值：

這個擬合值的變異有多大？最小平方法會給出每個 \( \hat{\beta}_j \) 的變異估計，以及係數之間的共變異。如果 \( \hat{\mathbf{C}} \) 是 \( 5 \times 5 \) 的係數共變異矩陣，且定義向量 \( \boldsymbol{\ell}_0^T = (1, x_0, x_0^2, x_0^3, x_0^4) \)，則擬合值的變異為：

把這個開根號就是 點標準誤。對每個 \( x_0 \) 重複此計算，就得到整條信賴區間曲線。

# 信賴區間的近似計算示範
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import numpy as np
from ISLP import load_data
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures

wage = load_data('Wage')
age = wage['age'].values.reshape(-1, 1)
wage_val = wage['wage'].values

# 擬合 degree-4 多項式
poly = PolynomialFeatures(degree=4, include_bias=False)
X_poly = poly.fit_transform(age)
model = LinearRegression().fit(X_poly, wage_val)

# 計算標準誤（近似：來自殘差變異）
age_grid = np.linspace(age.min(), age.max(), 200).reshape(-1, 1)
X_grid_poly = poly.transform(age_grid)
y_pred = model.predict(X_grid_poly)

# 殘差變異估計
n = len(wage_val)
residuals = wage_val - model.predict(X_poly)
sigma2_hat = np.sum(residuals**2) / (n - 5)  # 5 = number of coefficients

# 近似信賴區間（大樣本下寬度近似常數）
se_fit = np.sqrt(sigma2_hat) * np.ones_like(y_pred)

fig, ax = plt.subplots(figsize=(10, 6))
ax.scatter(age[::15], wage_val[::15], alpha=0.3, s=5, color='#58a6ff')
ax.plot(age_grid, y_pred, color='#f85149', linewidth=2, label='Degree-4 fit')
ax.fill_between(age_grid.ravel(), y_pred - 2*se_fit, y_pred + 2*se_fit,
                alpha=0.15, color='#58a6ff', label='Approx 95% CI')
ax.set(xlabel='Age', ylabel='Wage ($1,000)')
ax.legend()
plt.show()

課本 Figure 7.1（左）顯示：信賴區間在資料密集的區域（30–60 歲）較窄，在兩端（極年輕和極年長）較寬——因為兩端的樣本點少，估計自然較不確定。這跟選舉民調的邏輯一樣：樣本少，誤差範圍就大。

多項式 + 邏輯回歸：預測「高收入族群」

觀察 Figure 7.1（左）會發現一個有趣現象：薪資分布似乎來自 兩個不同族群——一群年收入低於 $250K 的「普通勞工」，和一小群超過 $250K 的「高收入者」（整個資料集只有 79 人）。

我們可以把薪資轉成二元變數（是否 > $250K），然後用 多項式邏輯回歸 來預測這個二元結果：

# 多項式邏輯回歸：預測薪資 > $250K 的機率
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import numpy as np
from ISLP import load_data
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import PolynomialFeatures

wage = load_data('Wage')
age = wage['age'].values.reshape(-1, 1)
high_earner = (wage['wage'].values > 250).astype(int)

# 擬合 degree-4 多項式邏輯回歸
poly = PolynomialFeatures(degree=4, include_bias=False)
X_poly = poly.fit_transform(age)
logreg = LogisticRegression(C=1e10, solver='lbfgs', max_iter=2000)
logreg.fit(X_poly, high_earner)

# 預測網格上的機率
age_grid = np.linspace(age.min(), age.max(), 200).reshape(-1, 1)
X_grid_poly = poly.transform(age_grid)
probs = logreg.predict_proba(X_grid_poly)[:, 1]

fig, ax = plt.subplots(figsize=(10, 6))
jitter_high = np.random.uniform(-0.02, 0.02, size=np.sum(high_earner == 1))
jitter_low  = np.random.uniform(-0.02, 0.02, size=np.sum(high_earner == 0))
ax.scatter(age[high_earner == 1], high_earner[high_earner == 1] + jitter_high,
           alpha=0.5, s=8, color='#f85149', label='> $250K')
ax.scatter(age[high_earner == 0], high_earner[high_earner == 0] + jitter_low,
           alpha=0.2, s=8, color='#8b949e', label='<= $250K')
ax.plot(age_grid, probs, color='#58a6ff', linewidth=2.5, label='Predicted probability')
ax.set(xlabel='Age', ylabel='Pr(Wage > $250K)')
ax.legend()
ax.set_ylim(-0.05, 1.05)
plt.show()

# 統計高收入者分布
n_high = np.sum(high_earner == 1)
print(f"High earners (>$250K): {n_high} / {len(high_earner)} ({n_high/len(high_earner)*100:.1f}%)")

課本 Figure 7.1（右）顯示：成為高收入者的機率在中年（約 40–50 歲）達到高峰，之後逐漸下降。但信賴區間相當寬——尤其是在右側（年長端）——原因是 樣本極度不平衡：雖然總樣本數 \( n = 3{,}000 \) 很大，但只有 79 人屬於高收入，導致係數估計變異很大。

多項式次數的選擇：太少？太多？

多項式迴歸最關鍵的選擇就是次數 \( d \)。用圖形化的方式幫助理解：

# 比較不同次數：使用交叉驗證來選擇最佳次數
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import numpy as np
from ISLP import load_data
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
from sklearn.model_selection import cross_val_score
from sklearn.pipeline import make_pipeline

wage = load_data('Wage')
age = wage['age'].values.reshape(-1, 1)
wage_val = wage['wage'].values

max_degree = 12
cv_means = []
cv_stds = []

for d in range(1, max_degree + 1):
    pipe = make_pipeline(PolynomialFeatures(degree=d, include_bias=False),
                         LinearRegression())
    scores = cross_val_score(pipe, age, wage_val, cv=10,
                             scoring='neg_mean_squared_error')
    cv_means.append(-scores.mean())
    cv_stds.append(scores.std())

fig, ax = plt.subplots(figsize=(9, 5))
ax.errorbar(range(1, max_degree + 1), cv_means, yerr=cv_stds,
            marker='o', color='#58a6ff', capsize=4)
best_d = np.argmin(cv_means) + 1
ax.axvline(best_d, color='#f85149', linestyle='--', alpha=0.5,
           label=f'Best d={best_d} (test MSE={cv_means[best_d-1]:.1f})')
ax.set(xlabel='Polynomial Degree', ylabel='10-Fold CV Mean MSE')
ax.legend()
plt.show()
print(f"Optimal degree d={best_d}, CV test MSE={cv_means[best_d-1]:.1f}")
print(f"d=1 test MSE={cv_means[0]:.1f}, d=12 test MSE={cv_means[11]:.1f}")

交叉驗證結果清楚展示偏差-變異權衡：低次數（\( d=1 \)）測試誤差高（高偏差、低變異），適中次數（約 \( d=4 \sim 7 \)）測試誤差最低，後續次數測試誤差再次上升（低偏差、高變異）——這就是過度擬合。

優缺點對照

類似方法比較

多項式迴歸的「基函數」觀點

多項式迴歸可以被視為 基函數方法（basis function approach） 的特例（§7.3）。我們選擇一組固定的函數 \( b_1(x), b_2(x), \dots, b_K(x) \)，然後擬合：

對多項式迴歸而言，基函數就是 \( b_j(x) = x^j \)。這個視角的好處是：一旦我們理解多項式只是其中一種基函數選擇，後續的階梯函數、樣條、小波等非線性方法就都掉進同一個框架——只是換了一組不同的 \( b_j(\cdot) \) 而已。整個第 7 章就是在討論不同基函數選擇的權衡。