7.8 Lab: 非線性建模（Lab: Non-Linear Modeling）

ISLP 7.8 pp. 309324 ★★★☆☆ 35
Python GAM pygam
7.7 GAM

7.8.1 與階梯函數（Polynomial Regression & Step Functions）

本節用 Wage 資料集示範如何擬合非線性模型。我們先從多項式迴歸開始——用 age 預測 wage，用 ANOVA 檢定選擇最佳多項式次數，再用階梯函數（step function）將連續變數離散化。

THEORY 多項式迴歸：y = 0 + 1x + 2x + ... + dxd。可用 ANOVA F-test 比較 nested models（ISLP 7.1）

SCENE ：薪資與年齡的非線性關係

工資隨年齡的變化不是一條直線——年輕時快速增長，中年達高峰，之後趨緩。多項式迴歸用一條彎曲的線捕捉這個「倒 U 型」模式，比直線更能反映真實的勞動市場。

# 多項式迴歸 + ANOVA 模型選擇
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
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import statsmodels.api as sm
from ISLP import load_data
from ISLP.models import ModelSpec as MS, poly, summarize
from statsmodels.stats.anova import anova_lm

Wage = load_data('Wage')
y = Wage['wage']
age = Wage['age']

# 擬合 4 次多項式
poly_age = MS([poly('age', degree=4)]).fit(Wage)
M = sm.OLS(y, poly_age.transform(Wage)).fit()
summarize(M)
print('---')
# ANOVA 比較 1~5 次多項式
models = [MS([poly('age', degree=d)]) for d in range(1, 6)]
Xs = [model.fit_transform(Wage) for model in models]
fits = [sm.OLS(y, X).fit() for X in Xs]
anova_res = anova_lm(*fits)
print(anova_res)
# ANOVA 結果：cubic (d=3) 足夠，4 次和 5 次不顯著
4 次多項式係數： coef std err t P>|t| intercept 111.7036 0.729 153.283 0.000 poly[0] 447.0679 39.915 11.201 0.000 poly[1] -478.3158 39.915 -11.983 0.000 poly[2] 125.5217 39.915 3.145 0.002 poly[3] -77.9112 39.915 -1.952 0.051 ANOVA 顯示 cubic 已是「夠用且不冗餘」的選擇

階梯函數（Step Functions）

階梯函數將連續變數切割成區間，每個區間配一個常數預測值——等於把「彎曲的關係」用階梯近似。

# 階梯函數：用 pd.qcut 將 age 切成 4 個區間
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
import matplotlib
matplotlib.use('Agg')
import statsmodels.api as sm
from ISLP import load_data
from ISLP.models import summarize

Wage = load_data('Wage')
y = Wage['wage']
age = Wage['age']

cut_age = pd.qcut(age, 4)
step_model = sm.OLS(y, pd.get_dummies(cut_age)).fit()
summarize(step_model)
coef std err t P>|t| (17.999, 33.75] 94.1584 1.478 63.692 0.0 (33.75, 42.0] 116.6608 1.470 79.385 0.0 (42.0, 51.0] 119.1887 1.416 84.147 0.0 (51.0, 80.0] 116.5717 1.559 74.751 0.0 高峰在 42~51 歲（$119,189），清楚展示「倒 U 型」薪資走勢

多項式迴歸

多項式迴歸

7.8.2 樣條（Splines）

樣條是分段多項式，在節點（knot）處平滑連接。核心優勢是局部彈性——改變一個點只影響附近區間，不會波及整條曲線。

# B-spline：3 個內部節點 + bs() helper
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
import matplotlib
matplotlib.use('Agg')
import statsmodels.api as sm
from ISLP import load_data
from ISLP.models import ModelSpec as MS, bs, summarize
from ISLP.transforms import BSpline

Wage = load_data('Wage')
y = Wage['wage']
age = Wage['age']

# B-spline 基礎矩陣（3 個內部節點 → 6 basis functions）
bs_obj = BSpline(internal_knots=[25, 40, 60], intercept=True).fit(age)
bs_matrix = bs_obj.transform(age)
print(f'B-spline shape: {bs_matrix.shape}')  # (3000, 7)

# 用 bs() helper + OLS 擬合
bs_age = MS([bs('age', internal_knots=[25, 40, 60])])
Xbs = bs_age.fit_transform(Wage)
M_bs = sm.OLS(y, Xbs).fit()
summarize(M_bs)

# 用 df= 自動選節點
bs_auto = BSpline(df=6).fit(age)
print(f'Auto knots: {bs_auto.internal_knots_}')  # [33.75, 42.0, 51.0] — 分位數
B-spline shape: (3000, 7) 6 個 spline coefficients（不含 intercept）: bs[0] 3.981 (p=0.751) — 不顯著 bs[1] 44.631 (p=0.000) — 中年族群薪資大幅上升 bs[2] 62.839 (p=0.000) — 45~50 歲達高峰 bs[3] 55.991 (p=0.000) bs[4] 50.688 (p=0.000) bs[5] 16.606 (p=0.385) — 尾端不顯著

自然樣條（Natural Splines）

自然樣條在邊界區域強制線性，避免邊界處的過度彎曲。

# 自然樣條（Natural Spline）
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
import matplotlib
matplotlib.use('Agg')
import statsmodels.api as sm
from ISLP import load_data
from ISLP.models import ModelSpec as MS, ns, summarize

Wage = load_data('Wage')
y = Wage['wage']
age = Wage['age']

ns_age = MS([ns('age', df=5)]).fit(Wage)
M_ns = sm.OLS(y, ns_age.transform(Wage)).fit()
summarize(M_ns)
coef std err t P>|t| intercept 60.475 4.708 12.844 0.000 ns(age, df=5)[0] 61.527 4.709 13.065 0.000 ns(age, df=5)[1] 55.691 5.717 9.741 0.000 ns(age, df=5)[2] 46.818 4.948 9.463 0.000 ns(age, df=5)[3] 83.204 11.918 6.982 0.000 ns(age, df=5)[4] 6.877 9.484 0.725 0.468 ← 最後一項不顯著
方法節點行為邊界行為自由度適用情境
無節點，全域高次時爆炸d+1簡單非線性
節點處 C 連續無約束K+4 (cubic)內部彈性需求高
同上 + 邊界線性強制線性K (內部節點)邊界預測穩定
全節點 + 懲罰由 控制等效 df自動平滑，免選節點

7.8.3 平滑樣條與 GAM（Smoothing Splines & GAMs）

平滑樣條在每個觀測點放節點，用懲罰項控制彎曲度。我們用 pygam 套件實作。

THEORY 平滑樣條：最小化 (yi f(xi)) + f(t) dt。 控制平滑： 0 內插所有點， 直線（ISLP 7.5） 
# 平滑樣條：不同  值的擬合效果
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
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from ISLP import load_data
from pygam import LinearGAM, s as s_gam

Wage = load_data('Wage')
y = Wage['wage']
age = Wage['age']
X_age = np.asarray(age).reshape((-1, 1))
age_grid = np.linspace(age.min(), age.max(), 100)

# 比較不同 lam (lambda) 值
fig, ax = plt.subplots(figsize=(8, 6))
ax.scatter(age, y, facecolor='gray', alpha=0.3)
for lam in np.logspace(-2, 6, 5):
    gam = LinearGAM(s_gam(0, lam=lam)).fit(X_age, y)
    ax.plot(age_grid, gam.predict(age_grid),
            label=f'{lam:.1e}', linewidth=2)
ax.set_xlabel('Age'); ax.set_ylabel('Wage')
ax.legend(title='lambda')
plt.savefig('/tmp/smooth_spline_lam.png', dpi=80, bbox_inches='tight')
plt.close()
print('Plot saved. Lambda range: 1e-2 to 1e6 shows wiggly to linear.')

也可用自由度（df）來控制，更直覺：

# 用自由度控制平滑程度（df+1 因為平滑樣條預設含 intercept）
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
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from ISLP import load_data
from pygam import LinearGAM, s as s_gam
from ISLP.pygam import approx_lam, degrees_of_freedom

Wage = load_data('Wage')
y = Wage['wage']
age = Wage['age']
X_age = np.asarray(age).reshape((-1, 1))
age_grid = np.linspace(age.min(), age.max(), 100)

gam = LinearGAM(s_gam(0, lam=0.6)).fit(X_age, y)
age_term = gam.terms[0]

fig, ax = plt.subplots(figsize=(8, 6))
ax.scatter(X_age, y, facecolor='gray', alpha=0.2)
for df in [1, 3, 4, 8, 15]:
    lam = approx_lam(X_age, age_term, df + 1)
    age_term.lam = lam
    gam.fit(X_age, y)
    ax.plot(age_grid, gam.predict(age_grid),
            label=f'df={df}', linewidth=2)
ax.set_xlabel('Age'); ax.set_ylabel('Wage')
ax.legend(title='Degrees of freedom')
plt.savefig('/tmp/smooth_spline_df.png', dpi=80, bbox_inches='tight')
plt.close()
print('df=4 gives a smooth, interpretable fit.')

多元 GAM（Additive Models with Several Terms）

GAM 的真正威力在於同時處理多變數，並對每個變數獨立指定平滑形式。以下對 age 和 year 用平滑樣條、education 用類別變數：

# pygam GAM：s(age) + s(year) + f(education)
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
import matplotlib
matplotlib.use('Agg')
from ISLP import load_data
from pygam import LinearGAM, s as s_gam, f as f_gam, l as l_gam
from ISLP.pygam import approx_lam, anova as anova_gam

Wage = load_data('Wage')
y = Wage['wage']
age = Wage['age']
high_earn = Wage['high_earn'] = (y > 250)

Xgam = np.column_stack([age,
                        Wage['year'],
                        Wage['education'].cat.codes])

# 擬合完整 GAM
gam_full = LinearGAM(s_gam(0) + s_gam(1, n_splines=7) + f_gam(2, lam=0))
gam_full = gam_full.fit(Xgam, y)

# 設定 age 和 year 各 4 個自由度
age_term = gam_full.terms[0]
age_term.lam = approx_lam(Xgam, age_term, df=4 + 1)
year_term = gam_full.terms[1]
year_term.lam = approx_lam(Xgam, year_term, df=4 + 1)
gam_full = gam_full.fit(Xgam, y)

# ANOVA 比較 year 的三種設定
gam_0 = LinearGAM(age_term + f_gam(2, lam=0)).fit(Xgam, y)
gam_linear = LinearGAM(age_term + l_gam(1, lam=0) + f_gam(2, lam=0)).fit(Xgam, y)
ar = anova_gam(gam_0, gam_linear, gam_full)
print('ANOVA for year effect:')
print(ar)
# 結論：year 線性顯著 (p~0.002)，但非線性不顯著 (p~0.435)
ANOVA for year effect: deviance df deviance_diff df_diff F pvalue 0 3714362.366 2991.004 NaN NaN NaN NaN 1 3696745.823 2990.005 17616.543 0.999 14.265 0.002 2 3693142.930 2987.007 3602.894 2.998 0.972 0.436 結論：year 的線性效果顯著（p=0.002），不需要非線性（p=0.436）

Logistic GAM：預測高薪者（wage > $250K）

# Logistic GAM：預測 wage > 250
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
import matplotlib
matplotlib.use('Agg')
from ISLP import load_data
from pygam import LogisticGAM, s as s_gam, f as f_gam, l as l_gam
from ISLP.pygam import approx_lam

Wage = load_data('Wage')
y = Wage['wage']
age = Wage['age']
high_earn = Wage['high_earn'] = (y > 250)

Xgam = np.column_stack([age,
                        Wage['year'],
                        Wage['education'].cat.codes])

# 先檢查 education 分布
ct = pd.crosstab(Wage['high_earn'], Wage['education'])
print('Cross-tabulation:')
print(ct)

# <HS 類別沒有任何高薪者 → 排除後重擬合
only_hs = Wage['education'] == '1. < HS Grad'
Wage_ = Wage.loc[~only_hs]
Xgam_ = np.column_stack([Wage_['age'],
                         Wage_['year'],
                         Wage_['education'].cat.codes - 1])
high_earn_ = Wage_['high_earn']

gam_logit = LogisticGAM(s_gam(0) + s_gam(1) + f_gam(2, lam=0))
gam_logit.fit(Xgam_, high_earn_)
print(f'Logistic GAM fitted. Accuracy: {gam_logit.accuracy(Xgam_, high_earn_):.4f}')
Cross-tabulation: 1. < HS Grad 2. HS Grad 3. Some College 4. College Grad 5. Advanced Degree False 267 967 643 680 415 True 0 14 19 38 45 <HS 類別 0 人高薪 → 排除後模型更合理

7.8.4 局部迴歸（Local Regression / LOWESS）

局部迴歸在每個預測點 x 附近取「鄰居視窗」，對視窗內做加權線性迴歸。視窗大小由 span 參數控制。

# LOWESS：不同 span 值的局部線性迴歸
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
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import statsmodels.api as sm
from ISLP import load_data

Wage = load_data('Wage')
y = Wage['wage']
age = Wage['age']

lowess = sm.nonparametric.lowess
fig, ax = plt.subplots(figsize=(8, 6))
ax.scatter(age, y, facecolor='gray', alpha=0.3)

for span in [0.2, 0.5]:
    fitted = lowess(y, age, frac=span, return_sorted=True)
    ax.plot(fitted[:, 0], fitted[:, 1],
            label=f'span={span}', linewidth=2)

ax.set_xlabel('Age'); ax.set_ylabel('Wage')
ax.legend()
plt.savefig('/tmp/lowess_span.png', dpi=80, bbox_inches='tight')
plt.close()
print('LOWESS: span=0.2 (wiggly), span=0.5 (smooth).')

SCENE ：房價估價中的局部迴歸

房地產估價師常用局部迴歸：要評估一間 80 坪房子的價值，只看附近（80~100 坪）的成交紀錄做加權平均，而不是用全市 20~200 坪的所有資料跑一條直線。這就是「局部」的精髓——近鄰比遠親更有參考價值。

非線性方法核心機制彈性來源優點缺點
高次多項式次數 d簡單、可微分全域震盪、邊界失控
分段常數區間數極簡、無假設不連續、資訊損失
分段多項式 + 節點節點數/位置局部彈性、連續需選節點
樣條 + 邊界線性內部節點數邊界穩定邊界較不彈性
全節點 + 懲罰 (等效 df)自動平滑、免選節點計算量較大
鄰域加權線性span（視窗比例）極度局部邊界稀疏、外推差
GAMfj(xj) 加總每變數獨立控制可解釋、靈活需選各項平滑度
「你不需要在直線和彎曲線之間二選一——你可以對 age 用彎曲線、對 year 用直線、對 education 用階梯。GAM 的本質就是把複雜問題拆成可獨立控制的簡單部分。」

Hermes 系統設計啟發：模組化的彈性控制

GAM 的設計哲學——「每個變數獨立選擇平滑度」——對 AI agent 架構有深刻啟發。在 Hermes 的 skill 系統中，我們也應該對不同任務類型採用不同的「彈性等級」：

就像 GAM 不加假設地給 year 多餘自由度只是在浪費參數，agent 系統也該對不同子任務動態分配推理資源，而不是所有任務都給完整的 toolset 和 token budget。