劉應興的部落格

在線性迴歸模型 Yi = β_0 + β_1 X_1i + ... + β_k X_ki + ε_i，或矩陣式 Y = Xβ + ε 中，最常被使用的可能是最小平方估計 (least square estimate, LSE)

minimize Q(β) = (Y-Xβ)'(Y-Xβ)    或    (Y-Xβ)'W(Y-Xβ)

無加權的最小平方法適用於 Cov(Y) = σ^2 I，加權最小平方法則適用於 Cov(Y) = σ^2V，其中 V 已知且可逆，並取 W = V^{-1}，結果，Gauss-Markov 定理告訴我們：

［定理］可估函數 λ'β 之 LSE 是 LUE 中變異數最小的；
　　　　可估函數 Λ'β 之 LSE 是 LUE 中共變異矩陣最小的。

但是：當迴歸模型諸解釋變數 X_j 相互間有高度相關，使得 (X'X)^{-1} 其中有些數值巨大，或計算困難，甚至 (X'X) 不可逆，稱之為線性重合或多重共線性 (Multicollinearity) 問題，結果是 β 的統計誤差 ( Cov(β^) = σ^2(X'X)^{-1} ) 擴大，或計算誤差大（矩陣 (X'X) 是 ill-conditioned）, 甚或不能計算。雖然理論上要估計 E[Y] = Xβ 不成問題，但就統計應用上之解釋與預測，需要的是估計迴歸係數 β 而非 Xβ 。

再者，「不偏性 (unbiasedness)」雖然是很理所當然而被初學統計者認為是一個估計量應具備的「優良」性質，熟悉統計基礎觀念之後卻發現：要求估計量不偏實在不是一個聰明的選擇。依統計決策理論的想法，不偏性只是對統計決議（估計量或預測量、檢定函數、信賴集合）的一個附加限制；統計決策需要考慮的主要目標是在設定之損失函數 (loss function) 儘量降低其平均值（風險）。在線性迴歸問題上，假設 β 是我們觀注的重心，採用最常見的平方誤差和損失函數 (sum of squared error loss):

l(b, β) = (b - β)'W(b - β)

式中 b 是 β 的估計，若限制在線性估計量中考慮，可寫 b = AY，則

R(AY, β) = E[ (b - β)' W (b - β) ] 
        = E[(b - E[b])' W (b - E[b])] + (E[b] - β)' W (E[b] - β)
        = tr(W Cov(b)) +  (E[b] - β)' W (E[b] - β)
        = σ^2 tr(WAVA') + β' (AX-I)' W (AX-I) β

若有不偏的限制，則 (AX-I) β = 0 對所有 β 都要成立，則 AX = I，如最小平方估計是

A = (X'X)^{-1}X'，或 A = (X'wX)^{-1}X'w

式中 w 是任意可逆矩陣，都滿足 AX = I；事實上使 AY 是 β 的不偏估計的充要條件是 A 為 X 的左逆矩陣，AX = I。不過，這裡我們不想限制估計量不偏，則風險函數 (risk function) R(AY, β) 包含「變異」的成分 σ^2 tr(WAVA') 和「偏誤」的成分 β' (AX-I)' W (AX-I) β，我們希望允許些許偏誤，可以降低整體風險。

舉個簡單的例子，群體平均數 μ 常以樣本平均（暫以 X 表示）估計，在一些悄況下它甚至是「最優」的，但僅限於在不偏估計量之中。如果取 Y = aX + b, 則偏誤為 (a-1)μ+b, 而變異為 a^2σ^2/n，總均方誤差 (MSE) 或風險函數為

R(Y, μ) = MSE = a^2(σ^2/n) + [(a-1)μ + b]^2

例如 b = 0, a=0.9, 則 MSE = 0.81(σ^2/n) + 0.01μ^2，如果 μ 相對於 σ/√n 不大，則 Y 比 X 有更小的均方誤差。如果我們有額外的訊息知道 μ 很可能在 ν 附近，取 b = (1-a)ν，則 Y 等於在 X 和 ν 之間做個平均，結果

MSE = a^2(σ^2/n) + (1-a)^2(ν-μ)^2

當然，在這個例子，收縮估計量 (shrnkage estimator) Y 只在參數值處於一個小範圍時其 MSE 比標準估計量 X （樣本平均數）來得小。但如 Stein (1955) 發現的：如果同時估計 p 個常態群體平均數，可以找到比採用各樣本平均數當估計量（在平方誤差和損失下）一致更優的估計量，正是一種收縮估計量。在單一常態群體變異數的估計，也可以找到收縮估計量 aS^2，其均方誤差比樣本變異數 S^2 一致地較小，只要 (n-3)/(n+1) < a < 1, 例如 a = (n-1)/(n+1), 即計算變異數估計量時分母用 n+1 取代 n-1，是較準確的估計。

線性迴歸模型迴歸係數 β 有瘤如群體平均數，最小平方估計有如樣本平均數。在一般模型，在平方誤差損失下，LSE 是線性不偏估計量 (LUE) 中最好的；在常態誤差模型，它又是所有不偏估計量中最好的。但對有偏估計，特別是從線性估計來看，偏誤 (AX-I) β 是無界的，其（加權）平方和當然也是無界的。這造成兩個問題：我們必須在某個參數設定值之下，才能談論風險值最小化問題；或者依大中取小決策準則，必須將β 及 σ^2 限制在有界範圍，才可能極小化在參數範圍內之最高風險。一個簡單的想法：如同平均數收縮估計那樣，調整 (X'X)^{-1} 中的 X'X，同時解決 X'X 不可逆問題：

b = (X'X + λI)^{-1}X'Y + α

通常 α 取 0；不過，如果有訊息可判定 β 應是在特定點 γ 附近，不妨取

b = (X'X + λI)^{-1}X'Y + λ(X'X + λI)^{-1}γ

稱之為「脊估計量 (ridge estimator)」, 或說這是在做「脊迴歸 (ridge regression)」。調整 λ 值觀察風險函數值在參數值重要範圍內的變化，可以人為選擇「最適當的」λ 值。脊估計有時候也可用

b = (X'X + λＤ)^{-1}(X'Y + λＤγ)

式中 D 為適當的可逆對角線矩陣。如果用加權最小平方準則，

b = (X'wX + λＤ)^{-1}(X'wY + λＤγ)

通常 Cov(Y) = Cov(ε) = σ^2 V 時，w = V^{-1}。

在前項線性（迴歸）模型，假設誤差項是常態 N(0, σ^2 V), 又假設 β 服從 N(γ, Σ)，則得後驗分布 (posterior distribution)

β ～ N( (X'V^{-1}X/σ^2 + Σ^{-1})^{-1}(X'V^{-1}Y/σ^2 + Σ^{-1}γ)

因此，脊迴歸估計相當於貝氏估計 (Bayesian estimate)：λＤ 對應 σ^2 Σ^{-1}, w 對應 V^{-1}；或 λＤ 對應 Σ^{-1}, w 對應 V^{-1}/σ^2。

脊迴歸另一想法，在機器學習 (machine learning) 領域稱為正則化 (regularization，或譯正規化；但後一譯名也被用來表示 normalization，而這有時也被中譯規格化), 是將二次損失函數，即最小平方準則加上懲罰項

Q(β) = (Y-Xβ)'(Y-Xβ) + λ Σ_j |β_j|^p  或 (Y-Xβ)'(Y-Xβ) + λ ||β||_p^p

如果 λ 不是事先指定的，這就像限制範圍的最小平方法

minimization Q(β) = (Y-Xβ)'(Y-Xβ)

subject to:  ||β||_p = c

通常 p 取 1 或 2，分別稱為 L1 或 L2 正則化，其（在 ML 中）的作用是採取高階多項式迴歸模型時，藉由正則化或迴歸係數範圍限制避免過度配適 (over fitting) 問題。以 L2 正則化來說，其解正是

b = (X'X + λI)^{-1}X'Y

讓我們將 Q(β) 修改一下：

Q(β) = (Y-Xβ)'w(Y-Xβ) + λ(β-γ)'D(β-γ)

則其極小化之解在

b = (X'wX + λＤ)^{-1}(X'wY + λＤγ)

在統計決策理論的想法中,

R(AY, β) = E[ (b - β)' W (b - β) ]
        = σ^2 tr(WAVA') + β' (AX-I)' W (AX-I) β

權量矩陣 W 和前面最小平方脊迴歸法的 w 是兩回事，後者是對不同資料點 (Yi, x_i) 的加權，源於 Yi 或隨機誤差項 ε_i 的異幅變異數 (heteroscedasticity) 及相互間存在相關，故如能知道 Cov(Y) = V 或 σ^2 V 而 V 是已知矩陣，則取 w = V^{-1} 以獲得某種意義上的「最佳」。前者，W 矩陣則是對 β 各成分的重視程度有所不同，如果取 W 為對角線矩陣，損失函數

l(b, β) = (b - β)'W(b - β) = Σ_j Wj (b_j - β_j)^2

因諸 β_j 是隨變數 x(j) 的數值單位而變的，因此應依其數值單位及決策者之重視程度而決定權量。將風險函數用代數式展開：

     R(AY, β) = σ^2 Σ_iΣ_jΣ_kΣ_l w{ij}a{jk}v{kl}a{il}
                      + Σ_{i,j} w{ij} (Σ_kΣ_l a(ik}x{kl}β_l - β_i)(Σ_kΣ_l a(jk}x{kl}β_l - β_j)

在 σ^2, β 已知時，它是矩陣 A 的函數，對 A 之元素 a{rs} 的偏微分是

    D{rs} = σ^2 Σ_iΣ_jΣ_kΣ_l w{ij}v{kl}(δ(jr)δ(ks)a{il} + a{jk}δ(ir)δ(ls))
                   + Σ_{i,j} w{ij} (Σ_kΣ_l δ(ir}δ(ks)x{kl}β_l)(Σ_kΣ_l a(jk}x{kl}β_l - β_j)
                   + Σ_{i,j} w{ij} (Σ_kΣ_l a(ik}x{kl}β_l - β_i)(Σ_kΣ_l δ(jr)δ(ks)x{kl}β_l)
          = σ^2 (Σ_iΣ_l w{ir}v{sl}a{il} + Σ_jΣ_k w{rj}v{ks}a{jk})
                  + Σ_j w{rj} Σ_l x{sl}β_l (Σ_kΣ_l a(jk}x{kl}β_l - β_j)
                   + Σ_i w{ir} Σ_l x{sl}β_l (Σ_kΣ_l a(ik}x{kl}β_l - β_i)

基於權量矩陣 W 和（誤差項）共變異矩陣 V 是對稱矩陣，上列結果簡化為

D{rs} = 2σ2 Σ_iΣ_l w{ri}a{il}v{ls}
               + 2 Σ_i w{ri} Σ_l x{sl}β_l (Σ_kΣ_l a(ik}x{kl}β_l - β_i)

第一項是 σ^2 WAV 的第 (r,s) 元素，第二項是 W (AX - I) ββ'X' 的 (r,s) 元素。故平穩點（偏導數均為 0）為

σ^2 ＷAV + W (AX - I) β β'X' = 0

假設 W 是可逆的，則可以從上列關於 A 的方程式中移除，則

A = β β'X (σ^2 V + Xββ'X')^{-1}

此式是在特定參數值 σ^2, β 之下計算，但 β 未知，而 A 就是用來估計 β 的；而且要計算 n×n 方陣的反矩陣，比 X'X 的階數大多了。

上列結果與權量 W 無關，何故？如果我們的損失函數改為

l(b, β) = (b - λ'β)^2

一個特殊惰形是只看某一個迴歸係數。假設 b 限制在 a'Y 形式的線性估計，則

R(a'Y, β) = E[(a'Y - E[a'Y])^2] + (E[a'Y] - λ'β)^2
        = σ^2 a'Va +  (a'Xβ - λ'β)^2

這是 a 的凸函數，其最低點發生在平穩點，即對 a 的偏導數（梯度）為 0 的地方

２σ^2 Va + ２(Xβ)(β'X'a - β'λ) = 0

其解為 a =  (σ^2 V + Xββ'X')^{-1}Xββ'λ。此結果套用在 β 的各成分，與前面整體考慮 β 的估計比較，可知欲整體 R(AY, β) 最小，就是使 β 各成分的 MSE 最小，因此 W 之加權不重要。

當 A = β β'X'(σ^2 V + Xββ'X')^{-1} 時，

(AX - I) β = β β'X (σ^2 V + Xββ'X')^{-1}Xβ - β

令 β'X (σ^2 V + Xββ'X')^{-1}Xβ = θ，則 R(AY, β) 中屬於偏誤的部分是 (θ-1)^2 β'Wβ。另方面，

AVA' = β β'X (σ^2 V + Xββ'X')^{-1} V (σ^2 V + Xββ'X')^{-1}X' β β'

把 β'X (σ^2 V + Xββ'X')^{-1} V (σ^2 V + Xββ'X')^{-1}X' β 用 φ 表示，則

tr(WAVA') = tr(Wφββ') = φ β'Wβ

因此，R(AY, β) = [φσ^2 + (θ-1)^2]β'Wβ，這算是在參數值為 (σ^2, β) 時風險函數的下限，而在有界參數值範圍，可得此下限之最大值。將這統計決策問題看成是兩人零和賽局 (zero-sum game)，此最大值即為此賽局的小中取大解 (maximin solution)。另方面，在設定參數值範圍，也可對每一線性估計量 AY 得其風險函數最大值，而後變動 A 得 β 之大中取小估計量 (minimax estimator)。一般而言大中取小解的值不小於小中取大值，適當條件下以上兩結果會相等，也就是說有一個均衡解。