劉應興的部落格

假設 Y 是一個置中的（centralized, 減去平均數或即期望值）反應變數，X 是一個「解釋變數」，假設也是置中的變數。不論 Y 與 X 間是何種關聯，我們可以假設 Y = βX + ε, 這就是標題所謂的「擬線性關係」。這裡 ε 不應稱為隨機誤差，比較適合稱為「殘差」；不過因產差一詞常用於以樣本估計一個群體模型後其隨機誤差項的估計，這裡如有必要，就暫且稱 ε 為「殘餘」或殘餘項。用在隨機向量時，如主成分分析就在取一各行為相互正交之單位向量構成的矩陣 Q 使 tr(Cov(Q'X))/tr(Cov(X)) 極大化，也就是 tr(Cov(X-QQ'X)) 極小化，或說是極小化殘餘項總變異數。

當 Y, X 都是實數值隨機變數時，最小平方法要求 Var(Y-βx) 最小。由於

Var(Y-βX) = Var(Y) - 2β Cov(Y,X) + β^2 Var(X)

可得 β 的最小平方解是

β = Cov(Y,X)/Var(X)

此時 ε = Y-βX 與 X 是零共變異 Cov(X,ε) = 0，零相關 Corr(X,ε) = 0；而 corr(βX,X) = 1, Cov(βX,X) = β Var(X) = Cov(Y,X)。從殘餘項 ε 與解釋變項 X 的零相關，我們可以說：Y 與 X 的擬線性關係，就是把 Y 分解為 X 方向和與 X 「垂直」方向的和，所以 βX 可以說就是 Y 在 X 方向的垂直投影。總結來說，Y 與 X 的擬線性關係可以從兩點來看：一是殘餘項變異數最小化；二是 Y 在 X 與其垂直方向的分解或即垂直投影。統計上簡單直線迴歸的估計即是上列結果的樣本版本。

從統計上或實務上來說，Y 可能與 X 並不具備有意義的關聯，因為「變數」其實是個表象，可能只是某個不可觀測變數的代理。這樣的觀念其實就是因子模型的單變量版本，不過因為代理的可觀測變數與不可觀測的底層變數之間是一一對應，而不像因子模型假設的多個變數共同表現少數的共同因子，因此稱之為「量測誤差模型 (measurement error model, 或 errors-in-variables model)」——當然，實際上量測誤差模型是包含反應變數：

X = U + δ
Y = βU + ε

反應變數 Y 不和 X 直接關聯，而是和 X 所代理的 U 有關聯。雖然 Y 依舊可以寫成和 X 呈線性關聯的樣子：

Y = β(X-δ) + ε = βX + (ε-βδ)

但是，

Cov(X,ε-βδ) = Cov(U,ε) + Cov(δ,ε) - β Cov(U,δ) - β Var(δ)

右邊前三項基本上假設是零相關甚至機率獨立的，那麼寫成 Y = βX + ε*, 則 ε* 不再和 X 保持正交（垂直）性；這也就是說：如果我們計算 Y 與 X 的擬線性關係，在 X 方向的係數不再是 β。因此，量測誤差模型必須有自己一套求解或估計方法。在 U, δ, ε 兩兩零相關的條件下，

Var(X) = Var(U) + Var(δ)
Var(Y) = β^2 Var(U) + Var(ε)
Cov(X,Y) = β Var(U)

有 4 個未知數：β, Var(U), Var(δ), Var(ε)，只有 3 個方程式，必須再加一條件才能計算所有未知數，或者是設定的 4 個未知數間的關係，或者考慮高階動差。注意 Y-βU = ε, 不能如前面 Y 對 X 找擬線性關係一般以極小化 Var(ε) 的方式來求得 β 的解。

假設反應變數 Y 仍是單變量，解釋變數 X 是多變量（隨機向量），則最小平方法要極小化

Var(Y-X'β) = Var(Y) - 2β' Cov(X,Y) + β' Cov(X) β

式中 β' 是列向量，是行向量 β 經矩陣轉置 (transpose) 而得；同樣列向量 X' 是行向量 X 的轉置。最小平方的解是

β = (Cov(X))^(-1) Cov(X,Y)

式中 Cov(X) 產生的線性空間是 X 所在的空間，而 Y = X'β + ε 同樣將 Y 分解為相互垂直的 X'β 和 ε 兩部分，X'β 就是 Y 在 X 伸展的空間上的垂直投影。

Cov(Y-X'β,X'β) = β' Cov(X,Y) - β' Cov(X) β = 0；

而且 X'β 在 X 伸展的空間上，而最小平方 || Y-X'β ||^2 的意思就是 Y 到 X 空間上一成員 X'β 距離最短，這就是垂直投影。

我們可以換個想法來看 β 的選取：如果我們選擇 β 是要使 Y 與 X'β 相關係數 

Corr(Y,X'β) = β'Cov(X,Y)/√(Var(Y) Var(X'β))

最大。由於相關係數是尺度不變的，也就是 β 沿同一方向漲縮一個倍數，上列相關係數不會改變，因此，問題可以改成：

在 Var(X'β) = c (某個常數）下極大化 Cov(Y,X'β)。

適當選取 c 可以得最適解 β = (Cov(X))^(-1) Cov(X,Y)。把原本 f(x)/g(x) 的極值問題化成在 g(x) = c 限制下 f(x) 的極值問題可行嗎？

Suppose f(cx) = c f(x) and g(cx) = c g(x) for any constant c.
Let x* be the optimal solution for problem
    max f(x),  s.t. g(x) = c .
Then x* is also the optimal solution for
    max f(x)/g(x).
[proof]
For any x ≠ x*, let c* = g(x), then
    g((c/c*)x) = c = g(x*), 
so f((c/c*)x) ≦ f(x*).
Therefore,,
    f(x)/g(x) = f((c/c*)x)/g((c/c*)x)
                  = f((c/c*)x)/g(x*)
                  ≦ f(x*)/g(x*)

因此看來把極大化相關係數改成適當限制條件之下極大化共變異數是可行的。

如果目標就是極大化 Cov(Y,X'β)，由於 

Cov(Y, cX'β) = c Cov(Y,X'β)

所以也需要一個限制式來限定 β 的大小，則前述 Var(X'β) = c 型的限制式得相同結果。問題是：我們的目標如果單純是共變異數最大化，對 β 規模的限制方法很多，一定要用 Var(X'β) = c 來限制嗎？如果限制條件改成 β'Wβ = c, 則 β = c' W^(-1) Cov(X,y),  其中 c' 是由 c 所決定的另一常數。因此，極大化 Cov(Y, X'β) 欲得到 β = (Cov(X))^(-1) Cov(X,Y) 的解，對 β 規模之限制條件必須採 Var(X'β) = c 的形式，也就是 β' Cov(X) β = c, 而具體的 c 選擇是

c = Cov(Y,X') (Cov(X))^(-1) Cov(Y,X)

也就是說，限制條件為

β' Cov(X) β = Cov(Y,X') (Cov(X))^(-1) Cov(Y,X)

假設 Y 有 q 個變數，X 有 p 個變數，我們希望分別做線性組合成單一變數後可以極大化相關係數：

max ρ = Corr(α'Y,β'X)
       = α' Cov(Y,X) β/√[α' Cov(Y) α)(β' Cov(X) β)]

為了簡化，我們取 Cov(Y), Cov(X) 的「平方根」，並令

a = (Cov(Y))^(1/2) α
b = (Cov(X))^(1/2) β

則

ρ = a' (Cov(Y))^(-1/2) Cov(Y,X) (Cov(X))^(-1/2) b/(||a|| ||b||)
   ≦ ||(Cov(Y))^(-1/2) Cov(Y,X) (Cov(X))^(-1/2) b||/||b||

不等式是應用 Cauchy 不等式，等號成立的條件是 (Cov(Y))^(-1/2) Cov(Y,X) (Cov(X))^(-1/2) b 與 a 同向；同理可以得不等式

ρ ≦ ||(Cov(X))^(-1/2) Cov(X,Y) (Cov(Y))^(-1/2) a||/||a||

等號成立的條件是 Cov(X))^(-1/2) Cov(X,Y) (Cov(Y))^(-1/2) a 與 b 同向。 令

W = (Cov(Y))^(-1/2) Cov(Y,X) (Cov(X))^(-1/2)

則 W' = Cov(X))^(-1/2) Cov(X,Y) (Cov(Y))^(-1/2)，前述 Cauchy 不等式即 a'Wb ≦ ||a||．||Wb|| 及 ||W'a||．||b||；而兩等號成立之條件就是存在 s > 0, t > 0, 使

Wb = sa, W'a = tb;  故  WW'a = sta, W'Wb = stb

故 a, b 分別是非負確定矩陣 WW', W'W 的特徵向量，且對應相同之特徵值；而由 a 可得 b 之方向，由 b 亦可得 a 之方向。不失一般性，取 ||a|| = 1 = ||b||, 則可得所要的最大相關係數 ρ。同時，由 a, b 可計算 α, β；再由 α, β 可得 V1 = α'Y, U1 = β'X, 而 ρ = Corr(U1,V1)。

由上面的敘述知 a, b 固然是不同矩陣的特徵向量，卻對應相同的特徵值。實際上我們可得對應矩陣的所有正特徵值 λ(i) 及對應的、相互正交的特徵向量 a(i), b(i)，由之所得到的組合變量 Ui, Vi 滿足

Cov(Ui,Uj) = δij = 1 if i = j;  = 0 if i ≠ j
Cov(Vi,Vj) = δij

由於 Wb(j) = s(j)a(j), 故

Cov(Ui,Vj) = ((Cov(Y))^(-1/2) a(i))'Cov(Y,X)(Cov(X))^(-1/2) b(j) 
           = a(i)' (Cov(Y))^(-1/2) Cov(Y,X) (Cov(X))^(-1/2) b(j) 
           = a(i)'Wb(j) = a(i)'s(j)a(j) = s(j)a(i)'a(j) = s(j)δij

因此，對應 WW'（或 W'W）建立之新變量 Ui, Vi 滿足：除對應相同正特徵值的成對 (Ui,Vi) 具有正相關係數值外，其餘配對均零相關。這些 Corr(Ui,Vi) 稱為 Y, X 兩組變數間的「典型相關係數 (canonical correlation coefficients)」；建構之新變量 Ui, Vi 等，稱為「典型變量 (canonical variate)」；利用樣本資料估計典型相關，並用以解釋資料所代表的實務問題，就稱為典型相關分析。以群體模型來說，WW' 或 W'W 的正特徵值個數等於

rank(WW') = rank(W) = rank(Cov(Y,X))

單一 Y 變數對一或多變數的 X 的複相關係數如果非 0，就是其唯一的典型相關係數。

在典型相關問題中，Y 和 X 的角色可以說是對等的；如果 Y 是反應變數，X 是解釋變數，我們比較喜歡或習慣用迴歸模型來表示或處理 Y 與 X 間的關係。前面「擬線性關係」可說是對 Y 與 X 強加的關係，比較偏向對 Y 的分解；而迴歸模型則是假設 Y 與 X 具有真正的統計關係，例如 Y = f(X,ε), Y = f(X) + ε, Y = X'β + ε 等，其中 ε 是隨機誤差項，通常假設它和解釋變項是獨立的。其中 Y = X'β + ε 稱為線性模型 (linear model) 或線性迴歸模型 (linear regression model)，這裡「線性」一詞不是對 X，而是對 β，因為其中 X 隨機向量的內容或許包含同一隨機變數的多個函數，例如 Y = β1*X + β2*X^2 + ε, 解釋變數其實只有一個，X；但解釋項卻有兩個，X 和 X^2。這模型對 Y 與 X 的關係是曲線，或非線性的；但對 β = (β1,β2)' 而言是線性的。對線性模型而言，最小平方法在 Cov(X) 可逆的情況下可得 β = (Cov(X))^(-1) Cov(X,Y)。樣本版的線性模型允許 X 是非隨機（已知）的，例如實驗設計模型；更允許誤差項在不同樣本點之間不同及不獨立，此時會改用加權最小平方法或一般化最小平方法，以得到 β 更好的估計。關於樣本的迴歸分析及推論這裡不細談，我們只關心群體模型，藉群體模型來了解這類統計方法的本質。

在 Y 只是一個（實數值）隨機變數時，誤差項 ε 也是單一隨機變數，因此我們只需要求 ε 與 X 獨立或僅僅零相關，最小平方法很合適。若 Y 是隨機向量呢？此時模型為

Y = B'X + ε 或 Y' = X'B + ε'

其中 Y 和 ε 是 q×1 隨機（行）向量，B 是 p×q 迴歸係數矩陣，X 是 p×1 隨機向量。上列模型右式的樣本版是

Y(n×q) = X(n×p)B(p×q) + ε(n×q)

括號中是資料或係數矩陣的大小，其中 n 是樣本大小，q 是反應變數個數，p 是解釋變數個數。

回到群體模型，仍假設 X, Y 均置中。如果 Cov(ε) = σ^2 Ω, 則

σ^2 Ω = Cov(ε) = Cov(Y-B'X)
  = B' Cov(X) B - Cov(Y,X) B - B' Cov(X,Y) + Cov(Y)

取  q(B) = B' (Cov(X))^(1/2) - Cov(Y,X) (Cov(X))^(-1/2), 則

σ^2 Ω = (q(B)) (q(B))' - Cov(Y,X) (Cov(X))^(-1) Cov(X,Y) + Cov(Y)

要「極小化」Cov(Y-B'X)，顯然需要 q(B) = 0, 也就是 B = (Cov(X))^(-1) Cov(X,Y)。

但是，如果對 ε 做一個線性變換，ε* = Ω^(-1/2) ε, 則得 Cov(ε*) = σ^2 I。這相當於

Y* = Ω^(-1/2)Y = Ω^(-1/2)B'X + ε*, Cov(ε*) = σ^2 I

此新模型似乎應得

ＢΩ^(-1/2) = (Cov(X))^(-1) Cov(X,Y) Ω^(-1/2)

所以 B 仍是 (Cov(X))^(-1) Cov(X,Y)。

以上最小平方法考慮「極小化」Cov(ε)；為何不考慮 E[ε'Ｗε]？此處假設 W 為正半確定 (positive semidefinite) 對稱方陣，則

E[ε'Wε] = E[(Y-B'X)'W(Y-B'X)] = E[tr(W(Y-B'X)(Y'-X'B))]
   = E[Y'WY - X'BWY - Y'WB'X + X'BWB'X]
   = E[tr(WYY'-WYX'B-WB'XY'+WB'XX'B)]
   = E[tr(W(YY'-YX'B-B'XY'+B'XX'B))]
   = tr(W(Cov(Y)-Cov(Y,X)B-B'Cov(X,Y)+B'Cov(X)B))

我們發現用此「期望二次式誤差」 結果只是把原來的誤差共變異矩陣左乘權量矩陣 W 然後取對角線和，因此 B 的最小平方解仍是一個「最適解」。

以上我們一直假設 Cov(X) 的反矩陣存在，或稱它是可逆的 (invertible)，或限定 Cov(X) 是「滿秩的 (full rank)」。如果 Cov(X) 不滿秩，也就不可逆，那麼前述 β 或 B 的公式就不適用了。設 Cov(X) = ω, 則 P(X in C(ω)) = 1, 意謂 X 就分布於 ω 各行所張開的空間之內。舉個例子，Z 是任意一個非退化的隨機變數，E[Z] = 0, Var(Z) = v。令 X = [Z Z]', 則 ω = Cov(X) 就是元素皆為 v 的 2×2 方陣，C(ω) 即 R^2 上 45° 線，X 就分布在這線上。如果 [W Z]' 是非退化隨機向量，X = [ W Z W+Z]', 則 ω = Cov(X) 的最後一行是前兩行之和，最後一列是前兩列之和。所以 C(ω) 是 WZU 空間中的平面 U = W + Z, 而 X 就分布在此平面上。對於此類 X，結合反應變數 Y 後 XY' 是什麼呢？XY' 各行即是 X 的不同倍數，也就是在 X 所在的空間上，換句話說，XY' 各行都在 C(ω) 上，因此 E[XY'] 各行也在 C(ω) 上。

雖然 ω 不可逆，但我們可取 ω 的一個「廣義反矩陣 (generalized invers)」υ = ω^- 滿足 ωvω = ω; 由於 ω 是對稱的、正半定的方陣，υ 也可取對稱的。這是可以做到的，設 A 是任意非 0 對稱正半定方陣，則可分解為 A = BB', 其中 B 是滿行秩，也就是 B (及 A) 的秩等於 B 的行數，令 Ｈ = (B'B)^(-1)B', 則 HB = I。取 K = A^(-) = H'H, 則 K 是對稱方陣，且

AKA = BB'H'HBB' = B(HB)'(HB)B' = BB' = A

也就是說 K 是 A 的一個對稱的廣義反矩陣。若 w 在 C(ω) 中, 則存在 b 使 w = ωb, 因此

ωυw = ωυωb = ωb = w

因此，ω = Cov(X) 可以寫成 ωυω, 而 Cov(X,Y) 可以寫成 ωυ Cov(X,Y)。於是，最小平方準則 Cov(εε') 可以寫成

Q(B) = B' ωυω B - Cov(Y,X) υω B - B' ωυCov(X,Y) + Cov(Y)
        = (B'ω-Cov(Y,X))υ(ωB-Cov(X,Y)) 
                          - Cov(Y,X) υ Cov(X,Y) + Cov(Y)

由於 υ 也是正半定，顯然 ωB = Cov(X,Y) 可以極小化 Q(B)，而 B = ((Cov(X))^-)Cov(X,Y) 是一最適解。由於廣義反矩陣不唯一，這最適解也是不唯一，在樣本版這問題就是「不可估的 (inestimable)」。但由於 X 在 C(ω) 中 (with probability 1), B'X 不受 υ 之不同選擇的影響，在樣本版本，它是一個可估函數。

雖然理論上即使 Cov(X) 不可逆，最小平方法也能適用；但實務上莫說非滿秩情況，Cov(X) 或其樣本版本的 X'X 存在特徵值偏低情形時，都會使計算發生困難，至少導致計算誤差偏高，在樣本版本估計量之理論誤差也偏高。因此，統計學者有了許多不同方式來解決這樣的問題。首先是 X 的正交化，這改變了（擬）線性式 B'X 的外形 B*'X*, 其中 X* = AX，新變數間具有正交性。這種方法特別是在多項式迴歸時常用，因為 x, x^2, x^3,... 之間會有高相關，改成正交多項式 p1(x), p2(x), p3(x),... 在迴歸係數 β 或 B 之計算上可以避免一些問題。主成分迴歸也是將 X 用相互正交的 X* 取代，是 X 的成員數量龐大或彼此高度相關時適用的方法

Y = B'X + E = B*'X* + E*

其中 X* 是 X 的前少數主成分，而 X* 之外的被併到誤差項。為了方便解釋，以 X* 表現的迴歸式可以再轉回用 X 表示，等於允許 B 有少量的偏誤換取計算上或估計上的方便和穩定。不過，這兩種方法都只考慮對 X 做正交化轉換以便計算或估計，卻未從 Y 與 X 的關係上考慮。

偏最小平方法或部分最小平方法 (partial least squares, 簡寫 PLS)，通過投影解釋變數 X 和反應變數 Y 到一個新空間來尋找一個線性迴歸模型。仍假設 X, Y 已置中，依維基上的說明，PLS 迴歸底層假設

Y = Q V + E
X = P U + D

而目標是 V, U 之間共變異量的極大化。最後可從 V 對 U 之迴歸關係可反推 Y 對 X 的迴歸式：

V = B'U + R  ==>  Y = Q[B'P'(X-D)+R]+E

從最後 V 和 U 的迴歸關係來看，有點量測誤差模型的味道，或說是量測誤差模型的多變量版本；而 Y 和 X 的分解又有因子模型的味道，但比因子模型多了要求 P, Q 各行正交（且單位化）的要求（参見中、英文維基）；如果分別取 V, U 是 Y, X 的主成分，則 PLS 又和主成分關聯上；而 PLS 要求 V 與 U 之極大共變異，卻又和典型相關扯上關係。但以上個人理解的正確嗎？PLS 究竟要如何進行？