8.3 Lab: 決策樹方法實戰

🧪 Lab 概覽

本 Lab 涵蓋五個核心實作：分類樹、迴歸樹、Bagging、隨機森林、Boosting 與貝氏加法迴歸樹（BART）。每一節都從資料前處理開始，經模型擬合、交叉驗證剪枝，到測試集評估，形成完整的工作流程。準備好你的 Python 環境，我們開始動手！

8.3.1 分類樹實作

資料準備與初始擬合

我們使用 Carseats 資料集，將連續變數 Sales 轉換為二元變數 High（銷量 > 8 為 "Yes"），然後以 DecisionTreeClassifier 擬合分類樹。

# Block 1: 分類樹初始擬合與訓練準確率
import numpy as np
import pandas as pd
from ISLP import load_data
from ISLP.models import ModelSpec as MS
from sklearn.tree import DecisionTreeClassifier as DTC
from sklearn.metrics import accuracy_score

Carseats = load_data('Carseats')
High = np.where(Carseats.Sales > 8, "Yes", "No")

# 建立模型矩陣（排除 Sales）
model = MS(Carseats.columns.drop('Sales'), intercept=False)
D = model.fit_transform(Carseats)
feature_names = list(D.columns)
X = np.asarray(D)

# 擬合分類樹：entropy 準則，最大深度 3
clf = DTC(criterion='entropy', max_depth=3, random_state=0)
clf.fit(X, High)

# 訓練準確率
train_acc = accuracy_score(High, clf.predict(X))
print(f"訓練準確率: {train_acc:.4f}")

訓練誤差率約 27.25%。sklearn 的決策樹直接使用 one-hot 編碼後的虛擬變數進行分割，而非像理論所述「自然地分割類別變數的 levels」。這點實際使用時需注意。

殘差 Deviance 計算

殘差 deviance 是分類樹的擬合品質指標，定義為終端節點中各類別比例乘積對數的負兩倍總和：

其中 \(n_{mk}\) 是第 \(m\) 個終端節點中第 \(k\) 類的觀測數，\(\hat{p}_{mk}\) 是該節點的預測機率。deviance 越小代表訓練擬合越好。

# Block 2: 計算殘差 deviance
import numpy as np
import pandas as pd
from ISLP import load_data
from ISLP.models import ModelSpec as MS
from sklearn.tree import DecisionTreeClassifier as DTC
from sklearn.metrics import log_loss

Carseats = load_data('Carseats')
High = np.where(Carseats.Sales > 8, "Yes", "No")
model = MS(Carseats.columns.drop('Sales'), intercept=False)
D = model.fit_transform(Carseats)
X = np.asarray(D)

clf = DTC(criterion='entropy', max_depth=3, random_state=0)
clf.fit(X, High)

resid_dev = np.sum(log_loss(High, clf.predict_proba(X)))
print(f"殘差 deviance: {resid_dev:.4f}")

查看決策樹結構（文字版）

使用 export_text 可以將整棵樹印成易讀的文字格式，包括每個分割的條件與終端節點的類別權重。

# Block 3: 匯出樹的文字結構
import numpy as np
import pandas as pd
from ISLP import load_data
from ISLP.models import ModelSpec as MS
from sklearn.tree import DecisionTreeClassifier as DTC, export_text

Carseats = load_data('Carseats')
High = np.where(Carseats.Sales > 8, "Yes", "No")
model = MS(Carseats.columns.drop('Sales'), intercept=False)
D = model.fit_transform(Carseats)
feature_names = list(D.columns)
X = np.asarray(D)

clf = DTC(criterion='entropy', max_depth=3, random_state=0)
clf.fit(X, High)

print(export_text(clf, feature_names=feature_names, show_weights=True))

從樹結構可以看出，ShelveLoc（貨架位置）是最重要的分割變數：貨架位置「好」的商店，銷量普遍較高。價格與廣告支出則在次級分割中扮演關鍵角色。

驗證集評估

直接使用訓練準確率會高估模型表現。我們使用 ShuffleSplit 進行一次性 hold-out 驗證：

# Block 4: ShuffleSplit 驗證
import numpy as np
import pandas as pd
from ISLP import load_data
from ISLP.models import ModelSpec as MS
from sklearn.tree import DecisionTreeClassifier as DTC
import sklearn.model_selection as skm

Carseats = load_data('Carseats')
High = np.where(Carseats.Sales > 8, "Yes", "No")
model = MS(Carseats.columns.drop('Sales'), intercept=False)
D = model.fit_transform(Carseats)
X = np.asarray(D)

clf = DTC(criterion='entropy', max_depth=3, random_state=0)
validation = skm.ShuffleSplit(n_splits=1, test_size=200, random_state=0)
results = skm.cross_validate(clf, D, High, cv=validation)
print(f"驗證準確率: {results['test_score'][0]:.4f}")

驗證準確率 68.5%，比訓練準確率 72.75% 低約 4%，符合預期——訓練集上總是有輕微的過擬合。

成本複雜度剪枝 (Cost-Complexity Pruning)

現在我們進行完整的剪枝流程：先分割訓練/測試集，在訓練集上用 10-fold CV 遍歷 ccp_alpha 路徑，選出最佳 \(\alpha\) 後在測試集上評估。

# Block 5: Cost-Complexity Pruning 完整流程
import numpy as np
import pandas as pd
from ISLP import load_data, confusion_table
from ISLP.models import ModelSpec as MS
from sklearn.tree import DecisionTreeClassifier as DTC
from sklearn.metrics import accuracy_score
import sklearn.model_selection as skm

Carseats = load_data('Carseats')
High = np.where(Carseats.Sales > 8, "Yes", "No")
model = MS(Carseats.columns.drop('Sales'), intercept=False)
D = model.fit_transform(Carseats)
feature_names = list(D.columns)
X = np.asarray(D)

# 50/50 split
X_train, X_test, High_train, High_test = skm.train_test_split(
    X, High, test_size=0.5, random_state=0
)

# 未剪枝的完整樹
clf = DTC(criterion='entropy', random_state=0)
clf.fit(X_train, High_train)
full_acc = accuracy_score(High_test, clf.predict(X_test))
print(f"未剪枝樹測試準確率: {full_acc:.4f}")

# Cost-complexity pruning with CV
ccp_path = clf.cost_complexity_pruning_path(X_train, High_train)
kfold = skm.KFold(10, random_state=1, shuffle=True)
grid = skm.GridSearchCV(clf,
    {'ccp_alpha': ccp_path.ccp_alphas},
    refit=True, cv=kfold, scoring='accuracy')
grid.fit(X_train, High_train)

best_ = grid.best_estimator_
n_leaves = best_.tree_.n_leaves
print(f"最佳 CV 準確率: {grid.best_score_:.4f}")
print(f"剪枝後葉子數: {n_leaves}")

# 測試集評估
test_acc = accuracy_score(High_test, best_.predict(X_test))
print(f"剪枝後測試準確率: {test_acc:.4f}")

# 混淆矩陣
confusion = confusion_table(best_.predict(X_test), High_test)
print("混淆矩陣:")
print(confusion)

剪枝後剩 30 片葉子，測試準確率 72%（比完整樹的 73.5% 略低）。這說明在此資料集上，剪枝只砍掉 5 片葉子，效果有限——但 CV 選模本身是無偏的程序，不同 random seed 下結果會有變異，這正是模型選擇的固有隨機性。

8.3.2 迴歸樹實作

迴歸樹的流程與分類樹幾乎相同，只是改用 DecisionTreeRegressor，且評估指標從準確率改為 MSE。我們使用 Boston 房價資料集。

# Block 6: 迴歸樹擬合、剪枝與評估
import numpy as np
import pandas as pd
from ISLP import load_data
from ISLP.models import ModelSpec as MS
from sklearn.tree import DecisionTreeRegressor as DTR
import sklearn.model_selection as skm

Boston = load_data("Boston")
model = MS(Boston.columns.drop('medv'), intercept=False)
D = model.fit_transform(Boston)
feature_names = list(D.columns)
X = np.asarray(D)
y = Boston['medv']

# 70/30 split
X_train, X_test, y_train, y_test = skm.train_test_split(
    X, y, test_size=0.3, random_state=0
)

# 初始迴歸樹（max_depth=3）
reg = DTR(max_depth=3)
reg.fit(X_train, y_train)

# Cost-complexity pruning
ccp_path = reg.cost_complexity_pruning_path(X_train, y_train)
kfold = skm.KFold(5, shuffle=True, random_state=10)
grid = skm.GridSearchCV(reg,
    {'ccp_alpha': ccp_path.ccp_alphas},
    refit=True, cv=kfold, scoring='neg_mean_squared_error')
grid.fit(X_train, y_train)

best_ = grid.best_estimator_
test_mse = np.mean((y_test - best_.predict(X_test))**2)
print(f"測試 MSE: {test_mse:.2f}")
print(f"測試 RMSE: {np.sqrt(test_mse):.2f}  (約 $5300)")
print(f"最佳葉子數: {best_.tree_.n_leaves}")

迴歸樹的 RMSE 約 $5,300，意思是模型預測的中位數房價與真實值平均差距約五千多美元。關鍵分割變數包括 lstat（低收入人口比例）和 rm（房間數），與課本 §8.1 的理論分析一致。

8.3.3 Bagging 與隨機森林

Bagging 是隨機森林的特例（max_features = p，即每次分割考慮全部變數）。我們在 Boston 資料上比較兩者：

# Block 7: Bagging、隨機森林與變數重要性
import numpy as np
import pandas as pd
from ISLP import load_data
from ISLP.models import ModelSpec as MS
from sklearn.ensemble import RandomForestRegressor as RF
import sklearn.model_selection as skm

Boston = load_data("Boston")
model = MS(Boston.columns.drop('medv'), intercept=False)
D = model.fit_transform(Boston)
feature_names = list(D.columns)
X = np.asarray(D)
y = Boston['medv']

X_train, X_test, y_train, y_test = skm.train_test_split(
    X, y, test_size=0.3, random_state=0
)

# Bagging: max_features = p (=12)
bag_boston = RF(max_features=X_train.shape[1], n_estimators=500,
                random_state=0)
bag_boston.fit(X_train, y_train)
y_hat_bag = bag_boston.predict(X_test)
bag_mse = np.mean((y_test - y_hat_bag)**2)
print(f"Bagging (500 trees) 測試 MSE: {bag_mse:.2f}")

# Random Forest: max_features = 6
RF_boston = RF(max_features=6, random_state=0)
RF_boston.fit(X_train, y_train)
y_hat_RF = RF_boston.predict(X_test)
rf_mse = np.mean((y_test - y_hat_RF)**2)
print(f"Random Forest (max_features=6) 測試 MSE: {rf_mse:.2f}")

# 變數重要性
feature_imp = pd.DataFrame(
    {'importance': RF_boston.feature_importances_},
    index=feature_names
).sort_values(by='importance', ascending=False)
print("\n變數重要性（前 5）：")
print(feature_imp.head(5))

幾個關鍵觀察：

• Bagging MSE=14.6 約為單棵剪枝後迴歸樹（28.07）的 一半——集成效果顯著
• 隨機森林（max_features=6）的 MSE=20.0，在此資料上反而不如 Bagging——這提醒我們「RF 不一定總是贏」
• lstat 和 rm 合計佔了 70% 的重要性，其他變數貢獻都在 6% 以下
• 增加樹的數量從 100 到 500，MSE 變化極小（14.63→14.61）：Bagging/RF 不會因樹太多而過擬合

8.3.4 Boosting 實作

Boosting 與 Bagging 的關鍵差異：順序學習——每棵樹擬合前一棵的殘差。我們使用 GradientBoostingRegressor，並觀察訓練/測試誤差隨迭代次數的變化軌跡。

# Block 8: Gradient Boosting 與 staged_predict 軌跡
import numpy as np
import pandas as pd
from ISLP import load_data
from ISLP.models import ModelSpec as MS
from sklearn.ensemble import GradientBoostingRegressor as GBR
import sklearn.model_selection as skm

Boston = load_data("Boston")
model = MS(Boston.columns.drop('medv'), intercept=False)
D = model.fit_transform(Boston)
X = np.asarray(D)
y = Boston['medv']

X_train, X_test, y_train, y_test = skm.train_test_split(
    X, y, test_size=0.3, random_state=0
)

# Boosting: λ=0.001
boost_boston = GBR(n_estimators=5000, learning_rate=0.001,
                   max_depth=3, random_state=0)
boost_boston.fit(X_train, y_train)

# 測試誤差軌跡
test_error = np.zeros_like(boost_boston.train_score_)
for idx, y_ in enumerate(boost_boston.staged_predict(X_test)):
    test_error[idx] = np.mean((y_test - y_)**2)

print(f"最終訓練誤差: {boost_boston.train_score_[-1]:.2f}")
print(f"最終測試誤差 (λ=0.001): {test_error[-1]:.2f}")

# 比較 λ=0.2
boost2 = GBR(n_estimators=5000, learning_rate=0.2,
             max_depth=3, random_state=0)
boost2.fit(X_train, y_train)
y_hat2 = boost2.predict(X_test)
mse2 = np.mean((y_test - y_hat2)**2)
print(f"測試 MSE (λ=0.2): {mse2:.2f}")

# 預測
y_hat_boost = boost_boston.predict(X_test)
mse_boost = np.mean((y_test - y_hat_boost)**2)
print(f"測試 MSE (λ=0.001): {mse_boost:.2f}")

Boosting 的測試 MSE=14.48，與 Bagging（14.61）相當。λ=0.001 和 λ=0.2 的結果幾乎一樣——這並不常見。一般來說，較小的 λ（較慢的學習）需要更多棵樹，但往往能達到更好的泛化；較大的 λ 收斂快但有過擬合風險。這個例子中兩者打成平手，提醒我們調參效果因資料集而異。

8.3.5 貝氏加法迴歸樹 (BART)

BART 是 §8.2 介紹的貝氏集成方法，使用 ISLP.bart.BART 實作。與隨機森林/Boosting 不同，BART 透過 MCMC 後驗抽樣來產生預測分佈，能自然地給出不確定性估計。

# Block 9: BART 擬合與預測
import numpy as np
import pandas as pd
from ISLP import load_data
from ISLP.models import ModelSpec as MS
from ISLP.bart import BART
import sklearn.model_selection as skm

Boston = load_data("Boston")
model = MS(Boston.columns.drop('medv'), intercept=False)
D = model.fit_transform(Boston)
X_raw = D  # 保留 DataFrame
y_raw = Boston['medv']

# BART 要求 float64 numpy array
X = X_raw.to_numpy(dtype=np.float64)
y = y_raw.to_numpy(dtype=np.float64)

X_train, X_test, y_train, y_test = skm.train_test_split(
    X, y, test_size=0.3, random_state=0
)

# BART 建構參數：num_trees + burnin + ndraw（非 sklearn 慣例！）
# 使用小 burnin/ndraw 加速示範（正式使用建議 burnin=100, ndraw=10+）
bart_boston = BART(num_trees=200, burnin=5, ndraw=15, random_state=0)
bart_boston.fit(X_train, y_train)

yhat_test = bart_boston.predict(X_test)
bart_mse = np.mean((y_test - yhat_test)**2)
print(f"BART 測試 MSE: {bart_mse:.2f}")

# 變數 inclusion 比例
var_inclusion = pd.Series(
    bart_boston.variable_inclusion_.mean(0),
    index=D.columns
).sort_values(ascending=False)
print("\n變數 inclusion 比例（前 5）：")
print(var_inclusion.head(5))

📊 五種樹狀方法比較