① 依序閱讀:每節先讀說明;遇到互動元件時,先預測結果,再操作驗證。 ② 對照來源:章節旁標示課本章節與講義主題,正文保留公式、條件與算例;計算細節、完整實作與證明預設收合,可按需要展開。來源標記可查原教材。 ③ 檢查理解:答錯時先看回饋,再回到相關說明;標為延伸的內容可留待第二輪。 ④ 最後翻關鍵詞彙卡自測術語,並用 重點速查與來源查閱。
教科書:An Introduction to Statistical Learning with Applications in Python(ISLP);
進階參考:The Elements of Statistical Learning(ESL)。章節旁的來源標記可點到本頁書目。
第 3 章的最小平方法很好用:沒有調整參數、有封閉解、係數可以直接解讀。講義強調線性模型容易解讀,預測表現也常常不錯,所以本章不放棄線性模型,而是把「最小平方擬合」換成其他擬合方法。這在 $n\approx p$ 或 $p>n$ 時特別重要。
講義列出四個考慮替代最小平方法的理由:
三類方法都透過接受一些偏差來降低變異。最小平方是無偏的, 但在 $p$ 大的時候,高變異會降低預測表現;只要願意接受一點偏差,測試誤差往往降得很明顯。 本章會比較這些方法如何取捨偏差與變異。
$$\text{RSS} = \sum_{i=1}^{n}\left(y_i - \beta_0 - \sum_{j=1}^{p}\beta_j x_{ij}\right)^2$$下表整理三類方法的差別,讀完整章後可回來比較。
| 留幾個變數 | 係數會剛好是 0 嗎 | 有調整參數嗎 | 座標系 | 主要弱點 | |
|---|---|---|---|---|---|
| 最小平方 | 全部 $p$ 個 | 幾乎不會 | 沒有 | 原始變數 | $p$ 大時變異很大 |
| 子集選擇 | $k$ 個(自己選) | 沒被選的就是 0 | $k$ | 原始變數 | 搜尋空間太大、不穩定 |
| Ridge(L2) | 全部 $p$ 個 | 不會 | $\lambda$ | 原始變數 | 不做變數選擇 |
| Lasso(L1) | 由 $\lambda$ 決定 | 會 | $\lambda$ | 原始變數 | 相關變數之間的選擇很不穩 |
| PCR / PLS | 全部 $p$ 個都參與 | 不會 | $M$ | 換過的座標 | 新座標不好解讀 |
為什麼 $p$ 逼近 $n$ 的時候,「把全部變數丟進去做最小平方」是個壞主意?
最直白的想法是把每一種變數組合都試一次,挑最好的。這就是最佳子集選擇(best subset selection)。$p$ 個變數有 $2^p$ 種組合(含空模型)。講義的演算法分三步:
可以。直接用交叉驗證比較全部 $2^p$ 個子集,在統計上合理,也可能找到比兩階段法更好的預測模型。教科書仍介紹兩階段,主要有兩個原因。
以 $p=10$ 為例,共有 $2^{10}=1024$ 個候選模型。
| 做法 | 第一階段 | 第二階段 | 總擬合次數(約) |
|---|---|---|---|
| 兩階段+$C_p$、AIC、BIC | 1024 個模型各擬合一次、算訓練 RSS | 在 11 個模型間加上懲罰項比較 | $1024$ |
| 直接對所有子集做 5 折 CV | 不需要 | 1024 個模型各做 5 折 | $1024\times5=5120$ |
兩階段搭配 $C_p$、AIC、BIC 時,可以避免對每個候選模型反覆做 CV,省下的計算最多。若第二階段用 CV,而且照規矩做(見下方),第一階段的搜尋也要在每一折的訓練部分重新執行,這時兩階段就省不了多少計算。
假設 $p=3$,只看單變數模型:
| 模型 | 訓練 RSS | CV 誤差 |
|---|---|---|
| $Y\sim X_1$ | 90 | 120 |
| $Y\sim X_2$ | 95 | 105 |
| $Y\sim X_3$ | 110 | 130 |
兩階段法的第一階段用訓練 RSS 選出 $X_1$(90 最小),第二階段只會看到這個單變數模型,不會再考慮 $X_2$。直接對所有模型做 CV,則會選到 CV 誤差較小的 $X_2$。也就是說,兩者一般不保證相同:
$$\underset{S:\,|S|=k}{\arg\min}\;\mathrm{RSS}(S)\;\ne\;\underset{S:\,|S|=k}{\arg\min}\;\mathrm{CV}(S)\quad\text{(一般情況)}.$$第一階段回答「每個大小下,哪些變數最能擬合訓練資料」:$M_k=\arg\min_{|S|=k}\mathrm{RSS}(S)$。第二階段回答「模型應該多複雜」:$k^*=\arg\min_k\text{準則}(M_k)$。這是方便、有系統的搜尋策略,方便搭配 AIC、BIC,但不是數學上非如此不可的規定。
講義接著指出,這套做法不限最小平方,邏輯斯迴歸等其他模型也能用;這時以離差(deviance)取代 RSS 當擬合量。
離差和負二倍最大對數概似成比例:$D=-2\log\hat L+\text{常數}$,常數只和資料有關、與模型無關。愈小代表擬合愈好,角色和 RSS 一樣。
對誤差變異數 $\sigma^2$ 已知的常態線性模型,$-2\log\hat L=\mathrm{RSS}/\sigma^2+\text{常數}$,離差就是 RSS 換了量尺。所以講義的兩階段法可以原封不動搬過去:第一階段在同大小模型間比離差;第二階段同樣不能用訓練離差(巢狀模型加變數不會讓它變大),要改用 AIC、BIC 或 CV。這時 AIC 寫成 $-2\log\hat L+2k$、BIC 寫成 $-2\log\hat L+k\log n$,$k$ 計入實際估計的參數;軟體的數值量尺見下一節的分類離差的數值。
講義指出最佳子集的兩個問題。計算上,$2^p$ 長得太快,$p$ 超過 40 左右就算不完。統計上,搜尋空間愈大,愈容易找到在訓練資料上看起來很好、對未來資料卻沒有預測力的模型;龐大的搜尋空間會造成過度擬合,也讓係數估計的變異變大。逐步法只探索小得多的模型集合,所以是有吸引力的替代方案。
可以,但要看第二階段用的是什麼。這個問題常被叫做 data snooping 或 selection-induced overfitting(因大量搜尋而產生的過度擬合),和刻意搜尋顯著 $p$ 值的 p-hacking 不是同一件事。
假設有 $p=100$ 個變數,全部都和 $Y$ 無關。第一階段在大量組合中搜尋,可能碰巧找到訓練 $R^2$ 很高(例如 0.8)的模型,那只是擬合到了雜訊。如果第二階段能可靠地估測試誤差,就會發現它的預測能力很差,甚至不如空模型 $M_0$,於是把它淘汰。
| 第二階段的方法 | 能否處理第一階段的過度擬合? |
|---|---|
| 母體的期望測試誤差(理想,實際上不知道) | 可以辨識每個候選模型的真實預測風險 |
| 獨立驗證集,或正確執行的 CV | 可以緩解,但仍有估計誤差,挑最小值也有選擇偏差 |
| $C_p$、AIC、BIC | 可能修正不足 |
$C_p$、AIC、BIC 的懲罰只看模型大小 $k$,沒有反映「從 $\binom{p}{k}$ 個候選中挑出最好的一個」所增加的選擇複雜度;它們的推導前提是模型事先固定。正確的 CV 要在每個訓練折內重新做最佳子集,見上方「為什麼要分兩階段」收合區。
這時第一階段的搜尋不會直接影響第 3 步的測試誤差。但第 2 步挑出的最小驗證誤差仍偏樂觀;若還想對選中變數的係數做顯著性檢定,也不能忽略前面的選擇過程(選後推論的問題)。
最佳子集在高維時容易對訓練資料過度擬合,透過適當的樣本外驗證可以大幅緩解;兩階段法本身不保證完全消除選擇偏差。
前向逐步選擇(forward stepwise selection)從空模型開始,每一步加入一個「讓訓練 RSS 降最多」的變數,直到所有變數都在模型裡。 後向逐步選擇(backward stepwise selection)反過來,從含全部 $p$ 個變數的模型開始,每次刪掉最沒用的一個。 兩者都只要擬合
$$1 + \sum_{k=0}^{p-1} (p-k) = 1 + \frac{p(p+1)}{2} \; \text{個模型}$$$p = 20$ 時這是 211 個,而最佳子集是 $2^{20}=$ 1,048,576 個。 代價是貪婪走法不保證找到最佳子集(訓練誤差最低的子集)。講義的例子:$p=3$ 時,最佳單變數模型含 $X_1$,最佳雙變數模型卻是 $\{X_2,X_3\}$;前向法的 $M_1$ 含 $X_1$,$M_2$ 也只能含 $X_1$ 再加一個,於是錯過最佳的雙變數模型。下面第二張圖用四變數示意把這件事畫出來。
講義提醒:前向逐步雖然只擬合 $1+p(p+1)/2$ 個模型,卻是由資料引導的搜尋(guided search),實際涵蓋的有效模型空間遠比這個數字大。
以 $p=4$ 為例,變數為 $X_1,X_2,X_3,X_4$。
加上空模型,一共擬合 $1+4+3+2+1=11$ 個模型,而全部可能的子集有 $2^4=16$ 個。
這 11 個模型由訓練資料決定,並非事先固定。換一批資料,可能先選 $X_1$ 再選 $X_3$,也可能先選 $X_4$ 再選 $X_2$。$p=4$ 時變數加入的順序共有 $4!=24$ 種;一次執行只沿其中一條路徑前進,但走哪條路本身就用了 $y$。每個子集都是某種加入順序的前段,所以在不同資料下,全部 16 個子集都可能被選中。
| 數什麼 | $p=4$ 的數量 |
|---|---|
| 一次執行實際擬合的模型 | 11 |
| 不同資料下可能被選中的模型 | 全部 16 個子集 |
| 可能的加入順序 | $4!=24$ |
因此,統計上的選擇複雜度不能只用 11 個模型衡量,這就是講義說的 effective model space。用逐步法挑出的模型,同樣需要上一段所說的樣本外驗證。
在最小平方下,如果特徵在樣本中彼此正交,前向逐步選擇和最佳子集選擇對每個大小 $k$ 都會選出相同的變數組合(忽略並列的情況)。條件是樣本中的正交,不只是母體中不相關。
模型含截距,特徵已中心化,而且 $X_i^\top X_j=0\ (i\ne j)$。令 $y_c$ 為中心化後的反應,對任意變數子集 $S$,
$$\mathrm{RSS}(S)=\mathrm{TSS}-\sum_{j\in S}\frac{(X_j^\top y_c)^2}{X_j^\top X_j}.$$右邊是每個特徵的貢獻直接相加:一個變數能讓 RSS 降多少,和其他變數有沒有在模型裡無關。要選 $k$ 個變數,最佳子集只要挑 $(X_j^\top y_c)^2/(X_j^\top X_j)$ 最大的 $k$ 個;前向逐步每一步挑的也正是剩下的變數中貢獻最大的那一個。兩者自然相同。
| 特徵 | 單獨加入時 RSS 的減少量 |
|---|---|
| $X_1$ | 100 |
| $X_2$ | 60 |
| $X_3$ | 30 |
最佳子集:$k=1$ 選 $X_1$,$k=2$ 選 $\{X_1,X_2\}$,$k=3$ 選全部。前向逐步:第一步選 $X_1$(RSS 減 100),第二步選 $X_2$(再減 60),第三步選 $X_3$(再減 30)。兩者完全相同。
若 $X_1$ 與 $X_2$ 高度相關,單獨用 $X_1$ 能大幅降低 RSS,但加入 $X_1$ 後,$X_2$ 能提供的額外資訊就很少;反而 $\{X_2,X_3\}$ 可能比 $\{X_1,X_3\}$ 更好。某個變數單獨表現很好,不代表它屬於最佳的多變數組合。前向逐步不能移除已選的變數,可能錯過這個組合;最佳子集則會搜尋到。
| 特徵的條件 | 兩種方法一定相同? |
|---|---|
| 樣本中完全正交 | 是(忽略並列) |
| 母體中不相關,但樣本有相關 | 不一定 |
| 樣本相關性很小 | 常常接近,但不保證 |
| 特徵明顯相關 | 不保證 |
即使母體滿足 $\operatorname{Cov}(X_1,X_2)=0$,有限樣本的 $\widehat{\operatorname{Cov}}(X_1,X_2)$ 通常也不是 0,所以母體不相關不能保證兩種演算法在訓練資料上給出相同結果。前向逐步與最佳子集的差異,來自變數之間的相依結構:單一變數的貢獻會隨已選的變數而改變。
讀逐步加入變數時,每個候選模型包含前一步已選的變數:CRBI → 加 Hits → 加 PutOuts…一路只加不減。這正是前向逐步的定義,也正是它可能錯過最佳子集的原因。
Ch06-varselect-lab-zh.ipynb
np.int64(7)
bic_b 是照 models2 的索引順序(1、2、…、19)收集的,0-based 的 argmin() = 7 指的是第 8 個元素,也就是 8 變數模型;而前向那邊(下一節)是 6 變數。同一份資料、同一個準則,兩個搜尋方向給出不同大小的模型——講義講的就是這件事。
Ch06-varselect-lab-zh.ipynb
| 搜尋法 | 要擬合幾個模型 | $p=20$ 時 | $n<p$ 能用嗎 | 保證找到最佳子集嗎 |
|---|---|---|---|---|
| 最佳子集 | $2^p$ | 1,048,576 | 可以(只比較可估計的小子集) | 會 |
| 前向逐步 | $1+p(p+1)/2$ | 211 | 可以 | 不保證 |
| 後向逐步 | $1+p(p+1)/2$ | 211 | 不能(要先擬合全模型) | 不保證 |
| 混合逐步 | 略多於逐步 | 約 211+ | 看實作 | 不保證(但比較接近) |
同一份資料上做最佳子集、前向逐步、後向逐步。大小同樣是 $k$ 的三個模型裡,哪一個的訓練 RSS 最小?
空模型只含截距,對每筆資料都預測同一個常數;即使沒有預測變數,仍然可以做預測。比較加入變數後的模型時,就以這個常數模型為基準。對平方損失,這個常數是訓練反應的平均值;預測新資料時也要沿用訓練平均,不能先看驗證資料的反應值,再用它們算平均。
對任意常數 $c$,先寫成 $y_i-c=(y_i-\bar y)+(\bar y-c)$,再展開平方。因為 $\sum_i(y_i-\bar y)=0$,交叉項消失:
$$\sum_i(y_i-c)^2=\sum_i(y_i-\bar y)^2+n(c-\bar y)^2.$$第二項非負,且只有 $c=\bar y$ 時為零;所以最小 RSS 等於 TSS。這也是含截距模型的 $R^2$ 使用 TSS 當基準的原因。
講義的 Credit 有十個原始預測變數,但三類別的變數需要兩個虛擬欄,因此迴歸設計矩陣可有十一個斜率欄。若把虛擬欄分別搜尋,「一個變數」可能指一個欄,而不是整個類別因素。分析前要決定,是逐欄選擇,還是讓整個因素一起進出模型?這會改變候選模型的範圍,也影響選出結果如何解讀。下方 Group Lasso 就是以組為單位選擇的做法。
Credit 圖的每一個點是一個候選模型,紅色前緣只表示同樣大小下的訓練贏家。自由度相同的模型可用 RSS 排序;模型大小不同時,就不能只依 RSS 選出最終模型。講義的 Hitters Lab 則用 Salary 當反應,先移除缺失觀測,再透過設計矩陣處理類別欄。不要把 Credit 的圖上模型大小,當作 Hitters 的結果。
前向法錯過 $\{X_2,X_3\}$ 的例子見正文;Backward 的限制方向相反:一旦刪掉某欄,單純後向法不會把它加回。混合逐步法允許加入之後再考慮刪除,但仍不是窮舉所有子集。逐步法的搜尋路徑由 $y$ 決定,選擇偏差的討論見上方「有效模型空間」與「第二階段的驗證」兩個收合區。
高維時,Forward 只在仍可識別的子模型中前進;含截距需要設計矩陣滿欄秩,$n$ 大於欄數才有正殘差自由度。Backward 先需要完整模型可估計;即使軟體用廣義逆算出一組係數,也要另外檢查係數是否唯一,以及常規推論的條件是否成立。
Lab 的自訂迴圈沒有列入空模型,因此搜尋範圍和講義的 $M_0,\ldots,M_p$ 不完全相同。ISLP 的完整路徑則包含空模型。讀 argmin 時,要先核對索引是「零起算位置」還是「變數個數」,不能直接把兩者畫上等號。
來源入口:Credit 資料字典、Hitters 資料字典、ISLP 逐步選擇、sklearn SequentialFeatureSelector。後者以交叉驗證分數逐步比較,不能未檢查 scoring 就宣稱重現 Lab 的 RSS 路徑。
講義延伸來源:math.stackexchange.com。
講義延伸來源:quantifyinghealth.com。
第二階段的問題是:怎麼比較「大小不同」的模型?訓練 RSS 與 $R^2$ 不能用, 因為它們單調偏好大模型。兩條路——間接(把訓練誤差加一個懲罰項,修掉它的樂觀偏差) 與直接(用驗證集或交叉驗證真的去估測試誤差,下一節講)。
常見的間接評估準則有四個。設模型有 $d$ 個預測變數,$\hat\sigma^2$ 是用全模型估出來的誤差變異數:
$$C_p = \frac{1}{n}\left(\mathrm{RSS} + 2 d \hat\sigma^2\right), \qquad \mathrm{BIC} = \frac{1}{n}\left(\mathrm{RSS} + \log(n)\, d \hat\sigma^2\right)$$ $$\text{Adjusted } R^2 = 1 - \frac{\mathrm{RSS}/(n-d-1)}{\mathrm{TSS}/(n-1)}$$兩個公式長得幾乎一樣,只差在每多一個變數要罰多少:$C_p$ 罰 $2\hat\sigma^2$,
BIC 罰 $\log(n)\,\hat\sigma^2$。整數樣本數 $n\ge8$ 時 $\log n>2$,所以BIC 罰得比 $C_p$ 重,偏好更小的模型。
$n$ 愈大差距愈明顯:Credit 的 $n = 400$,$\log 400 \approx 5.99$,BIC 的懲罰約是 $C_p$ 的三倍。
$C_p$ 的兩個要點:所有候選模型共用同一個 $\hat\sigma^2$,通常由含全部 $p$ 個預測變數的全模型估計;$2d\hat\sigma^2$ 是在修正訓練 RSS 對測試誤差的樂觀偏差。加一個變數,訓練 RSS 必須降超過 $2\hat\sigma^2$,$C_p$ 才會變小。
我們假設 $Y=f(X)+\varepsilon$、$\operatorname{Var}(\varepsilon)=\sigma^2$。$\sigma^2$ 是資料生成過程的雜訊變異數,不會因為我們換了模型而改變,所以所有候選模型用同一個估計值:
$$\hat\sigma^2=\frac{\mathrm{RSS}_{\text{full}}}{n-p-1}.$$如果每個候選模型都用自己的 RSS 估 $\sigma^2$,較小的模型因為欠擬合,會把模型偏差也算成雜訊,估計值偏大。用全模型當共同基準,前提是全模型足以描述真實的平均反應函數;否則 $\hat\sigma^2$ 也會有偏差。
訓練 RSS 因為擬合到雜訊,期望比 $n\sigma^2$ 少 $d\sigma^2$;拿同一組係數去預測新反應,參數估計的不確定性又讓測試 RSS 多出 $d\sigma^2$。兩半合起來是 $2d\sigma^2$。完整推導(含模型有近似偏差的情形),以及 AIC、BIC 的懲罰從哪裡來,見本節後面的三個準則的懲罰項推導。
$d$ 變成 $d+1$,懲罰項增加 $2\hat\sigma^2$。只有當
$$\mathrm{RSS}_{\text{old}}-\mathrm{RSS}_{\text{new}}>2\hat\sigma^2$$時 $C_p$ 才會下降。例如 $\hat\sigma^2=4$,新變數必須讓 RSS 降超過 8。假設某個新變數讓 RSS 降了 12:
| $\hat\sigma^2$ | 新增一個變數的懲罰 | RSS 降 12,值得加入嗎? |
|---|---|---|
| 1 | 2 | 值得 |
| 4 | 8 | 值得 |
| 9 | 18 | 不值得 |
| 16 | 32 | 不值得 |
直覺是:雜訊愈大,即使新變數完全沒用,也可能因為剛好擬合到隨機波動而讓 RSS 明顯下降,所以要求更大的改善才相信它。在誤差獨立、同變異的常態線性模型下,考慮一個事先指定、不在原模型欄空間內的新變數;若原模型已經正確、新變數的真實係數為 0,加入它造成的下降滿足
$$\frac{\mathrm{RSS}_{\text{old}}-\mathrm{RSS}_{\text{new}}}{\sigma^2}\sim\chi^2_1,$$它本身就隨 $\sigma^2$ 放大(分母是母體的 $\sigma^2$;從許多候選中挑出 RSS 降最多的那個變數,下降量就不再服從這個分布)。所以 $C_p$ 的門檻其實是在看相對於雜訊的改善:$\Delta\mathrm{RSS}/\sigma^2>2$。這也表示「$\sigma^2$ 愈大,一定選愈少變數」並不嚴格成立,因為 RSS 的下降量本身也會隨雜訊改變。
第一階段固定 $d$,懲罰項對同大小的模型都相同,最小化 $C_p$ 等同於最小化訓練 RSS。第二階段比較不同的 $d$,$2d\hat\sigma^2$ 才發揮作用,避免只因變數變多、訓練 RSS 變小就一直加變數。
statsmodels 的 .aic 用的是 log-likelihood 版,數值與 $C_p$ 不同,所以下圖畫成獨立一條線;為什麼成正比見下方「懲罰項各從哪裡來」收合區,軟體量尺見「先分清估計量與比較量尺」收合區。
| 準則 | 式子(最小平方) | 每多一個變數罰多少 | $n=400$ 時的懲罰 | Credit 選出 |
|---|---|---|---|---|
| $C_p$ | $(\mathrm{RSS}+2d\hat\sigma^2)/n$ | $2\hat\sigma^2$ | $2\hat\sigma^2$ | 6 個 |
| AIC(課本) | $\mathrm{RSS}+2d\hat\sigma^2$ | $2\hat\sigma^2$ | 與 $C_p$ 成正比 | 同 $C_p$ |
| BIC | $(\mathrm{RSS}+\log(n)\,d\hat\sigma^2)/n$ | $\log(n)\,\hat\sigma^2$ | $5.99\hat\sigma^2$ | 4 個 |
| 調整後 $R^2$ | $1-[\mathrm{RSS}/(n-d-1)]/[\mathrm{TSS}/(n-1)]$ | 分母的 $n-d-1$ 變小 | 很輕 | 7 個 |
| 交叉驗證 | 直接估測試誤差 | 不需要 $\hat\sigma^2$ | — | 6 個 |
AIC 與 BIC 的目標不同:AIC 著重預測準確度(漸近效率),BIC 著重找出真實模型(選模一致性)。真實模型沒有模型偏差,但在有限樣本下,它的結構不一定最容易估計,也不一定預測得最好。
假設候選模型為 $M_1,\ldots,M_K$,各有真實的預測風險 $R(M_k)$。若知道真實分布,就能挑出風險最小的 $M_{\text{oracle}}=\arg\min_kR(M_k)$。實際上只能用 AIC 選出 $\hat M_{\text{AIC}}$。漸近效率的典型意思是
$$\frac{R(\hat M_{\text{AIC}})}{\min_kR(M_k)}\xrightarrow{P}1\qquad(n\to\infty),$$也就是資料愈多,AIC 所選模型的風險愈接近候選中最好的那個。這裡的風險要適當定義(例如 KL 散度或超額預測風險),而且這是特定理論條件下的性質,不是所有選模問題都成立。AIC 設計的目標就是近似模型對新資料的期望 KL 差異,所以它偏向預測表現。
| AIC | BIC | |
|---|---|---|
| 主要目標 | 預測能力 | 辨識真實模型 |
| 漸近性質(特定條件下) | 漸近效率 | 選模一致性 |
| 對微弱訊號 | 較願意保留 | 較容易排除 |
| 理論優勢 | 真實函數無法用有限維模型完整表示時 | 真實模型固定、有限維且在候選清單中時 |
例如真正的模型只含 $X_1$,候選卻還有 $X_2,X_3$。即使樣本數愈來愈多,AIC 仍有正的機率選到多餘變數,所以它沒有一般意義下的選模一致性;BIC 在真實模型屬於候選清單等條件下,選中它的機率會趨近 1。反過來,若真正的函數非常複雜(例如無限多項的多項式),沒有任何有限維候選能完全表示它,選模的目標就變成找偏差–變異取捨好的近似模型,這時 AIC 的漸近效率比較有意義。
假設真實模型為 $Y=2X_1+0.05X_2+\varepsilon$,$\varepsilon\sim N(0,1)$,$n=100$,$X_1,X_2$ 獨立且已標準化。$X_2$ 的係數 0.05 不是零,所以真實模型包含 $X_2$。比較 $X_2$ 對預測風險的額外貢獻(其他部分相同):
| 包含 $X_2$(真實結構) | 不包含 $X_2$ | |
|---|---|---|
| 偏差 | 0 | 有偏差 |
| 額外的變異 | 約 $\sigma^2/n=0.01$ | 0 |
| 額外的平方偏差 | 0 | $0.05^2=0.0025$ |
$0.0025<0.01$:不放 $X_2$ 反而可能有較小的測試誤差。估計這個微弱訊號付出的變異,比忽略它造成的偏差還大。
當 $n\to\infty$,$\operatorname{Var}(\hat\beta_2)\approx1/n\to0$,估計 $X_2$ 的成本消失,忽略它的偏差卻一直存在。所以真實模型固定、有限維且正確指定時,樣本夠大之後選到真實模型才會佔優勢,這也是 BIC 選模一致性的價值。
要分清楚:若連真正的參數值都知道,真實模型在平方誤差下當然預測最好。這裡說的是「結構選對了,參數仍要用有限樣本估計」。預測誤差同時取決於偏差與變異;有限樣本下,略有偏差但變異小得多的簡單模型,可能預測得更好。AIC 的漸近效率也不代表在每個資料集上 AIC 的測試誤差都比 BIC 小,兩者的理論優勢建立在不同的假設與目標上。
三個準則都是「訓練擬合+隨模型變大的懲罰」,公式很像,懲罰的理由卻不同:
共同的起點是:加變數不會讓訓練 RSS 變大,也不會讓最大對數概似變小,但這不代表預測變好。以下依序推導三個懲罰。
記號。$n$ 是訓練觀測數;$y=(y_1,\ldots,y_n)^\top$ 是訓練資料的反應向量;$X$ 是候選模型的 $n\times r$ 設計矩陣(每列一筆觀測,含截距欄)。$d$ 是候選模型的預測變數個數,$r$ 是迴歸係數個數(含截距時 $r=d+1$),$k$ 是所有實際估計的參數個數(例如再加上誤差變異數)。$\mathrm{RSS}=\|y-\hat y\|^2$ 是訓練殘差平方和,$\|\cdot\|$ 是向量長度,$E[\cdot]$ 是對資料的隨機性取期望。
先固定 $X$ 與候選模型,不用 $y$ 搜尋欄位。設 $y=\mu+\varepsilon$:$\mu=E[y]$ 是真實的平均反應向量(不一定落在 $X$ 的欄空間內,也就是模型可以有偏差),$\varepsilon$ 是誤差向量,$E[\varepsilon]=0$、$\operatorname{Cov}(\varepsilon)=\sigma^2I$,$\sigma^2$ 是雜訊變異數,$I$ 是 $n\times n$ 單位矩陣。OLS 擬合值 $\hat y=Hy$,其中帽子矩陣 $H=X(X^\top X)^{-1}X^\top$ 把 $y$ 投影到 $X$ 的欄空間;它對稱冪等($H^\top=H$、$H^2=H$),跡 $\operatorname{tr}(H)$(對角元素和)等於 $r$。
訓練 RSS。$y-\hat y=(I-H)\mu+(I-H)\varepsilon$,交叉項期望為 0,且 $E\|(I-H)\varepsilon\|^2=\sigma^2\operatorname{tr}(I-H)$:
$$E[\mathrm{RSS}]=\|(I-H)\mu\|^2+(n-r)\sigma^2.$$第一項是近似偏差,第二項是雜訊。擬合時有 $r\sigma^2$ 的雜訊被模型「吃掉」,所以比 $n\sigma^2$ 少。
新反應。在相同 $X$ 上取一組新的反應 $y^{\text{new}}=\mu+\varepsilon^{\text{new}}$,$\varepsilon^{\text{new}}$ 與 $\varepsilon$ 獨立、分布相同;測試 RSS 定義為 $\|y^{\text{new}}-\hat y\|^2$,也就是用訓練得到的 $\hat y$ 去預測新反應的誤差。$y^{\text{new}}-\hat y=(I-H)\mu+\varepsilon^{\text{new}}-H\varepsilon$,三項之間的交叉項期望都是 0:
$$E\|y^{\text{new}}-\hat y\|^2=\|(I-H)\mu\|^2+n\sigma^2+r\sigma^2.$$多出的 $r\sigma^2$ 來自估計參數的變異。兩式相減,偏差項正好抵消:
$$E[\text{測試 RSS}]=E[\mathrm{RSS}]+2r\sigma^2.$$每多估一個參數,訓練 RSS 相對於預測誤差的樂觀偏差就多 $2\sigma^2$:一半來自擬合雜訊,一半來自估計的不確定性。用所有模型共用的 $\hat\sigma^2$ 取代 $\sigma^2$、再除以 $n$,就是 $C_p=(\mathrm{RSS}+2d\hat\sigma^2)/n$(截距的 $2\hat\sigma^2/n$ 各模型相同,省略與否見$C_p$ 的風險目標與常數)。$\hat\sigma^2$ 不偏就能在期望中取代 $\sigma^2$,不需要和 RSS 獨立。若模型本身是用同一份 $y$ 搜尋出來的,$H$ 已經依賴 $y$,這個推導就不再成立。
設 $g$ 是實際產生資料的未知機率密度(真實分布)。候選模型是一族密度 $f(y\mid\theta)$:$\theta$ 是模型的參數向量,$f(y\mid\theta)$ 是參數取 $\theta$ 時模型給資料 $y$ 的機率密度,簡寫為 $f_\theta$。令 $Y$ 是從 $g$ 抽出的一筆新觀測,$E_g$ 表示對 $Y\sim g$ 取期望。兩個分布的差距用 KL 散度衡量:
$$D_{\mathrm{KL}}(g\,\|\,f_\theta)=E_g\!\left[\log\frac{g(Y)}{f(Y\mid\theta)}\right]=E_g[\log g(Y)]-E_g[\log f(Y\mid\theta)].$$第一項和候選模型無關,所以比較模型只要看第二項:期望對數概似愈大,離真實分布愈近。$g$ 未知,參數只能用訓練資料估。記 $\ell(\theta)=\log f(y\mid\theta)$ 為訓練資料的對數概似,$\hat\theta=\arg\max_\theta\ell(\theta)$ 為最大概似估計。我們想知道 $\hat\theta$ 在另一組同樣大小、獨立的新資料 $y^{\text{new}}$ 上的對數概似 $\ell_{\text{new}}(\hat\theta)=\log f(y^{\text{new}}\mid\hat\theta)$,手上卻只有 $\ell(\hat\theta)$。$\hat\theta$ 本來就是讓 $\ell$ 最大的值,所以 $\ell(\hat\theta)$ 偏高。
樂觀偏差約為 $k$。假設模型正確指定、參數個數 $k$ 固定,並滿足一般的正則條件。令 $\theta_0$ 為真實參數(模型正確時 $g=f_{\theta_0}$),$I(\theta_0)=-E_g\bigl[\partial^2\log f(Y\mid\theta)/\partial\theta\,\partial\theta^\top\bigr]_{\theta=\theta_0}$ 為單筆觀測的費雪資訊矩陣(對數密度曲率的期望)。記 $Q=n(\hat\theta-\theta_0)^\top I(\theta_0)(\hat\theta-\theta_0)$,它衡量 $\hat\theta$ 離 $\theta_0$ 多遠;大樣本下 $Q$ 近似自由度 $k$ 的卡方分布 $\chi^2_k$,所以 $E[Q]\approx k$。在兩個最大值附近各做二階展開:
因為 $E[\ell(\theta_0)]=E[\ell_{\text{new}}(\theta_0)]$,兩式相減得
$$E[\ell(\hat\theta)]-E[\ell_{\text{new}}(\hat\theta)]\approx k.$$結構和 $C_p$ 一樣:一半是擬合到訓練雜訊,一半是參數估計的不確定性。所以 $\ell(\hat\theta)-k$ 才是新資料對數概似的近似不偏估計,乘上 $-2$:
$$\mathrm{AIC}=-2\ell(\hat\theta)+2k.$$模型指定錯誤時,$k$ 要換成一個跡項 $\operatorname{tr}(J^{-1}K)$(令 $\theta^\ast$ 為使 KL 散度最小的參數值,$J=-E_g[\nabla^2\log f(Y\mid\theta^\ast)]$ 是單筆對數密度的負期望二階導數矩陣,$K=\operatorname{Var}_g[\nabla\log f(Y\mid\theta^\ast)]$ 是一階導數(分數函數)的變異數矩陣;這個版本稱為 TIC);只有模型正確時 $J=K$,跡才化簡成 $k$。
高斯線性迴歸的寫法。$x_i^\top$ 是 $X$ 的第 $i$ 列,$\beta$ 是 $r$ 維係數向量,$\theta=(\beta,\sigma^2)$。各觀測獨立且 $y_i\sim N(x_i^\top\beta,\sigma^2)$ 時,
$$\ell(\beta,\sigma^2)=-\frac n2\log(2\pi\sigma^2)-\frac{\mathrm{RSS}}{2\sigma^2}.$$$\sigma^2$ 也要估時,最大概似估計是 $\hat\sigma^2=\mathrm{RSS}/n$,代回得 $-2\ell(\hat\theta)=n\log(\mathrm{RSS}/n)+n[1+\log(2\pi)]$。參數有 $d$ 個斜率、截距與 $\sigma^2$,$k=d+2$;後兩個每個模型都有,比較時可以略去:
$$\mathrm{AIC}=n\log(\mathrm{RSS}/n)+2d+\text{常數}.$$statsmodels 的 .aic、.bic 用的就是這種對數概似量尺;常數與參數個數的算法見先分清估計量與比較量尺。
若 $\sigma^2$ 已知,或所有模型共用同一個固定值,只有 $r$ 個迴歸係數要估,$k=r$,而且
$$-2\ell(\hat\beta)=\frac{\mathrm{RSS}}{\sigma^2}+n\log(2\pi\sigma^2).$$於是
$$\mathrm{AIC}=\frac{\mathrm{RSS}+2r\sigma^2}{\sigma^2}+\text{常數},$$和 $C_p$ 只差一個正倍數與常數,排序完全相同。這時 AIC 的「$\approx k$」變成等號:第一部分的 $2r\sigma^2$ 除以 $\sigma^2$ 正好是 $2r$。這不是巧合:高斯誤差下,負對數概似就是換了量尺的平方誤差。課本的 AIC 公式 $(\mathrm{RSS}+2d\hat\sigma^2)/(n\hat\sigma^2)$ 就是這個版本,把 $\sigma^2$ 換成共用的 $\hat\sigma^2$(省略共同的截距項)。若各模型各自重估 $\sigma^2$,用的是上面 $n\log(\mathrm{RSS}/n)$ 的版本,數值不同,排序通常接近,但不保證相同。
BIC 不估預測誤差,它問的是:哪個模型最可能產生這份資料?這裡 $M$ 代表一個候選模型(例如「用某一組 $d$ 個預測變數的線性迴歸」),$y$ 是觀測到的訓練反應向量。$P(M)$ 是看資料之前對 $M$ 的先驗機率,$P(M\mid y)$ 是看到 $y$ 之後的後驗機率;$\pi(\theta\mid M)$ 是模型 $M$ 內參數 $\theta$ 的先驗密度,$f(y\mid\theta,M)$ 是模型 $M$ 在參數 $\theta$ 下產生 $y$ 的機率密度(也就是概似)。由貝氏定理,
$$P(M\mid y)\propto P(y\mid M)\,P(M),\qquad P(y\mid M)=\int f(y\mid\theta,M)\,\pi(\theta\mid M)\,d\theta.$$$P(y\mid M)$ 稱為邊際概似(模型證據),是概似依先驗分布加權的平均。各模型的先驗機率相同時,比後驗機率就是比邊際概似。
Laplace 近似。$n$ 大時,概似集中在 $\hat\theta$ 附近:$\ell(\theta)\approx\ell(\hat\theta)-\tfrac12(\theta-\hat\theta)^\top J_n(\theta-\hat\theta)$,$J_n=-\nabla^2\ell(\hat\theta)$ 是觀測資訊矩陣(對數概似在最大值處的負二階導數),$n$ 大時 $J_n\approx nI(\hat\theta)$。先驗在這個小範圍內近乎常數,積分變成高斯積分:
$$P(y\mid M)\approx e^{\ell(\hat\theta)}\,\pi(\hat\theta\mid M)\,(2\pi)^{k/2}\,|J_n|^{-1/2},\qquad |J_n|\approx n^k\,|I(\hat\theta)|.$$$|\cdot|$ 是矩陣的行列式,$k$ 是模型 $M$ 的參數個數。
取對數:
$$\log P(y\mid M)\approx\ell(\hat\theta)-\frac k2\log n+\underbrace{\log\pi(\hat\theta\mid M)+\frac k2\log(2\pi)-\frac12\log|I(\hat\theta)|}_{\text{不隨 }n\text{ 增長}}.$$略去不隨 $n$ 增長的項,再乘上 $-2$:
$$\mathrm{BIC}=-2\ell(\hat\theta)+k\log n.$$直覺:後驗在每個參數方向上只集中在寬度約 $1/\sqrt n$ 的範圍,先驗卻把機率攤在固定大小的範圍。參數愈多,先驗裡真正和資料吻合的比例愈小(約 $n^{-k/2}$),邊際概似就被這個比例拉低。複雜模型的最大概似雖然高,對整個參數空間平均起來卻不一定高。
高斯線性迴歸同樣化成 $\mathrm{BIC}=n\log(\mathrm{RSS}/n)+d\log n+\text{常數}$;$\sigma^2$ 已知時是 $\mathrm{RSS}/\sigma^2+r\log n+\text{常數}$,代入共用的 $\hat\sigma^2$ 就得到課本的 $(\mathrm{RSS}+\log(n)\,d\hat\sigma^2)/n$(差一個正倍數與常數)。所以 $\log n$ 不是從預測誤差的樂觀偏差推出來的,而是對參數積分的結果。
| $C_p$ | AIC | BIC | |
|---|---|---|---|
| 想估計的量 | 固定設計下的測試 RSS | 新資料的期望對數概似(KL 散度的相對大小) | 邊際概似 $P(y\mid M)$ |
| 懲罰的來源 | 訓練 RSS 的樂觀偏差 $2r\sigma^2$ | 最大對數概似的樂觀偏差 $k$ | Laplace 近似的 $\tfrac k2\log n$ |
| 理論工具 | 投影矩陣的跡 | 最大概似估計的漸近常態 | Laplace 近似 |
| 每多一個參數的懲罰($-2\ell(\hat\theta)$ 量尺) | $2$ | $2$ | $\log n$ |
| 主要假設 | 固定設計、誤差同變異、模型不是用 $y$ 搜尋出來的 | 模型正確指定、大樣本 | 大樣本、先驗在 $\hat\theta$ 附近平滑且為正 |
目標與漸近性質(漸近效率、選模一致性)的對照見上方漸近效率與選模一致性。總結:$C_p$ 和 AIC 的 2 不是隨手挑的懲罰係數,而是訓練擬合相對於新資料的期望樂觀程度,可以算出來;BIC 的 $\log n$ 回答的是另一個問題(哪個模型最可能產生資料)。兩類準則的差別在目標,不只是罰多罰少。
調整後 $R^2$ 的精神:$R^2$ 問的是「加入預測變數後,訓練 RSS 比只用平均數預測時下降了多少比例」,只要加變數就不會下降。調整後 $R^2$ 改問「扣掉估計參數消耗的自由度之後,殘差變異數 $\mathrm{RSS}/(n-d-1)$ 是否真的比只用平均數時的 $\mathrm{TSS}/(n-1)$ 小」,加入的變數若沒有帶來足以抵銷自由度成本的 RSS 改善,它就會下降;不過無用的變數也可能因樣本波動剛好讓它上升。
設 $Y_i=\mu+\varepsilon_i$,$\operatorname{Var}(\varepsilon_i)=\sigma^2$。$\mu$ 已知時可用 $\frac1n\sum_i(Y_i-\mu)^2$ 估 $\sigma^2$;$\mu$ 未知、改用 $\bar Y$ 估時,殘差必須滿足 $\sum_i(Y_i-\bar Y)=0$,$n$ 個殘差只有 $n-1$ 個能自由變動。可以證明
$$E\Bigl[\sum_{i=1}^n(Y_i-\bar Y)^2\Bigr]=(n-1)\sigma^2,$$所以除以 $n-1$ 才得到 $\sigma^2$ 的不偏估計。
模型有 $d$ 個預測變數和一個截距,OLS 估 $d+1$ 個參數,殘差 $e$ 要滿足 $d+1$ 條線性限制 $X^\top e=0$,剩下 $n-d-1$ 個自由度。線性模型正確指定、設計矩陣滿秩、誤差變異數相同時,
$$E[\mathrm{RSS}]=(n-d-1)\sigma^2\quad\Longrightarrow\quad\hat\sigma^2=\frac{\mathrm{RSS}}{n-d-1}.$$除以自由度,是在修正「用同一批資料估參數,讓 RSS 系統性偏小」的現象。參數愈多,RSS 本來就會下降,這是預期中的事,不該每降一點都算成模型進步。
| $R^2$ | 調整後 $R^2$ | |
|---|---|---|
| 比較基準 | 只用平均數預測 | 只用平均數預測 |
| 衡量方式 | 平方誤差總和 | 自由度修正後的變異數 |
| 增加預測變數 | 不會下降 | 可能下降 |
| 主要目的 | 樣本內解釋程度 | 考慮模型複雜度後的改善 |
分母 $\mathrm{TSS}/(n-1)$ 是只用平均數預測時的樣本變異數;分子 $\mathrm{RSS}/(n-d-1)$ 是用了 $d$ 個預測變數之後的殘差變異數估計。因此
$$R^2_{\text{adj}}>0\iff\frac{\mathrm{RSS}}{n-d-1}<\frac{\mathrm{TSS}}{n-1},$$也就是模型的殘差均方比只用平均數的樣本變異數還小。這仍是訓練資料上的比較,不保證測試誤差比較小;自由度修正也沒有處理「從很多候選中搜尋變數」造成的過度擬合。新增一欄何時讓調整後 $R^2$ 上升,見下方的門檻證明。
形式上只差在每多一個變數的懲罰(上方公式與表格);BIC 的懲罰在 $n\ge8$ 時較重,所以比較同一組候選模型時,最優模型的變數數目不會較大(並列時採相同的簡約規則)。
背後的動機不同:$C_p$/AIC 想預測得準,BIC 想找出正確的模型。三個係數各自怎麼推導出來,見上方「懲罰項各從哪裡來」收合區;目標的比較見「漸近效率、選模一致性」收合區。
實務上的建議:目的是預測就用 CV(或 AIC/$C_p$),目的是「哪些變數真的有關」就參考 BIC。而且兩者都要求你估得出 $\hat\sigma^2$——這在 $p>n$ 的時候無法按全模型的常規公式計算(最後一節會講)。
$C_p$、AIC、BIC 都有明確的推導目標(見上方收合區)。調整後 $R^2$ 沒有這種東西。它就是「把 RSS 除以 $n-d-1$ 而不是 $n$」這個直覺上合理的修正。
它的直覺(上方主文)是對的,問題是 RSS 要降多少才算值得,沒有像 $C_p$ 那樣從預測誤差推出來的刻度,所以它的懲罰強度是任意的——在 Credit 上它比 $C_p$ 還鬆(選 7 個 vs 6 個)。
np.float64(3812.2130777982)
argmin() 回傳 5,而 bic_f 是 0-based 的 list,所以指的是清單的第 6 個元素——6 變數模型,BIC = 3812.21。注意下一格的 models1.model[5] 用的是 pandas 的標籤索引,取到的是 5 變數模型(輸出列出 5 個特徵)。同一個 5,兩種索引意義完全不同,讀別人的程式碼時要分清楚。
Ch06-varselect-lab-zh.ipynb
這是 adjust_r2() 的定義(就是式 6.4)加上在最佳子集前緣上取最大值。Hitters 上調整後 $R^2$ 選出 11 個變數;BIC 選 6 個。兩個準則差了 5 個變數,反映了它們對模型複雜度採用不同的懲罰強度。
Ch06-varselect-lab-zh.ipynb
比較準則前須使用同一候選集合並確認量尺。本節 BIC/$C_p$ 的懲罰比較採共用 $\hat\sigma^2$ 的設定;常數能否省略見下方「$C_p$ 的風險目標與常數」,軟體各自估變異數的概似版本見「先分清估計量與比較量尺」。
同一份資料上,BIC 選出的模型大小,跟 $C_p$ 選出的比起來?
先固定 $X$ 及候選線性模型,不用 $y$ 搜尋欄位。令 $H$ 為其投影矩陣,$r=\operatorname{rank}(H)$ 計入截距,並假設 $y=\mu+\varepsilon$、$\operatorname{Cov}(\varepsilon)=\sigma^2I$。以相同 $X$ 上一組獨立新反應 $y^{new}$ 定義固定設計預測風險,完整的 $C_p$ 型估計為
$$\widehat R=\frac{\mathrm{RSS}+2r\hat\sigma^2}{n}.$$這是對訓練反應重抽的期望風險,不是某次實現的測試 MSE。若所有模型都含截距、$d$ 只數斜率,$r=d+1$;比較大小時可省掉共同的 $2\hat\sigma^2/n$,但省略後不能再宣稱數值本身是完整風險的無偏估計。
若每個大小的最佳模型也是用同一份 y 搜尋出來的,$H$ 已依賴 $y$,單純代入 $r$ 不會自動修正搜尋的樂觀偏差。
$2r\sigma^2$ 這個樂觀偏差的推導(含近似偏差項為何相消),見上方三個準則的懲罰項推導第一部分。
獨立二元觀測的飽和模型逐筆令 $p_i=y_i$,以 $0\log0=0$ 慣例,其對數概似為0。故Bernoulli離差為 $D=-2\sum_i[y_i\log p_i+(1-y_i)\log(1-p_i)]$。對單筆二元標籤,可用 2 * log_loss(y, predict_proba(X), normalize=False) 計算;輸入必須是機率,不能把log機率再次當機率,也不能漏掉總和與平均的$n$倍差別。
這是固定模型的擬合量;有懲罰的係數估計不等於無懲罰MLE,不能直接把任意兩個正則化模型的離差差值當成一般卡方概似比檢定。
離差評估的是模型為已觀察標籤分配了多少機率。預測相同類別的兩個模型,若一個機率接近 0.5、另一個非常接近正確標籤,分類錯誤率可以相同而離差不同;錯誤且過度自信時,對數損失會很大。二元逐筆 Bernoulli 的飽和概似為零,故離差是總對數損失的兩倍;合併成分組二項計數時,飽和項通常不再是零,要重新用完整定義。
sklearn 的 predict_proba 欄序由 classes_ 決定;log_loss 的平均或總和由 normalize 決定。相同樣本數時兩者排序一致,數值仍相差 $n$ 倍。若要求 AIC/BIC 的常規最大概似版本,也要檢查邏輯斯模型是否取消預設正則化,而非直接把懲罰擬合當 MLE。
來源:log_loss 定義、Binomial deviance。
講義延伸來源:stackoverflow.com。
只要兩個模型使用同一批觀測、相同平方損失,而舊模型是新模型的子模型,新模型可以把新增係數設為零,得到舊模型的擬合。新模型再最佳化後不可能更差。理由是新增變數後,舊模型仍是一個可行解;這個論證不需要誤差常態,但只適用於訓練 RSS,不能推出測試誤差也會下降。移除不同的缺失列或改成正則化目標後,不能直接套用這個比較。
AIC 與 $C_p$ 何時排序相同、各模型重估 $\sigma^2$ 時的 $n\log(\mathrm{RSS}/n)$ 版本,推導都在上方三個準則的懲罰項推導。statsmodels 的 .aic、.bic 屬於重估 $\sigma^2$ 的對數概似尺度,並不是把課本的 $C_p$ 數值改名。只在同一資料、同一反應及相容概似設定下比較;模型用的列數改變,AIC 數字就失去原本的可比性。是否把共同截距與誤差變異數算進 $k$,影響數值常數;所有模型採相同慣例時,共同項不影響排名。
固定 $n$ 與 TSS,最大化調整後 $R^2$ 等於最小化 $\mathrm{RSS}/(n-d-1)$。舊模型有 $d$ 個斜率,再加一個欄後,必須讓 RSS 的下降抵得過少一個殘差自由度,調整後 $R^2$ 才會上升。所以 RSS 降了一點,還不夠讓調整後 $R^2$ 偏好加入新欄。
令 $R_0$ 是舊 RSS,$R_1$ 是新 RSS,舊殘差自由度為 $\nu=n-d-1$,新自由度為 $\nu-1>0$。調整後 $R^2$ 增加等價於
$$\frac{R_1}{\nu-1}\lt\frac{R_0}{\nu}\quad\Longleftrightarrow\quad R_0-R_1\gt\frac{R_0}{\nu}.$$用新模型殘差變異數寫出的單欄部分 F 值為 $(R_0-R_1)/(R_1/(\nu-1))$;上述條件也等價於 $F>1$。這只是選模門檻,不等於常用顯著水準下的假設檢定門檻。
CV 的選擇目標由 scoring 決定:用 MSE、$R^2$ 或分類對數損失可能得到不同的模型。
延伸對照:OLSResults.aic、OLSResults.bic、調整後 $R^2$ 定義。
講義延伸來源:stats.stackexchange.com。
講義延伸來源:stats.stackexchange.com。
講義延伸來源:stats.stackexchange.com。
講義延伸來源:en.wikipedia.org。
$C_p$、AIC、BIC 都需要估 $\hat\sigma^2$,而在 $p$ 接近 $n$ 的時候那個估計本身就不可靠。 交叉驗證不需要 $\hat\sigma^2$、不需要知道自由度、也不需要假設模型是對的。 它直接估測試誤差。所以只要算得動,CV 是首選。
做法就是第 5 章那一套,只是把「模型」換成「模型大小」:對每個大小 $k$, 在每一折的訓練部分做子集選擇、在驗證折上算誤差。 變數選擇必須只使用各折的訓練資料。這正是第 5 章 P06 講的那個錯誤。
one-SE 規則裡的標準誤,通常用 $K$ 個折的誤差估計:$\widehat{\mathrm{SE}}=\mathrm{sd}(E_1,\ldots,E_K)/\sqrt K$。各折的訓練資料重疊、誤差並不獨立,所以它只是粗略的不確定性指標。
設第 $j$ 折的驗證誤差為 $E_j$(例如該折的 MSE),CV 估計是 $\bar E=\frac1K\sum_jE_j$(各折大小不同時用講義第 5 章的加權式)。把 $E_1,\ldots,E_K$ 當成 $K$ 個觀測,平均數的標準誤為
$$\widehat{\mathrm{SE}}(\bar E)=\frac{1}{\sqrt K}\sqrt{\frac{1}{K-1}\sum_{j=1}^K(E_j-\bar E)^2}.$$Lab 用 NumPy 的 std(),預設 ddof=0(除以 $K$);改成 ddof=1 除以 $K-1$,誤差棒會略長。這個公式把各折當成獨立,但任兩折的訓練集共用約 $(K-2)/K$ 的資料,$E_j$ 彼此相關,所以真實的變異通常被低估。它適合拿來看「曲線在哪裡平坦、差異是否小於誤差棒」,不能當成嚴格的信賴區間。
單次 CV 的結果還受「這一次怎麼切」影響。RepeatedKFold(n_splits=K, n_repeats=R) 用 $R$ 種不同的隨機切分各做一次 $K$ 折,得到 $R$ 個 CV 估計 $\bar E^{(1)},\ldots,\bar E^{(R)}$。它們的平均比單次 CV 穩定;它們的標準差描述的是切分造成的波動,不是對新資料的抽樣變異,因為 $R$ 次用的是同一份資料。若把 $R\times K$ 個折誤差全部當獨立觀測、再除以 $\sqrt{RK}$,會更嚴重地低估不確定性。
程式假設已備妥 X、y;比較不同模型大小時,對每個大小重複同一段計算,並讓每個模型用同一組切分,差異才不會混入切分的運氣。分類問題可改用 RepeatedStratifiedKFold。
Pipeline 是關鍵:把標準化包進去,GridSearchCV 就會在每一折的訓練部分重新 fit scaler,而不是用全部資料的平均與標準差。這一步做錯,CV 誤差會偏低而且不會報錯。
Ch06-varselect-lab-zh.ipynb
三個理由。第一,運算成本:最佳子集有 $2^p$ 個模型,每個都跑 k-fold 是 $k \cdot 2^p$ 次擬合;而 $C_p$ 只要一次擬合加一個修正項。第二,這三個準則會出現在別人的論文與報表裡,你得看得懂。第三,它們示範了「懲罰複雜度」這個想法,目標的差別見上一節。
但如果你只是要挑一個模型來預測,而且資料量算得動:用 CV。
會,而且那是刻意的。它的立場是:在最佳 CV 誤差加一個標準誤的範圍內,偏好簡單。
這個偏好有實際理由——變數少的模型更好解釋、需要蒐集的資料更少、上線後更不容易因為某個變數的定義改變而壞掉。如果你的目標純粹是最小化預測誤差而完全不在意這些,那就選最低點;但要記得那個「最低」在下一份資料上很可能換人。
設第 $b$ 折、大小 $k$ 的保留資料 MSE 為 $E_{bk}$,先算 $K$ 折的平均,再找最小者 $k^*$。one-SE 規則把可接受區間的上緣設為最佳模型的平均誤差加上最佳點的標準誤,最後取仍在區間內的最小 $k$:
$$\hat k=\min\{k:\bar E_k\le \bar E_{k^*}+\widehat{\mathrm{SE}}(\bar E_{k^*})\}.$$標準誤怎麼估、它為什麼只是粗略指標,以及 RepeatedKFold 的做法,見上方的CV 標準誤收合區。
講義 Credit 的驗證切分用四分之三訓練、四分之一驗證,CV 用十折;Hitters Lab 則示範五折及 20% 的單次驗證。二者都比較模型大小,但數據與分折設定不同,不能拿誤差值相互核對。
講義延伸來源:www.cs.cmu.edu。
子集選擇是離散的:一個變數只有納入或排除兩種選擇。這讓它變異很大—— 資料動一點,選出來的變數就換一批。另一條路是留下全部變數,但把係數往零壓。
Ridge 在原本的 RSS 上加一個 L2 懲罰:
$$\sum_{i=1}^{n}\left(y_i - \beta_0 - \sum_{j=1}^{p}\beta_j x_{ij}\right)^2 + \lambda \sum_{j=1}^{p} \beta_j^2 \;=\; \mathrm{RSS} + \lambda \|\beta\|_2^2$$$\lambda = 0$ 就是最小平方;$\lambda \to \infty$ 所有係數趨近 0。 注意截距 $\beta_0$ 不罰——罰它等於在乎你把 $y$ 的原點放在哪裡。
np.float64(160.42371017725918)
ElasticNet 的 l1_ratio=0 就是純 Ridge。係數的 L2 範數會隨 $\lambda$ 單調下降。這個數字就是上圖 x 軸比值的分子。
Ch06-varselect-lab-zh.ipynb
RidgeCV 把「掃 $\lambda$ + 交叉驗證」包成一步。alphas 接收一組候選值;排列順序本身不表示 RidgeCV 使用暖啟動。這張 lab 卡把 scaler 放在內建 CV estimator 外層,因此 scaler 會先對傳入 pipeline 的全部資料擬合,RidgeCV 才在縮放後的資料內切折。
Ch06-varselect-lab-zh.ipynb
StandardScaler 放在 RidgeCV 或 ElasticNetCV 外層,
只能保證 scaler 不會看到 pipeline 之外的外部測試集;內建 estimator 切 CV 折以前,資料已經縮放完畢。
要讓每一折都只用該折訓練部分估平均與標準差,讓 GridSearchCV 搜尋一條包含
StandardScaler 與 Ridge 的完整 pipeline。
GridSearchCV 會在每一折複製並擬合整條 pipeline,因此該折的 scaler
只看該折訓練資料。這段概念碼假設已備妥外層訓練資料 X_train、Y_train;
三個 $\alpha$ 是示意候選值,這裡以 5 折平均 MSE 選擇。外部測試集留到選模完成後再使用。
把某個預測變數的單位從「元」改成「千元」(數值除以 1000)。Ridge 的預測會變嗎?
先用訓練資料把每個 $x_{ij}$ 與 $y_i$ 減去各自的平均,截距就可以分開處理。斜率最小化 $\sum_{i=1}^n\big(y_i-\sum_{j=1}^p\beta_jx_{ij}\big)^2+\lambda\sum_{j=1}^p\beta_j^2$,對每個 $\beta_j$ 偏微分設為 0,得到 $p$ 條線性方程(Ridge 的正規方程):
$$\sum_{k=1}^p\Big(\sum_{i=1}^nx_{ij}x_{ik}\Big)\beta_k+\lambda\beta_j=\sum_{i=1}^nx_{ij}y_i,\qquad j=1,\ldots,p.$$$\lambda\gt0$ 時,不論原始變數是否共線、甚至 $p\gt n$,這組方程都有唯一解。截距用 $\hat\beta_0=\bar y-\sum_j\hat\beta_j\bar x_j$ 還原;截距不能跟斜率一起用同一個 $\lambda$ 懲罰。Ridge 在各主成分方向上怎麼收縮,見 PCR 節的Ridge 與 PCR 的對照。
目標對 $\beta_j$ 的偏導數是 $-2\sum_ix_{ij}\big(y_i-\sum_kx_{ik}\beta_k\big)+2\lambda\beta_j$,設為零整理即得上面的方程。
唯一性:任取不全為零的 $v_1,\ldots,v_p$,沿著方向 $(v_1,\ldots,v_p)$ 的二階導數為
$$2\Big[\sum_{i=1}^n\Big(\sum_{j=1}^px_{ij}v_j\Big)^2+\lambda\sum_{j=1}^pv_j^2\Big]\gt0 .$$第一項可能為 0(例如變數共線時),但 $\lambda\sum_jv_j^2$ 一定為正,所以目標是嚴格凸函數,最小值只有一個。$\lambda=0$ 時少了這一項,共線就會讓解不唯一;最小平方在 $p\gt n$ 時的問題就出在這裡。
Ridge 的缺點是:它留下全部 $p$ 個變數。$p=500$ 的時候你得到 500 個很小但不為零的係數, 模型一點也不好解釋。Lasso 把 L2 換成 L1:
$$\mathrm{RSS} + \lambda \sum_{j=1}^{p} |\beta_j| \;=\; \mathrm{RSS} + \lambda \|\beta\|_1$$只差一個平方,行為完全不同:$\lambda$ 夠大時 Lasso 會把一部分係數壓成剛好 0, 於是它同時做了收縮與變數選擇。這種解叫做稀疏(sparse)。
講義從三個角度回答「為什麼 Lasso 會歸零、Ridge 不會」:限制式的幾何、正交設計下的閉式解,以及貝氏解讀。
同一個收縮問題有兩種寫法。限制式(primal form)先給一個預算 $s$,只在預算內找 RSS 最小的係數;懲罰式(Lagrange form)把預算改成加在 RSS 後面的代價,由 $\lambda$ 控制。講義的 Lasso 示意圖(RSS 等高線碰上預算 $s$ 的限制區域)畫的就是限制式。
$$\underbrace{\min_{\beta}\ \sum_{i=1}^n\Big(y_i-\beta_0-\sum_{j=1}^p\beta_jx_{ij}\Big)^2\ \ \text{s.t.}\ \sum_{j=1}^p|\beta_j|\le s}_{\text{限制式}} \qquad\Longleftrightarrow\qquad \underbrace{\min_{\beta}\ \mathrm{RSS}+\lambda\sum_{j=1}^p|\beta_j|}_{\text{懲罰式}}$$Ridge 把 $|\beta_j|$ 換成 $\beta_j^2$,限制區域從菱形變成圓。最佳子集也能寫成限制式,只是預算數的是「有幾個係數不為 0」:
$$\min_{\beta}\ \mathrm{RSS}\ \ \text{s.t.}\ \sum_{j=1}^pI(\beta_j\ne0)\le s .$$這種限制稱為 $L_0$ 限制(常被叫做 $L_0$ 範數,但它不滿足範數的齊次性,嚴格說不是範數)。所以三個方法可以排成 $L_0$(最佳子集)、$L_1$(Lasso)、$L_2$(Ridge),差別只在預算怎麼量。$s$ 與 $\lambda$ 可以互相對應,但一般不是一對一。
以 Lasso 為例,令 $P(\beta)=\sum_j|\beta_j|$(Ridge 換成 $\sum_j\beta_j^2$)。
從懲罰式到限制式。固定 $\lambda\gt0$,設 $\hat\beta^\lambda$ 是懲罰式的解,令 $s=P(\hat\beta^\lambda)$。任取滿足 $P(\beta)\le s$ 的 $\beta$,因為 $\hat\beta^\lambda$ 讓懲罰式最小,
$$\mathrm{RSS}(\beta)+\lambda P(\beta)\ \ge\ \mathrm{RSS}(\hat\beta^\lambda)+\lambda s\ \ge\ \mathrm{RSS}(\hat\beta^\lambda)+\lambda P(\beta),$$兩邊消去 $\lambda P(\beta)$ 得 $\mathrm{RSS}(\beta)\ge\mathrm{RSS}(\hat\beta^\lambda)$,所以 $\hat\beta^\lambda$ 也是預算 $s$ 下的限制式解。
從限制式到懲罰式。RSS 與 $P$ 都是凸函數,限制式的解滿足 KKT 條件:存在 Lagrange 乘數 $\lambda\ge0$,使同一個解也讓 $\mathrm{RSS}+\lambda P(\beta)$ 最小。若預算已經放得下最小平方解($s\ge P(\hat\beta^{\mathrm{OLS}})$),限制沒有作用,對應的是 $\lambda=0$。
為什麼不是一對一。上面的對應依資料而定:所有夠大的 $s$ 都對應到 $\lambda=0$;解若不唯一,同一個 $\lambda$ 也可能對到不只一個 $s$。所以只能說兩種寫法描述同一族解,不能把某個 $\lambda$ 換算成固定的 $s$。
$L_0$ 的懲罰式。最佳子集也有懲罰式 $\mathrm{RSS}+\lambda\sum_jI(\beta_j\ne0)$。前面選模準則的 $C_p$ 就是這個形式:$n\,C_p=\mathrm{RSS}+2\hat\sigma^2d$,等於取 $\lambda=2\hat\sigma^2$;BIC 對應 $\lambda=\hat\sigma^2\log n$。因為 $L_0$ 不是凸的,KKT 的反方向不成立:懲罰式只會選到「RSS 對模型大小 $d$ 曲線」下凸包上的大小,限制式能選的大小則更多。
講義延伸來源:stats.stackexchange.com。
三種限制式的計算難度差很多:$L_0$ 是非凸的組合搜尋,最難;Lasso 把它換成凸的 $L_1$,可以有效率地求解;Ridge 最容易,有閉式解,選 $\lambda$ 時還能用留一捷徑或廣義交叉驗證(GCV),不必真的重擬合。
窮舉法要列出所有不超過 $s$ 個變數的子集,對每個子集做一次最小平方,再比較 RSS。$s=p$ 時就是子集選擇一節的 $2^p$ 個模型:
| $p$ | $s$ | 要檢查的子集數 $\sum_{k=0}^{s}\binom pk$ |
|---|---|---|
| 20 | 10 | 616,666(約 62 萬) |
| 30 | 15 | 614,429,672(約 6.1 億) |
| 50 | 10 | 13,432,735,556(約 134 億) |
| 50 | 50 | $2^{50}\approx1.13\times10^{15}$ |
一般情況下,最佳子集選擇屬於 NP-hard 問題,但這不表示所有情況都算不出來:$p$ 小時可以窮舉,分支界定法與混合整數最佳化也能處理部分中大型問題,只是最壞情況下仍無法有效率地保證全域最佳。$\sum_jI(\beta_j\ne0)$ 是不連續、非凸的函數;兩個各只有一個非零係數的解取平均,非零個數會變成兩個,限制區域因此不是凸集合。
$\sum_j|\beta_j|$ 是凸函數,Lasso 的目標因此是凸的,任何局部最小值都是全域最小值。它沒有閉式解,求解方法見本節後面的座標更新與 LAR。
Ridge 的懲罰式可微,對 $\beta_j$ 偏微分設為 0,得到 $p$ 條線性方程,$\lambda\gt0$ 時解唯一(見 Ridge 節的正規方程)。解出的擬合值是 $y$ 的線性組合,$\hat y_i=\sum_{k=1}^nh_{ik}(\lambda)\,y_k$,其中權重 $h_{ik}(\lambda)$ 只由 $x$ 與 $\lambda$ 決定。第 5 章的留一捷徑對 Ridge 同樣成立:
$$\mathrm{CV}_{(n)}(\lambda)=\frac1n\sum_{i=1}^n\left(\frac{y_i-\hat y_i}{1-h_{ii}(\lambda)}\right)^2 .$$截距不懲罰時,$h_{ii}(\lambda)$ 包含截距貢獻的 $1/n$,加總為 $\sum_ih_{ii}(\lambda)=1+\mathrm{df}(\lambda)$,其中 $\mathrm{df}(\lambda)$ 是斜率部分的有效自由度(不含截距,與 ESL 的慣例相同)。GCV 再把每個 $h_{ii}(\lambda)$ 換成它們的平均:
$$\mathrm{GCV}(\lambda)=\frac1n\sum_{i=1}^n\left(\frac{y_i-\hat y_i}{1-\{1+\mathrm{df}(\lambda)\}/n}\right)^2 .$$兩者都只要在每個 $\lambda$ 擬合一次,不必真的做 $n$ 次或 $K$ 次重擬合。$\mathrm{df}(\lambda)$ 用主成分寫出來的形式見 PCR 節的Ridge 與 PCR 的對照。sklearn 的 RidgeCV 在預設 cv=None 時就是用這種高效率的留一計算。Lasso 的擬合值不是 $y$ 的固定線性組合(哪些係數為 0 也由 $y$ 決定),所以沒有同樣精確的捷徑。
回到幾何。RSS 的等高線是以 $\hat\beta^{\mathrm{OLS}}$ 為中心的同心橢圓。解就是橢圓第一次碰到限制區域的那一點。菱形有尖角而且尖角落在座標軸上,碰撞很容易發生在尖角,那裡有一個座標是 0。圓沒有尖角,碰到座標軸是機率 0 的事。
尖角為什麼這麼容易被碰到?看菱形 $|\beta_1|+|\beta_2|\le1$ 的尖角 $P=(1,0)$。交於 $P$ 的兩條邊,外法線分別是 $(1,1)$ 與 $(1,-1)$;夾在這兩個方向之間的所有方向,都可以當作 $P$ 的支撐法線,合起來稱為 $P$ 的法線扇形(normal cone)。等高線是圓時,只要 $\hat\beta^{\mathrm{OLS}}$ 落在從 $P$ 張開的這個扇形裡,Lasso 解就是 $P$,$\beta_2$ 剛好為 0。扇形是一整塊區域,所以很多不同的最小平方解,都會對應到同一個稀疏解。圓在每一點只有一個法線方向,沒有扇形可以收留這麼多 OLS 解。
取 $\mathrm{RSS}=(\beta_1-2)^2+(\beta_2-0.4)^2$,最小平方解是 $(2,0.4)$,限制為 $|\beta_1|+|\beta_2|\le1$。Lasso 解是 $P=(1,0)$。
Lasso:沿菱形的邊離開 $P$。想讓 $\beta_2$ 從 0 增加一點點 $\varepsilon\gt0$,菱形限制迫使 $\beta_1$ 至少降到 $1-\varepsilon$。RSS 的變化為
$$\begin{aligned}\Delta\mathrm{RSS}&=(1-\varepsilon-2)^2+(\varepsilon-0.4)^2-\big[(1-2)^2+0.4^2\big]\\&=1.2\,\varepsilon+2\varepsilon^2\gt0 .\end{aligned}$$$\beta_2$ 確實更接近最小平方的 0.4,但要犧牲 $\beta_1$,讓它離 2 更遠;兩者相抵後 RSS 反而上升,所以解留在尖角。
扇形條件。把 $\hat\beta^{\mathrm{OLS}}-P$ 寫成兩條外法線的組合:$\hat\beta^{\mathrm{OLS}}-P=a(1,1)+c(1,-1)$,扇形內就是 $a\ge0$、$c\ge0$。解得 $a=\tfrac12(\hat\beta_1-1+\hat\beta_2)$、$c=\tfrac12(\hat\beta_1-1-\hat\beta_2)$,所以條件是
$$|\hat\beta_2^{\mathrm{OLS}}|\le\hat\beta_1^{\mathrm{OLS}}-1 .$$本例 $\hat\beta_1^{\mathrm{OLS}}=2$,只要 $|\hat\beta_2^{\mathrm{OLS}}|\le1$,例如改成 $-0.5$ 或 $0.7$,解都仍是 $(1,0)$。嚴格的說法是:在 $P$ 點,RSS 下降最快的方向 $-\nabla\mathrm{RSS}=2(\hat\beta^{\mathrm{OLS}}-P)$ 落在法線扇形內;不是 RSS 梯度本身落在扇形內。
Ridge:沿圓周離開 $(1,0)$。限制換成 $\beta_1^2+\beta_2^2\le1$。把 $\beta_2$ 增加 $\varepsilon$,只需讓 $\beta_1=\sqrt{1-\varepsilon^2}\approx1-\varepsilon^2/2$。同一個 RSS 的變化變成
$$\Delta\mathrm{RSS}=4-4\sqrt{1-\varepsilon^2}-0.8\,\varepsilon\approx-0.8\,\varepsilon+2\varepsilon^2\lt0\quad(\varepsilon\text{ 很小}),$$RSS 會下降,所以 Ridge 不會停在 $(1,0)$。
| Lasso(菱形) | Ridge(圓) | |
|---|---|---|
| $\beta_2$ 增加 $\varepsilon$ | $\beta_1$ 至少減少 $\varepsilon$ | $\beta_1$ 只需減少約 $\varepsilon^2/2$ |
| 離開座標軸的代價 | 一階,$O(\varepsilon)$ | 二階,$O(\varepsilon^2)$ |
尖角容易產生稀疏解,原因就在這裡:Lasso 離開座標軸必須付一階代價,只要 RSS 的一階好處不夠大,就寧願讓 $\beta_2=0$;Ridge 在座標軸附近只付二階代價,只要最小平方解的 $\beta_2\ne0$,稍微離開座標軸幾乎總是划算。
考慮最簡單的情況:$n=p$,$X$ 是對角元素為 1 的對角矩陣,並省略截距。最小平方問題變成逐個座標最小化 $\sum_{j=1}^p(y_j-\beta_j)^2$,解為 $\hat\beta_j=y_j$。加上兩種懲罰後:
$$\text{Ridge:}\ \hat\beta_j^R=\frac{y_j}{1+\lambda},\qquad \text{Lasso:}\ \hat\beta_j^L=\begin{cases}y_j-\lambda/2, & y_j>\lambda/2,\\ y_j+\lambda/2, & y_j<-\lambda/2,\\ 0, & |y_j|\le\lambda/2.\end{cases}$$Ridge 把每個係數按同一比例往 0 縮,$y_j$ 再小也只是變小,不會等於 0。Lasso 則是軟門檻(soft-thresholding):每個係數往 0 平移 $\lambda/2$,絕對值不到 $\lambda/2$ 的直接變成 0。例如 $\lambda=2$、$y_j=0.8$:Ridge 給 $0.8/3\approx0.27$,Lasso 給 $0$;$y_j=3$ 時 Ridge 給 $1$,Lasso 給 $2$。講義指出,一般的 $X$ 沒有這麼乾淨的公式,但「Ridge 等比例收縮、Lasso 軟門檻」的主要想法仍近似成立。
Ridge。逐座標最小化 $(y_j-\beta_j)^2+\lambda\beta_j^2$。微分得 $-2(y_j-\beta_j)+2\lambda\beta_j=0$,所以 $\beta_j=y_j/(1+\lambda)$;二階導數 $2(1+\lambda)>0$,是唯一最小值。
Lasso。逐座標最小化 $g(\beta_j)=(y_j-\beta_j)^2+\lambda|\beta_j|$。分三種情況:
三種情況正好把實數線分完,合併就是 $\hat\beta_j^L=\operatorname{sign}(y_j)\,(|y_j|-\lambda/2)_+$。關鍵在第三種情況:$|\beta_j|$ 在 0 有尖角,次梯度是一整段區間 $[-1,1]$,只要 $|y_j|$ 落在 $\lambda/2$ 以內,0 就滿足最優條件。$\beta_j^2$ 在 0 是平滑的,導數只有 0 一個值,只有 $y_j=0$ 時解才會是 0。一般 $X$ 怎麼求解,見本節後面的座標更新。
貝氏觀點下,後驗正比於概似乘以先驗:$p(\beta\mid X,Y)\propto f(Y\mid X,\beta)\,p(\beta)$,並假設 $p(\beta)=\prod_{j=1}^pg(\beta_j)$。若 $g$ 是平均為 0 的常態分布,後驗眾數(給定資料下最可能的 $\beta$)就是 Ridge 的解;若 $g$ 是平均為 0 的雙指數(Laplace)分布,後驗眾數就是 Lasso 的解。兩種先驗的尺度參數都對應到 $\lambda$。
選 Laplace 先驗,等於事前就認為多數係數應該很接近 0。變異數相同時,Laplace 在 0 附近放的機率密度比常態高,尾部又比常態厚:大部分係數被拉到接近 0,少數係數仍可以保持很大。常態先驗的尾部衰減得快,比較傾向讓每個係數都中等大小。
Laplace 密度在 0 有尖峰,取對數後就是 $-|\beta_j|/b$ 這個在 0 有尖角的懲罰,所以後驗眾數可能剛好停在 0,得到稀疏解;常態先驗在 0 是平滑的圓頂,後驗眾數通常只會被拉近 0。這和解釋二的「尖角」是同一件事。下圖可以切換成負對數密度,直接看到兩種先驗各自對應的懲罰形狀。
獨立常態誤差、固定誤差變異數與不懲罰截距的設定下,高斯斜率先驗的 MAP 對應 Ridge,Laplace 斜率先驗的 MAP 對應 Lasso。其 RSS 目標的 $\lambda$ 分別為 $\sigma^2/\tau^2$ 與 $2\sigma^2/b$。完整常數對應見下方證明。
獨立 $N(\beta_0+\sum_j\beta_jx_{ij},\sigma^2)$ 反應的負對數概似,除去與 $\beta$ 無關的常數後為 $\mathrm{RSS}/(2\sigma^2)$。各斜率獨立 $N(0,\tau^2)$ 的負對數先驗為 $\sum_j\beta_j^2/(2\tau^2)$;負對數後驗相加後乘以 $2\sigma^2$,成為 $\mathrm{RSS}+(\sigma^2/\tau^2)\sum_j\beta_j^2$。
Laplace 密度為 $e^{-|\beta_j|/b}/(2b)$,負對數先驗為 $\sum_j|\beta_j|/b$ 加常數,因此乘以 $2\sigma^2$ 後 $\lambda=2\sigma^2/b$。截距使用不懲罰的設定;$\sigma$、$\tau$、$b$ 在此視為固定。這是後驗眾數,後驗平均須積分計算,一般不會因 Laplace 先驗精確歸零。
Laplace 是連續分布,先驗本身沒有在 $\beta=0$ 放正機率的點質量,所以不能解讀成「每個係數事前有某機率恰為零」。原作者教材:Bayesian Regression。
講義延伸來源:ekamperi.github.io。
array([-210.01008773, 243.4550306 , 0. , 0. ,
0. , 97.69397357, -41.52283116, -0. ,
0. , 39.62298193, 205.75273856, 124.55456561,
-126.29986768, 15.70262427, -59.50157967, 75.24590036,
21.62698014, -12.04423675, -0. ])這張卡用的是 Lab 的 Hitters 資料,和上圖的 Credit 不同。注意輸出裡有剛好等於 0 的係數。這是 Lasso 跟 Ridge 最實際的差別:你可以直接說「這些變數被模型丟掉了」。
來源:Ch06-varselect-lab-zh.ipynb
講義用兩張圖,把 Lasso 與 Ridge 的平方偏差、變異與測試 MSE 畫在一起。第一張的模擬資料中 45 個預測變數都和反應有關,Ridge 的測試 MSE 略低;第二張只有兩個變數有關,Lasso 明顯較好(何時選哪一個,見下方 Q&A「Lasso 一定比 Ridge 好嗎?」)。
兩張圖的橫軸都不是 $\lambda$,而是訓練資料上的 $R^2$。兩個方法的 $\lambda$ 量尺不同,同一個數字代表的收縮程度不一樣,直接用 $\lambda$ 對齊沒有意義。訓練 $R^2$ 是兩者共同的擬合指標:$\lambda$ 變大、正則化變強,訓練 $R^2$ 就下降。用它當橫軸,可以比較兩個方法在相同訓練擬合程度下的測試 MSE;它也可以當作模型複雜度的間接指標,沿著正則化路徑看偏差與變異怎麼交換。
對同一個方法,$\lambda$ 從大到小時,係數的限制愈來愈寬,訓練 RSS 不會上升,訓練 $R^2$ 因此單調增加:
| 訓練 $R^2$ | 正則化 | 有效的模型複雜度 |
|---|---|---|
| 低 | 強($\lambda$ 大) | 較低 |
| 高 | 弱($\lambda$ 小) | 較高 |
限制是:訓練 $R^2$ 不等於正式的有效自由度。兩個方法在相同訓練 $R^2$ 下,用掉的自由度不一定相同,所以跨方法時只能說「訓練擬合程度相同」,不能說「模型複雜度完全相同」。Ridge 的有效自由度 $\mathrm{df}(\lambda)$ 見三種限制式的計算成本。
不一定,取決於真實的係數結構。
如果真的只有少數幾個變數有用(真實係數稀疏),Lasso 佔優勢。它能把沒用的丟掉,省下估計它們所花的變異。如果很多變數都有小小的貢獻(係數密集),Ridge 通常較好——Lasso 會硬把一些真的有用的變數歸零,付出偏差的代價。
而你事前並不知道是哪一種。實務作法是兩個都跑、用 CV 比較,或直接用 elastic net(同時放 L1 與 L2 懲罰)讓資料自己決定混合比例。
不能。Lasso 的變數選擇沒有附帶 p 值,而且「先用資料選變數、再對選出來的變數做 t 檢定」會嚴重高估顯著性。那是同一份資料用了兩次,與第 5 章的 CV 誤用有相同的問題。
要做選後推論(post-selection inference)需要專門的方法(selective inference、debiased lasso 等),已經超出本課範圍。可以依實際的選擇結果說:「在這個 $\lambda$ 之下,Lasso 保留了這些變數」,而不是「這些變數顯著」。
Lasso 可能只保留其中一個,選擇也可能隨資料的小變動而改變;但這不是因為多留一個變數就必然多付一份 L1 懲罰。
若兩欄完全相同,係數 (2,0) 和 (1,1) 的預測相同,L1 懲罰也都等於 2,所以解的分配可能不唯一。高度相關但不完全相同時,微小差異及求解方式會影響最後保留哪一欄。解讀結果時,不能直接把被歸零的那欄判為沒有作用。
Ridge 的平方懲罰在這個例子分別為 4 和 2,因此偏好分散係數;只有完全相同且同尺度的兩欄,對稱的 Ridge 解才精確平分。Elastic net 可促進相關欄一起保留;group lasso 則對事先指定的群組做整組選擇。
為什麼 L1 約束區域(菱形)會讓解落在座標軸上,而 L2(圓)不會?
一般的 $X$ 沒有閉式解。套件常用座標下降:固定其他係數,只對 $\beta_j$ 求解,結果是解釋二的軟門檻:令部分殘差 $r_i^{(-j)}=y_i-\sum_{k\ne j}\beta_kx_{ik}$,先對 $\sum_ix_{ij}r_i^{(-j)}$ 取門檻 $\lambda/2$ 的軟門檻,再除以 $\sum_ix_{ij}^2$,即 $\hat\beta_j=\operatorname{sign}(u_j)\,(|u_j|-\lambda/2)_+\big/\sum_ix_{ij}^2$,其中 $u_j=\sum_ix_{ij}r_i^{(-j)}$;逐一更新直到收斂。不同套件若用 $\mathrm{RSS}/(2n)$,相同 $\lambda$ 的數字不能直接比較。
講義延伸來源:dafriedman97.github.io。
講義延伸來源:stats.stackexchange.com。
講義附錄只列出這兩個方法的連結,這裡各用一句話說明。LAR 從與殘差相關最大的特徵開始,沿活躍特徵的「等角」方向前進,其他特徵的相關追平時再加入;加上「係數走到零就移出」的規則,就能算出整條 Lasso 路徑(sklearn 的 lars_path、LassoLars)。
Group Lasso 以 $\lambda\sum_g\sqrt{p_g}\|\beta_g\|_2$ 懲罰事先指定的各組係數,讓同一類別變數的虛擬欄或同一組基底整組進出模型;選入的組內係數不必稀疏。
講義延伸來源:stats.stackexchange.com。
講義延伸來源:stats.stackexchange.com。
講義延伸來源:stats.stackexchange.com。
$\lambda$ 是調整參數(tuning parameter),不是從資料估出來的參數。 選它的方法跟第 5 章選多項式次數完全一樣:掃一排候選值,用交叉驗證比較。
Credit 資料上 Lasso 的 CV 選出 $\lambda=$ 364.76, 在那個 $\lambda$ 之下 11 個變數全部仍保留。這份資料的係數並不稀疏, 所以在這個設定下保留了全部變數。 Lasso 是否產生零係數,取決於資料與 $\lambda$ 的選擇。
為什麼用對數等距?正則化的效果取決於 $\lambda$ 的相對大小,而不是絕對增量:同樣加 1,$\lambda$ 從 0 到 1 會讓係數縮掉一半,從 99 到 100 幾乎沒有差別。等距網格會在模型幾乎不變的大 $\lambda$ 區域放太多點;對數網格讓每個數量級分到差不多的點數,各種正則化尺度都涵蓋得到。實作上用 np.logspace(-2, 5, 100) 或 10**np.linspace(a, b, 100)。
沿用 Lasso 節「解釋二」的對角設計(推導),Ridge 解為 $\hat\beta_j(\lambda)=y_j/(1+\lambda)$,所以係數保留比例是
$$r(\lambda)=\frac{\hat\beta_j(\lambda)}{y_j}=\frac1{1+\lambda}.$$| $\lambda$ | 係數保留比例 $r(\lambda)$ |
|---|---|
| 0 | 100% |
| 0.1 | 90.9% |
| 1 | 50.0% |
| 10 | 9.1% |
| 100 | 1.0% |
從 $\lambda=0$ 到 $1$,保留比例掉了一半;從 $\lambda=99$ 到 $100$,只從 1.00% 變成 0.99%。若搜尋範圍是 $[0.001,100]$,等距取 10 點時相鄰兩點差約 11.1,除了起點 0.001,其餘 9 點都落在 $\lambda\gt11$、保留比例不到 9% 的區域;對數等距時相鄰兩點的比值固定約 3.6,有 6 個點落在 $\lambda\lt1$。畫在 $\log\lambda$ 上,$r(\lambda)$ 是一條 S 形曲線,變化分散在好幾個數量級;畫在原尺度上,變化幾乎全擠在 $\lambda\lt10$ 的一小段。實際資料中,模型開始明顯改變的 $\lambda$ 事先不知道落在哪個數量級,所以在對數尺度上等距取點,讓每個數量級都分到點,也避免在模型幾乎不變的區域過度取樣。
np.float64(114690.73118253677)
mse_path_ 的形狀是($\lambda$ 個數 $\times$ 折數),.mean(1) 先對折取平均、再取最小值,就是 CV 曲線的最低點。這張 lab 卡同樣把 scaler 放在 ElasticNetCV 外層;若要在每一折內重新估縮放統計量,改用上面的 GridSearchCV(Pipeline(...)) 結構。
Ch06-varselect-lab-zh.ipynb
用 CV 選好 $\lambda$ 之後,最終交出的模型應該是?
sklearn Ridge 最小化 $\mathrm{RSS}+\alpha\|\beta\|_2^2$,因此與本章 Ridge 的 $\lambda$ 相同。Lasso 的目標是 $\mathrm{RSS}/(2n)+\alpha\|\beta\|_1$,乘以 $2n$ 後才對到本章 $\lambda=2n\alpha$。ElasticNet 用混合比 $\rho$,目標是
所以 $\rho=0$ 時對到 Ridge 的 $\alpha_{\text{Ridge}}=n\alpha_{\text{ElasticNet}}$;$\rho=1$ 對到 Lasso。不要把 R glmnet 的 alpha(混合比)與 sklearn 的 alpha(懲罰強度)混在一起。Lab 在 Ridge 路徑示範特別乘上樣本數以對照 ElasticNet,後面的手選網格改用 Ridge 自己的尺度。這個倍數來自目標函數的尺度換算,不能解讀成資料量增加後就應任意加強正則化。
GridSearchCV 會最大化 scoring,所以 neg_mean_squared_error 要取負號才是 MSE。若未指定 scoring,回歸估計器通常使用 $R^2$。$R^2$ 與 MSE 在單一固定測試集可作等價排序,但跨折平均時,每折的 TSS 不同,平均 $R^2$ 給的權重不同,未必選到相同 $\lambda$。Lab 畫兩種指標,就是要分清這個差別。
網格搜尋在訓練資料的內層折選 $\lambda$,再用整個訓練資料重擬合。best_score_ 是被選中的內層分數,不能當成獨立測試效能。Lab 先示範保留 25%、其餘 75% 內做五折選 $\lambda$,這是單次保留驗證加內部 CV;後面再用外層五折、內層五折的雙層迴圈,每一個外層折都重新選 $\lambda$,得到五個完整調參程序的保留分數。兩個版本都把選參數與最終評估分開,差別在外層只評估一次或輪流評估五折。
將 RidgeCV 或 ElasticNetCV 放在事先標準化整個 X 的後面,並不會自動讓標準化也在其內部每一折重新估計。要嚴格避免前處理看見保留折,使用對完整 Pipeline(StandardScaler, Ridge) 的 GridSearchCV。PCR 的 PCA 與 PLS 的方向估計也遵守同一原則。
查用法:ElasticNet 目標、GridSearchCV、RidgeCV。成分數或 $\lambda$ 選完後,最終模型仍要在允許使用的完整訓練資料重新擬合;獨立測試集只做最後評估。
第三條路:不動變數,改換座標。先把 $p$ 個相關的預測變數壓成 $M$ 個互不相關的方向($M\ll p$),再對這 $M$ 個方向做最小平方。
主成分迴歸(principal components regression, PCR)用的方向就是第 12 章的主成分: 第一主成分是 $X$ 變異最大的方向,第二個是在跟第一個正交的條件下變異最大的方向,依此類推。
$$Z_m = \sum_{j=1}^{p} \Phi_{jm} X_j, \qquad y_i = \theta_0 + \sum_{m=1}^{M} \theta_m z_{im} + \varepsilon_i$$$M$ 是調整參數,同樣用 CV 選。$X$ 滿欄秩且 $M=p$ 時 PCR 就退回最小平方(只是換了個座標); $M$ 小的時候,被丟掉的那些低變異方向就是被收縮掉的部分。
$Z_m=\sum_j\Phi_{jm}X_j$ 可以讀成投影:把第 $i$ 筆資料的 $p$ 個預測變數值 $(x_{i1},\ldots,x_{ip})$ 投影到方向 $(\Phi_{1m},\ldots,\Phi_{pm})$ 上,得到一個數,就是得分 $z_{im}=\sum_j\Phi_{jm}x_{ij}$。取 $M$ 個方向,就把每筆資料從 $p$ 個數換成 $M$ 個新特徵。PCA 的方向讓投影後的得分變異最大;這和「重建誤差最小」是同一件事,也和第 4 章 Fisher 判別的寫法形式類似,推導見本節末的收合區。
中心化之後,若方向滿足 $\sum_{j=1}^p\Phi_{jm}^2=1$,$z_{im}=\sum_j\Phi_{jm}x_{ij}$ 就是第 $i$ 筆資料沿該方向的正交投影座標:投影點的第 $j$ 個座標是 $z_{im}\Phi_{jm}$,剩下的殘差 $x_{ij}-z_{im}\Phi_{jm}$ 與方向垂直,因為
$$\sum_{j=1}^p\big(x_{ij}-z_{im}\Phi_{jm}\big)\Phi_{jm}=z_{im}-z_{im}\sum_j\Phi_{jm}^2=0 .$$方向長度不是 1 時,$z_{im}$ 只是一般的線性組合,數值會隨長度等比例縮放;這也是 PCA 要加上 $\sum_j\Phi_{jm}^2=1$ 的原因,否則把 $\Phi$ 放大就能讓變異無限變大。若 $M$ 個方向彼此正交($\sum_j\Phi_{jm}\Phi_{jl}=0$,$m\ne l$),每筆資料在這 $M$ 個方向上的投影點為 $\sum_{m=1}^Mz_{im}\Phi_{jm}$,等於把 $p$ 維資料投影到 $M$ 維子空間。
pca__n_components 這種雙底線寫法是 Pipeline 的參數命名慣例:步驟名 + __ + 該步驟的參數名。標準化同樣包在 pipeline 裡,PCA 必須在標準化之後做——不然變異數大的變數會主宰第一主成分。
Ch06-varselect-lab-zh.ipynb
array([0.3831424 , 0.21841076])
explained_variance_ratio_ 就是第 12 章的 PVE。累加起來可以回答「留幾個主成分才夠」。
Ch06-varselect-lab-zh.ipynb
什麼資料對 PCR 有利?講義用一組模擬資料示範:$X$ 的前五個主成分包含關於 $Y$ 的全部資訊。此時 PCR 的測試 MSE 在 $M=5$ 有明顯的最小值;Ridge 與 Lasso 也都比最小平方好,但 PCR 和 Ridge 略勝 Lasso。要注意,只有前五個主成分和 $Y$ 有關,不代表只有前五個原始變數和 $Y$ 有關。
講義與課本只說明反應完全由前五個主成分決定,沒有列出全部設定。下面用講義的符號寫出一種符合這個描述的生成方式:先產生主成分得分,讓 $Y$ 只依賴前 5 個得分,再把得分「轉回」原始預測變數。
第一步:產生主成分得分。對 $i=1,\ldots,n$、$m=1,\ldots,p$,獨立抽 $z_{im}\sim N(0,\sigma_m^2)$,並取 $\sigma_1^2\gt\sigma_2^2\gt\cdots\gt\sigma_p^2\gt0$,讓前面的成分變異較大。
第二步:用正交方向產生原始預測變數。選一組彼此正交的單位方向 $\Phi_{jm}$,滿足 $\sum_{j=1}^p\Phi_{jm}\Phi_{jl}=1$($m=l$)或 $0$($m\ne l$),令
$$x_{ij}=\sum_{m=1}^p z_{im}\Phi_{jm},\qquad j=1,\ldots,p .$$此時 $\operatorname{Cov}(X_j,X_k)=\sum_m\sigma_m^2\Phi_{jm}\Phi_{km}$,原始變數通常彼此相關;母體的主成分方向正是 $\Phi_{\cdot m}$,變異依序為 $\sigma_m^2$。正交性也讓得分可以反解回來:$z_{im}=\sum_j\Phi_{jm}x_{ij}$。
第三步:讓 $Y$ 只依賴前 5 個得分。
$$y_i=\theta_0+\sum_{m=1}^{5}\theta_m z_{im}+\varepsilon_i,\qquad \varepsilon_i\sim N(0,\sigma^2)\ \text{且與得分獨立}.$$前五個成分包含全部可預測的訊號,第 6 到第 $p$ 個成分與 $Y$ 沒有額外關係。代回 $z_{im}=\sum_j\Phi_{jm}x_{ij}$,就得到講義降維一節的係數關係:
$$y_i=\theta_0+\sum_{j=1}^p\beta_jx_{ij}+\varepsilon_i,\qquad \beta_j=\sum_{m=1}^{5}\theta_m\Phi_{jm}.$$$\Phi_{jm}$ 一般不為 0,所以 $\beta_1,\ldots,\beta_p$ 通常全部不為 0。主成分空間裡的稀疏,不等於原始變數空間裡的稀疏。
| 方法 | 在這個模擬中的情況 |
|---|---|
| PCR($M=5$) | 保留全部有效訊號,丟掉其餘 $p-5$ 個無關方向 |
| Ridge | 把變異小的方向縮得多、變異大的方向縮得少(見下方對照),效果接近 PCR |
| Lasso | 在原始變數上做選擇;真實的 $\beta_j$ 並不稀疏,它的長處用不上 |
有限樣本中,對 $X$ 算出的樣本主成分只會近似 $\Phi_{\cdot m}$,所以 PCR 用前五個樣本主成分時,只能近似地抓到全部訊號。若要在樣本中精確做到,可以先對已產生的 $x_{ij}$ 做 PCA,再用樣本得分代入第三步產生 $y_i$。反過來,若真實係數在原始變數上就稀疏(例如只有兩個 $\beta_j$ 不為 0),Lasso 才會佔優。
把資料改用主成分得分當座標後,兩種方法的擬合值都能寫成各主成分方向的貢獻相加,差別只在每個方向乘上的收縮因子:PCR 對前 $M$ 個方向乘 1、其餘乘 0;Ridge 對第 $m$ 個方向乘 $d_m^2/(d_m^2+\lambda)$,其中 $d_m^2=\sum_iz_{im}^2$ 隨 $m$ 遞減。變異大的方向幾乎不縮,變異小的方向縮得最多,Ridge 因此像是把 PCR 的「留或丟」換成平滑的漸進收縮(ESL 第 3.4.1、3.5.1 節)。
計算細節。
設 $x_{ij}$、$y_i$ 都已中心化,且 $n\gt p$、沒有完全共線。令 $\Phi_{jm}$($m=1,\ldots,p$)是全部 $p$ 個主成分方向,彼此正交且長度為 1;得分 $z_{im}=\sum_j\Phi_{jm}x_{ij}$。主成分得分兩兩不相關:$\sum_iz_{im}z_{il}=0$($m\ne l$),並記 $d_m^2=\sum_{i=1}^nz_{im}^2$,依 $m$ 遞減。$d_m^2$ 就是 $n$ 倍的第 $m$ 個特徵值 $\ell_m$(用 $\frac1n$ 的慣例,見本節末的收合區);模擬資料裡的母體變異 $\sigma_m^2$ 是它的母體對應。
換座標。任何線性預測都能改寫:令 $\theta_m=\sum_j\Phi_{jm}\beta_j$,因為方向正交,
$$\sum_{j=1}^p\beta_jx_{ij}=\sum_{m=1}^p\theta_mz_{im},\qquad \sum_{j=1}^p\beta_j^2=\sum_{m=1}^p\theta_m^2 .$$第二式表示旋轉不改變係數的長度,所以 Ridge 的懲罰在新座標下形式不變。
Ridge 拆成 $p$ 個一維問題。展開平方,交叉項因得分正交而消失:
$$\sum_i\Big(y_i-\sum_m\theta_mz_{im}\Big)^2+\lambda\sum_m\theta_m^2=\sum_iy_i^2+\sum_{m=1}^p\Big[(d_m^2+\lambda)\theta_m^2-2\theta_m\sum_iz_{im}y_i\Big].$$每個 $\theta_m$ 各自最小化方括號,得
$$\hat\theta_m^{\mathrm{ridge}}=\frac{\sum_iz_{im}y_i}{d_m^2+\lambda}=\frac{d_m^2}{d_m^2+\lambda}\,\hat\theta_m,\qquad \hat\theta_m=\frac{\sum_iz_{im}y_i}{d_m^2},$$其中 $\hat\theta_m$ 是 $y$ 對第 $m$ 個得分做簡單迴歸的係數,也就是 PCR 在該方向的係數。兩種方法的擬合值因此是
$$\hat y_i^{\mathrm{ridge}}=\sum_{m=1}^p\frac{d_m^2}{d_m^2+\lambda}\,\hat\theta_mz_{im},\qquad \hat y_i^{\mathrm{PCR}}=\sum_{m=1}^M\hat\theta_mz_{im}.$$PCR 的收縮因子是 1($m\le M$)或 0($m\gt M$);Ridge 的收縮因子 $d_m^2/(d_m^2+\lambda)$ 介於 0 與 1 之間、隨 $m$ 遞減。原始係數用 $\beta_j=\sum_m\Phi_{jm}\theta_m$ 換回。
有效自由度。由上式,中心化後的斜率部分 $\hat y_i^{\mathrm{ridge}}=\sum_k h_{ik}(\lambda)y_k$,$h_{ik}(\lambda)=\sum_m z_{im}z_{km}/(d_m^2+\lambda)$。對角元素加總得
$$\mathrm{df}(\lambda)=\sum_{i=1}^nh_{ii}(\lambda)=\sum_{m=1}^p\frac{d_m^2}{d_m^2+\lambda},$$從 $\lambda=0$ 時的 $p$ 平滑降到 $\lambda\to\infty$ 時的 0;PCR 用 $M$ 個方向時自由度就是 $M$。這就是前面 GCV 用到的 $\mathrm{df}(\lambda)$;加回截距那一份 $1/n$,才是留一捷徑裡的 $h_{ii}(\lambda)$。這套推導是 ESL 第 3.4.1 節(式 3.47)與第 3.5.1 節(圖 3.17)的純量寫法。
不算。每一個主成分都是全部 $p$ 個原始變數的線性組合,所以就算只留 $M=2$,最終模型仍可能需要全部原始欄位;個別 loading 可為零。
這是 PCR 跟 Lasso 最重要的差別:Lasso 給你「這 5 個變數有用、其餘丟掉」;PCR 給你「這 2 個方向有用」,而每個方向都攪拌了所有變數。後者在解釋上通常更難。你得去看 loadings 才知道那個方向大致代表什麼。
因為主成分是「變異最大的方向」,而變異數的大小完全取決於單位。
Credit 資料裡 Limit 的數值是幾千、Cards 是個位數。不標準化的話,第一主成分幾乎就等於 Limit 本身——這主要反映單位造成的尺度差異,不能直接當成重要性。第 12 章的 USArrests biplot 也呈現這個現象。
PCR 用 $M=2$ 個主成分。最終模型用到了幾個原始預測變數?
以下採每列一筆觀察的中心化矩陣 $X\in\mathbb R^{n\times p}$。若 $X=UDV^T$,則樣本共變異數 $S=X^TX/(n-1)=V\Lambda V^T$,其中 $\Lambda=D^2/(n-1)$(對角元素是特徵值 $\ell_j$,此處 $\Lambda$ 與懲罰參數 $\lambda$ 無關)。主成分方向是 $V$ 的欄,得分為 $XV=UD$;主成分方向滿足 $Sv=\ell v$;最大變異與重建誤差的關係見下方證明。講義也使用每欄一筆觀察的排列;對照公式時要連同矩陣一起轉置。
保留正特徵值方向後,白化得分 $Z=XV\Lambda^{-1/2}$ 的樣本共變異數為單位矩陣。ZCA 白化再乘上 $V^T$ 回到原特徵座標;保留全部正秩方向且滿秩時為 $XV\Lambda^{-1/2}V^T$。白化會把低變異方向放大,因此小特徵值可能放大雜訊;普通 PCR 不需要白化,也不能把白化當成保證改善預測的步驟。
若 $X=UDV^T$,則 $S=VD^2V^T/(n-1)$,特徵值為 $d_j^2/(n-1)$,得分為 $UD$。只保留正特徵值,令 $Z=XV_r\Lambda_r^{-1/2}$,其樣本共變異數為 $\Lambda_r^{-1/2}V_r^TSV_r\Lambda_r^{-1/2}=I_r$。ZCA 的回投影 $ZV_r^T$ 在 $p$ 維中的共變異數為 $V_rV_r^T$,只有 $r=p$ 時才是 $I_p$。
設 $x_{ij}$ 已中心化,樣本共變異數記為 $S_{jk}=\frac1n\sum_{i=1}^nx_{ij}x_{ik}$。講義附錄定義共變異數矩陣時用 $\frac1{n-1}$,推導最大變異時用 $\frac1n$;兩者只差常數倍,方向不變,本節用 $\frac1n$。方向 $(\Phi_{11},\ldots,\Phi_{p1})$ 上的得分 $z_{i1}=\sum_j\Phi_{j1}x_{ij}$,其樣本變異數除以方向長度的平方為
$$R(\Phi_{\cdot1})=\frac{\frac1n\sum_iz_{i1}^2}{\sum_j\Phi_{j1}^2}=\frac{\sum_{j=1}^p\sum_{k=1}^p\Phi_{j1}\Phi_{k1}S_{jk}}{\sum_{j=1}^p\Phi_{j1}^2}.$$這個比值稱為 Rayleigh quotient。它不受方向長度影響,所以「在 $\sum_j\Phi_{j1}^2=1$ 下最大化得分變異」等於最大化 $R(\Phi_{\cdot1})$。最大值是共變異數的最大特徵值,達到最大值的方向就是第一主成分方向;第二主成分在與第一方向正交的條件下再最大化,依此類推。
與第 4 章 Fisher 判別的連結。第 4 章的 Fisher 準則同樣在找投影方向 $a$,最大化的是「投影後的類間散布 ÷ 類內散布」。兩者分子分母都是「方向的二次式」,Fisher 的比值稱為廣義 Rayleigh quotient:分母換成類內散布,而不是方向的長度。PCA 不看類別,只問投影後整體散得多開;Fisher 問的是投影後類別分得多開。
最大變異等於最小重建誤差。用單位方向的投影點 $z_{i1}\Phi_{j1}$ 重建第 $i$ 筆資料,殘差與方向垂直(見本節前面的投影座標),所以每筆資料都滿足畢氏定理:
$$\sum_{j=1}^px_{ij}^2=z_{i1}^2+\sum_{j=1}^p\big(x_{ij}-z_{i1}\Phi_{j1}\big)^2 .$$對 $i$ 加總,左邊是固定的總平方和,與方向無關;右邊第一項是 $n$ 倍的得分變異,第二項是重建誤差。兩者相加為定值,所以讓得分變異最大的方向,必定讓重建誤差最小。取 $M$ 個正交方向時同理,總平方和等於 $M$ 個得分的平方和加上重建誤差。
為了簡潔,暫時把 $\Phi_{j1}$ 寫成 $\Phi_j$。在 $\sum_j\Phi_j^2=1$ 下最大化 $\sum_j\sum_k\Phi_j\Phi_kS_{jk}$(講義附錄的作法)。引入 Lagrange 乘數 $\ell$,對 $\sum_j\sum_k\Phi_j\Phi_kS_{jk}-\ell\big(\sum_j\Phi_j^2-1\big)$ 對 $\Phi_j$ 偏微分設為 0:
$$\sum_{k=1}^pS_{jk}\Phi_k=\ell\,\Phi_j,\qquad j=1,\ldots,p .$$所以極值方向必須是共變異數的特徵向量,$\ell$ 是對應的特徵值(這裡不用 $\lambda$,以免和懲罰參數混淆)。兩邊乘以 $\Phi_j$ 再對 $j$ 加總,得 $\sum_j\sum_k\Phi_j\Phi_kS_{jk}=\ell\sum_j\Phi_j^2=\ell$:在特徵向量上,得分變異恰好等於特徵值。因此最大值是最大特徵值 $\ell_1$,第一主成分就是對應的特徵向量。
要確認它是全域最大,把任一單位方向寫成正交特徵向量 $v_{\cdot1},\ldots,v_{\cdot p}$ 的組合 $\Phi_j=\sum_rc_rv_{jr}$,$\sum_rc_r^2=1$。代入後得分變異為 $\sum_r\ell_rc_r^2\le\ell_1\sum_rc_r^2=\ell_1$。
講義用城市人口 pop 與廣告花費 ad 說明 PCA。先減去每一欄平均,第一方向沿點雲最大的變動方向,第二方向與其垂直。每個點投影到第一方向後的座標是得分,原點到該點的距離不是得分。旋轉座標後,第一得分與原來兩欄通常都有較強關係,第二得分表示第一方向未保留的差異。
這個投影最小化的是各點到子空間的垂直平方距離,並不是把 ad 當反應、只最小化其縱向殘差的 OLS 直線。中心化決定方向要穿過平均點;縮放決定哪一種距離/變異受到重視。若原量尺本來具有共同科學意義,可以不標準化;若人口與費用量尺差異很大,標準化讓兩欄各以一個標準差衡量。PCA 不要求變數有相同物理單位,但要明確選定比較尺度。
本文用 $n\times p$ 的 $X$(每列一筆觀測),講義附錄用 $p\times n$ 的 $A=X^\top$(每欄一筆)。所以其 $S=AA^\top/(n-1)$ 對應本文的 $X^\top X/(n-1)$,講義得分 $U^\top A$ 對應本文 $XV$。閱讀 $U$、$V$ 時,要連同資料矩陣的排列一起看。
不是。對相同的保留方向做無懲罰 OLS,各得分的縮放會被迴歸係數抵消,擬合不變;白化的定義與限制見上方共變異數、SVD 與白化。
Lab 的 explained_variance_ratio_ 是 PCA 對 $X$ 的變異解釋比例;不能讀成 Salary 的解釋比例。少數 PC 捕捉大部分 $X$ 變異,仍可能漏掉對 $y$ 很重要的低變異方向。反例就是 PCR 節的前五個主成分模擬:那組資料刻意把反應資訊全放在前五個主成分,一般資料沒有這個保證。$M$ 的選擇應使用預測 CV,而非只以累積 $X$ 變異達到某百分比決定。
PCR 的原始係數限制在前 $M$ 個方向張成的子空間。它保留的是主成分方向,而沒有直接選出一組原始欄。Lab 將 StandardScaler、PCA 與 LinearRegression 放在同一 Pipeline 搜尋 $M$,讓每折重新學方向。只含截距的空模型可另行計算,不能把某個 API 的零成分是否可用當成空模型不存在的理由。
對照:PCA API 與得分/解釋變異、PCR 與 PLS 官方例子;完整非監督式解讀見非監督式學習章。
講義延伸來源:stats.stackexchange.com。
PCR 挑方向時完全沒看 $y$。若想讓 $y$ 也參與方向的選擇,可以使用 PLS。
偏最小平方(partial least squares, PLS)修掉這一點:第一個方向的權重直接取 $\Phi_{j1} \propto$ 每個 $X_j$ 與 $y$ 的簡單線性迴歸係數, 也就是跟 $y$ 相關性愈強的變數,權重愈大。 取完第一個方向後,把各變數對它迴歸取殘差,再在殘差上重複同樣的步驟得到第二個方向。
兩種方法的差別,寫成第一個方向的目標式最清楚。設各 $X_j$ 已標準化,方向長度固定為 1:
$$\text{PCR:}\ \max_{\sum_j\Phi_{j1}^2=1}\ \operatorname{Var}\Big(\sum_{j=1}^p\Phi_{j1}X_j\Big),\qquad \text{PLS:}\ \max_{\sum_j\Phi_{j1}^2=1}\ \operatorname{Corr}^2\Big(Y,\sum_{j=1}^p\Phi_{j1}X_j\Big)\operatorname{Var}\Big(\sum_{j=1}^p\Phi_{j1}X_j\Big).$$PCR 只要新特徵散得開;PLS 還要求它和 $Y$ 有關,化簡後等於找與 $Y$ 共變異平方最大的線性組合。上面講義用簡單迴歸係數當權重的作法,正好就是這個最佳化的解;推導與後續方向的算法見本節末的收合區。
上一節的互動示範了這個差別:拖滑桿改變 $y$ 的方向,PLS1 跟著轉、PC1 不動。
| 方向怎麼挑 | 有沒有用到 $y$ | 做變數選擇嗎 | 何時較好 | |
|---|---|---|---|---|
| PCR | $X$ 變異最大的方向 | 沒有 | 沒有 | $X$ 的主要變異方向剛好跟 $y$ 有關時 |
| PLS | 單位權重下與 $y$ 共變異最大的方向 | 有 | 沒有 | $X$ 的大變異方向跟 $y$ 無關時 |
| Ridge | 不換座標,全部收縮 | 有(擬合時) | 沒有 | 很多變數都有小貢獻 |
| Lasso | 不換座標,部分歸零 | 有(擬合時) | 有 | 真實係數稀疏 |
實務上 PLS 的表現常常跟 PCR 差不多,有時還略差。原因是: 它雖然降低了偏差(方向跟 $y$ 有關),但也增加了變異(方向是估出來的,而且用到了 $y$)。 ISLP §6.3.2 的結論就是這樣:PLS 沒有一致地贏過 PCR 或 Ridge。
PLSRegression 的 n_components 跟 PCR 的 n_components 是同一個角色的調整參數,一樣用 CV 選。注意它預設 scale=True,已經幫你標準化了。
Ch06-varselect-lab-zh.ipynb
PLS 的方向用到 $y$,所以成分數 $M$ 與標準化都要在訓練折內選定,否則驗證誤差會過度樂觀。
PLS 相對於 PCR 的差別是什麼?
設每個 $x_{ij}$ 已中心化並標準化($\frac1n\sum_ix_{ij}^2=1$),$y_i$ 已中心化。記 $c_j=\operatorname{Cov}(X_j,Y)=\frac1n\sum_ix_{ij}y_i$。為了簡潔,下面先把第一個方向的權重 $\Phi_{j1}$ 寫成 $\Phi_j$。新特徵 $Z=\sum_j\Phi_jX_j$ 與 $Y$ 的共變異是 $\operatorname{Cov}(Y,Z)=\sum_j\Phi_jc_j$,所以第一個 PLS 方向解
$$\max_{\sum_j\Phi_j^2=1}\Big(\sum_{j=1}^p\Phi_jc_j\Big)^2\quad\Longrightarrow\quad \Phi_{j1}=\frac{c_j}{\sqrt{\sum_kc_k^2}}\ \ (\text{正負號皆可}).$$另一方面,$Y$ 對單一 $X_j$ 做簡單迴歸的斜率是 $\operatorname{Cov}(X_j,Y)/\operatorname{Var}(X_j)$。標準化後 $\operatorname{Var}(X_j)=1$,斜率正好是 $c_j=\operatorname{Corr}(X_j,Y)\operatorname{SD}(Y)$。因此講義「用簡單迴歸係數當權重」與最佳化的解只差一個正規化常數,權重也與相關係數成正比。
注意目標是共變異而不是相關:$\operatorname{Cov}(Y,Z)=\operatorname{Corr}(Y,Z)\operatorname{SD}(Z)\operatorname{SD}(Y)$,PLS 同時考慮與 $Y$ 的關係及新特徵的變動量。若沒有先標準化,各斜率是 $\operatorname{Corr}(X_j,Y)\,\operatorname{SD}(Y)/\operatorname{SD}(X_j)$,不再與相關係數共用同一個倍數,所以講義的「權重與相關成比例」要連同標準化條件一起讀。
由 Cauchy–Schwarz 不等式,
$$\Big(\sum_{j=1}^p\Phi_jc_j\Big)^2\le\Big(\sum_{j=1}^p\Phi_j^2\Big)\Big(\sum_{j=1}^pc_j^2\Big)=\sum_{j=1}^pc_j^2 .$$只要不是所有 $c_j$ 都為 0,等號成立若且唯若 $\Phi_j$ 與 $c_j$ 成比例,再由 $\sum_j\Phi_j^2=1$ 決定比例常數,就得到上面的 $\Phi_{j1}$。簡單迴歸斜率 $\hat b_j=\sum_ix_{ij}y_i/\sum_ix_{ij}^2$;代入相關係數 $r_j=\sum_ix_{ij}y_i/\sqrt{\sum_ix_{ij}^2\sum_iy_i^2}$,即得 $\hat b_j=r_j\,\mathrm{SD}(Y)/\mathrm{SD}(X_j)$。
講義延伸來源:stats.stackexchange.com。
講義延伸來源:stats.stackexchange.com。
講義說後續方向是「取殘差後重複」。用純量寫出來(ESL 第 3.5.2 節的 PLS 演算法),每一輪都是「最大化共變異 → 從 $X$ 扣掉已找到的方向 → 重複」:
把每一輪的方向換回原始變數的權重 $\alpha_j$($Z_m=\sum_j\alpha_jX_j$,與上面作用在殘差上的 $\Phi_{jm}$ 不是同一組數字),第 $m$ 個方向等於在 $\sum_j\alpha_j^2=1$ 且「與前面的得分都正交」的條件下,最大化 $\operatorname{Cov}^2\big(Y,\sum_j\alpha_jX_j\big)$(ESL 第 3.5.2 節)。所以整個演算法就是把上面的共變異目標一輪一輪地解下去。
兩個細節。第一,殘差化後的 $x_{ij}^{(2)}$ 不再有單位變異數,此時 $Y$ 對它做簡單迴歸的斜率是 $\operatorname{Cov}/\operatorname{Var}$,不再等於共變異。要維持與最佳化的等價,權重應該用殘差與 $Y$ 的共變異(如上),或明確說明每一輪重新標準化並採用哪一種 PLS 版本。第二,有些演算法同時把 $y_i$ 換成殘差 $y_i-\hat y_i^{(m-1)}$。因為 $x_{ij}^{(m)}$ 與前面所有得分正交,而 $\hat y_i^{(m-1)}$ 是那些得分的線性組合,$\sum_ix_{ij}^{(m)}y_i=\sum_ix_{ij}^{(m)}\big(y_i-\hat y_i^{(m-1)}\big)$,兩種做法得到相同的方向。
預測新資料時,沿用訓練時估好的平均、尺度、各輪權重 $\Phi_{jm}$ 與殘差化係數 $a_{jm}$:新觀測依序算 $z_{\mathrm{new},m}=\sum_j\Phi_{jm}x_{\mathrm{new},j}^{(m)}$,再更新 $x_{\mathrm{new},j}^{(m+1)}=x_{\mathrm{new},j}^{(m)}-a_{jm}z_{\mathrm{new},m}$,不能用新資料重算方向。
PLS 的方向需要已觀測的訓練反應,PCR 的方向只看 $X$(第一方向的推導見PLS 的共變異目標)。
PLS 第一方向不等於直接求 OLS 的係數:OLS 同時解 $p$ 條正規方程,一次處理所有變數之間的相關;PLS 從各個 $\operatorname{Cov}(X_j,Y)$ 組成的方向開始,再逐步消去已解釋的成分。對完整秩方向使用足夠多成分,PLS 可重現 OLS 擬合;保留少數方向才施加低維限制。PLS 會形成新方向並逐步扣除已解釋的部分,不能只把各欄單變數回歸的預測值相加。
「看了 $y$」可避免 PCR 忽略低變異預測方向,但也可能把訓練雜訊當訊號。增加 $M$ 仍有偏差與變異取捨,所以 Lab 對 PLSRegression(n_components=M) 做 CV,而不是預設 $M$ 越大越好。其 scale=True 預設會在擬合時縮放 $X$ 與 $y$,若另加 StandardScaler 要分清每一步處理什麼;預測會回到 $y$ 的原量尺。
比較 PCR、PLS、Ridge 或 Lasso 時,使用同樣資料、同樣保留折與相同 MSE 定義。講義圖中的方向關係和 Lab 的 Hitters 最佳成分數,不構成某方法在所有資料上比較好的證據。
基因資料、文字資料、感測器資料常常是 $p$ 遠大於 $n$: 幾萬個變數、幾百個樣本。這時候傳統最小平方的唯一係數解與一般推論不再適用,需要額外的結構或選解規則。
當 $p\ge n$,且含截距的設計矩陣有 $n$ 個獨立的列時,最小平方能 完美擬合訓練資料:$\mathrm{RSS}=0$;$y$ 非常數時訓練 $R^2=1$。 這只說明模型能插值,不代表一定能或一定不能預測新資料。係數通常不唯一, 必須說明如何選解,再用獨立資料評估;不能拿訓練滿分當作泛化證據。
新增變數也會增加估計的不確定性。
每個新增變數都需要花自由度去估它的係數,而估計本身帶進變異。一個跟 $y$ 完全無關的變數,它的係數估計值是一個隨機的小數字,通常會偏離 0。那個隨機性直接進到預測裡。
若加入的變數沒有提供新的預測資訊,加入更多候選變數後,訓練誤差下降,測試誤差卻可能上升。因此要用獨立評估判斷新增變數是否值得保留。
這一章的三類工具全部都是為此設計的:子集/逐步選擇(明確減少變數)、收縮(Ridge 與 Lasso,Lasso 還順便選變數)、降維(PCR 與 PLS)。它們的共同點是降低變異,代價是接受一點偏差。
在 $p\gg n$ 的場景,Lasso 與 elastic net 特別受歡迎,因為稀疏假設常常合理(幾萬個基因裡有關的通常只有少數幾個),而且結果可以直接列出「這幾個變數」。
講義先用 $n=2$、截距加一個斜率示範:兩點可被直線完全通過,但沒有剩餘殘差自由度可估誤差變異。問題不是常態假設突然消失,而是訓練資料已不足以同時估模型與剩餘誤差。一般含截距模型的殘差自由度是 $n-\operatorname{rank}([1,X])$;若飽和且滿列秩,$\mathrm{RSS}=0$,通常使用 $\mathrm{RSS}/(n-p-1)$ 的估變異公式不再可用。
下一張圖固定 $n=20$,加入與反應沒有關係的特徵:訓練 $R^2$ 趨近一、訓練 MSE 趨近零,測試 MSE 卻變大。這不是新增資訊的功效,而是更多自由度能追蹤樣本中的偶然變動。訓練完美不能證明擷取到真實關係,低維與高維都需要獨立的預測評估。
講義的 Lasso 比較另固定 $n=100$,只有二十個真正有用特徵,再改變總 $p$ 與懲罰強度。即使真正訊號沒有增加,更多無關欄仍增加誤選、估計與調參的困難。這個例子說明,除了採用正則化,也必須選好 $\lambda$。「Lasso 可用於 $p>n$」並不表示「任何 $\lambda$ 都能克服高維」。
高維下,若多欄高度相關,很多不同係數組合可能給出相似預測。Lasso 挑出其中一欄,不等於其他欄沒有作用;Ridge 分散係數,不等於每欄皆有獨立證據。若選完變數後,直接套用事先固定模型的 p 值,就漏算了用同一份反應資料挑選欄位的影響,產生選擇後推論問題。可以先以外部資料驗證預測穩定性,但預測驗證本身不提供因果結論。
資訊準則若依賴全模型殘差估 $\sigma^2$,在這裡失去原本的計算前提;CV 則能比較可擬合的 Ridge、Lasso、PCR、PLS,但搜尋、縮放與方向估計仍必須在折內完成。評估模型選擇程序時,再用外層保留資料將選擇與報告分開。
下面幾題取自 ISLP §6.6 的課後習題,題號都對得回課本。先自己想過再點選項; 每個選項都有解說;需要完整解答時,可查閱下面3個站。
最佳子集、forward stepwise、backward stepwise 各自選出的「$k$ 個變數的最佳模型」中,測試誤差最低的會是哪一個?
相對於最小平方,Lasso 的特性是?
Ridge 的懲罰參數 $\lambda$ 從 0 開始往上增加。訓練 RSS 會怎麼變化?
課本第 9 題在 College 資料上比較最小平方、Ridge、Lasso、PCR、PLS 的測試誤差。預期會看到什麼?
考前可用這一頁複習各種方法的差別。
| 做什麼 | 調整參數 | 會歸零嗎 | 運算成本 | 何時最適合 | |
|---|---|---|---|---|---|
| 最佳子集 | 搜遍 $2^p$ 個模型 | 模型大小 $k$ | 會(沒選就是 0) | $2^p$,窮舉在 $p>40$ 左右不可行 | $p$ 很小 |
| Forward/backward | 貪婪地一次加/減一個 | 模型大小 $k$ | 會 | $1+p(p+1)/2$ | $p$ 中等 |
| Ridge(L2) | 全部係數一起收縮 | $\lambda$ | 通常不會;不是變數選擇方法 | 很快,有封閉解 | 很多變數都有小貢獻 |
| Lasso(L1) | 收縮 + 部分歸零 | $\lambda$ | 會 | 很快(座標下降) | 真實係數稀疏 |
| PCR | 換成 $M$ 個非監督式方向 | $M$ | 不選原始欄,個別 loading 可為零 | 很快 | $X$ 的大變異方向跟 $y$ 有關 |
| PLS | 換成 $M$ 個監督式方向 | $M$ | 不選原始欄,個別 loading 可為零 | 很快 | $X$ 的大變異方向跟 $y$ 無關 |
| 名稱 | 式子 | 備註 |
|---|---|---|
| $C_p$ | $\frac{1}{n}(\mathrm{RSS} + 2d\hat\sigma^2)$ | 式 6.2 |
| AIC | $-2\hat\ell+2k$ | 最大概似版本;共同固定誤差變異數時與 $C_p$ 型排序一致 |
| BIC | $\frac{1}{n}(\mathrm{RSS} + \log(n)\,d\hat\sigma^2)$ | 式 6.3,$\log n > 2$ 時罰得比 $C_p$ 重 |
| 調整後 $R^2$ | $1 - \frac{\mathrm{RSS}/(n-d-1)}{\mathrm{TSS}/(n-1)}$ | 式 6.4 |
| Ridge | $\mathrm{RSS} + \lambda\sum_j \beta_j^2$ | 式 6.5,對應限制式 $\sum_j\beta_j^2 \le s$(見 Lasso 節的限制式與懲罰式) |
| Lasso | $\mathrm{RSS} + \lambda\sum_j |\beta_j|$ | 式 6.7,對應限制式 $\sum_j|\beta_j| \le s$ |
| 主成分方向 | $Z_m = \sum_j \Phi_{jm} X_j$ | 式 6.16,$\Phi$ 由 $X$ 的變異決定 |
講義把兩個係數放在同一平面,以完整最小平方的 $(4,2)$ 作共同終點,讓不同複雜度的解形成路徑。最佳子集只能在坐標軸上的受限模型或完整模型之間選擇;Ridge 平滑地收縮,Lasso 可以沿軸停留到某係數進入;PCR 由保留幾個主成分決定可達子空間,PLS 的方向也使用 $y$。這張圖比較解的幾何,沒有表示所有方法都沿同一條直線縮回原點,也沒有提供測試風險的排名。
詞彙卡整理本章講義與 ISLP 第 6 章中常用的統計或機器學習專有名詞,正面附英文名稱。 先看正面的術語,試著說出定義,再翻面核對。可以洗牌複習,也可以搭配本章例子練習解釋。