劉應興的部落格

先前在「細談：樣本空間與隨機變數」一文我們提及除了實數值和複數值隨機變數之外，可以有向量值隨機變數（隨機向量）、矩陣值隨機變數（隨機矩陣）等等，統計上就有一個重要的機率分布：Wishart 分布，就是涉及一種隨機矩陣，scatter matrix: 設 X1,...,Xr 是 i.i.d. N(0,V) 的 R^p 隨機向量，則隨機矩陣 S = Σ_i XiXi' 的分布稱為具自由度 r 的 Wishart 分布，以 W_p(V,r) 或 W(V,p,r) 表示。

設 V 是正確定的 (positive definite), r ≧p（英文版 wiki）或 r > p（中文版 wiki），則 S 有機率密度

    f(S) = |S|^{(r-p-1)/2}e^{-tr(V^(-1)S)/2}/{2^(rp/2)|V|^(r/2)Γ_p(r/2)}

其中

            Γ_p(x) = ∫_{S > 0} e^{-tr(S)}|S|^{x-(p+1)/2} dS

稱為多變量 gamma 函數（其實只有單變量 x，但定義式涉及矩陣變量 S），式中 S > 0 表示 S 是正確定矩陣。多變量 gamma 函數和普通 gamma 函數有下列關係：

Γ_p(x) = π^{p(p-1)/4} Π_i Γ(x+(1-j)/2)

如 p = 2 時 Γ_2(x) = √π Γ(x) Γ(x-1/2)；又可得遞迴關係 Γ_{p+1}(x) = π^(p/2) Γ(x-p/2) Γ_p(x)。

隨機矩陣 S 是對稱 p 階矩陣，有 p(p+1)/2 個變量；隨機向量 Xi 值域在 R^p。因此，要由 X1,...,Xr 的 p.d.f. 導得 S 的 p.d.f., 需要 rp ≧ p(p+1)/2, 即 r ≧ (p+1)/2。但這並不夠，以 p=3, r=2, V=I 為例，兩個 X 向量等於 6 個實變量，矩陣變量 S 也有有 6 個（數學）獨立變量（數學上的獨立只是表示其中任一個變量不能用其他變量的函數表示）：Sij, i≦j，按說兩者之間可以相互求解；但經計算發現 S 變量對諸 X 的 Jacobian 行列式值為 0, S 矩陣的行列式值亦為 0。事實上 S = Σ_i XiXi' 而 Xi 為 p×1 行向量，則需 r≧p 才能得到正確定的 S，所以 S 存在 p.d.f. 的條件是 r≧p。

設 S = Σ XiXi' 其中 Xi, i=1,...,r, i.i.d. N(0,V); A 是全列秩 (full row rank) m×p 矩陣。則 ASA' = Σ(AXi)(AXi)', AXi, i=1,...,r, i.i.d. 服從 N(0,AVA')。因此，若 S 服從 W(V,p,r), 則 ASA' 服從 W(AVA',m,r)。因 V 是正確定，存在一個可逆的矩陣 R 使 V = RR', 可稱 R 是 V 的一個「平方根」。令 A = R^(-1), 則 AVA = I, 故 Wishart 可以「標準化」成 W(I,p,r), 其 p.d.f. 是

f(S) = |S|^{(r-p-1)/2}e^{-tr(S)/2}/{2^(rp/2)Γ_p(r/2)

它是 r 個標準化 p 階常態隨機向量自乘和的分布，即 S = Σ ZiZi' 的分布。當 p = 1 時，W(I,1,r) 即 χ^2(r)。而對應任意正確定 V，W(V,p,r) 也可以由 W(I,p,r) 得到。

首先，我們來看 p = 2 的情形。S 的（數學）獨立元素有 S11, S22 和 S12。由於 V = I, 給定 S11=u 時 S12 與 S22* = S22-S12^2/S11 相互獨立，前者服從 N(0,u) 而後者服從 χ^2(r-1)， 所以 S22*=v*, S12=w 的條件 j.p.d.f. 是

g(v*,t|u) = v*^{(r-3)/2}e^(-v*/2)/{2^((r-1)/2)Γ((r-1)/2)}
                          ．e^{-w^2/(2u)}/√(2πu)

轉回 S22 = S22*+S12^2/S11, 則得給定 S11=u 時 S22=v, S12=w 的條件 j.p.d.f. 是

h(v,w|u) = g(v,w|u)|J| 
   = (v-w^2/u)^{(r-3)/2}e^{-(v-w^2/u)/2}/{2^((r-1)/2)Γ((r-1)/2)}
               ．e^{-w^2/(2u)}/√(2πu)
   = u^(-1/2)(v-w^2/u)^{(r-3)/2}e^(-v/2)/{2^(r/2)√πΓ((r-1)/2)}

因 S11 服從 χ^2(r)，故 S 的 j.p.d.f. 為

f(u,v,w) = (uv-w^2)^{(r-3)/2} e^{-(u+v)/2}/{2^r√πΓ(r/2)Γ((r-1)/2)}
   = |S|^{(r-3)/2} e^{-tr(S)/2}/{2^{(2r)/2}Γ_2(r/2)}

正是前述形式。

假設 p = 1, ..., k 時我們已得前述 Wishart 分布 p.d.f.，則當 p = k+1 時多出了 Sp1, ..., Spk, Spp 這些 S 的元素。我們把 S 的前 k×k 方陣記作 U, Spp 記作 V, Sp1,...,Spk 聯合記作 W。其中 Spt = Σ_i Xpi Xti = Stp，則在給定 U 之下，W 與 V*=V-W'U^(-1)W 相互獨立，W 服從 N(0,U)，V* 服從 χ^(r-k)。仿前述 p=2 方法，可得 p=k+1 時 W(I,r,p) 如前述形式。

假設 Xi 服從一般 N(0,V), 即使 V 並非正確定，只是非負確定 (non-negative definite), S = Σ Xi Xi' 之分布的 p.d.f. 可能不存在，但仍可宣稱是一個 Wishart 分布。若 V 的秩小於其階數，存在一列滿秩矩陣 A 使 AVA' 為正確定。反之，一個有 p.d.f. 的 Wishart 分布，經一個非列滿秩的矩陣做前述變換，結果成為一個不存在 p.d.f. 的 Wishart 分布。又，若 A 是列向量（單列矩陣），則 ASA' = Σ (AXi)^2 而 AXi 是 N(0,AVA')  隨機變數，故 S/(AVA') 服從 χ^2(r)，特例是 S 的每個主對角線元素都是尺度調整過的卡方 (scaled chi-squared) 變量。

欲證前述一般 W(V,r,p) 的結果，我們也可從特性函數 (ch.f.) 或動差母函數 (m.g.f.) 來看。為了避免複數的問題，我們考慮 m.g.f.

m(T) = E[e^{tr(TS)}]

其中 T 可以取上三角或下三角矩陣，因 S 是對稱矩陣，只需考慮其上三角或下三角（含主對角線）；但這相當於以 (T+T')/2 代替 T，因此也可要求 T 是對稱矩陣，只是其非主對角線元素用 t(ij)/2 而非 t(ij)。

隨機矩陣 X 的聯合 p.d.f. 為

p(x) = (2π)^(-rp/2) |V|^(-r/2) e^{- Σ Xi'V^(-1)Xi/2}

故

E[e^{tr(TS)}] = ∫...∫ C(V) e^{- Σ Xi'V^(-1)Xi/2 + tr(TS)} dx
  = ∫...∫ C(V) e^{-tr(V^(-1)XX'/2)+tr(TXX')} dx
  = ∫...∫ C(V) e^{-tr((V^(-1)-2T)XX')/2} dx
  = ∫...∫ C(V) e^{-tr((V^(-1)(I-2TV)XX')/2} dx
  = C(V)/C(V(I-2TV)^(-1))

故得 m(T) = |I-2TV|^(-r/2)。

令 X = [X1,...,Xr], 諸 Xi 是 N(0,V) 隨機行向量，行間相互機率獨立，則 Wishart 矩陣可表示為 S = XX'。考慮非負確定矩陣 A1,...,Ak, 設 I = A1+...+Ak, Ai^2=Ai,  AiAj=0 當 i≠j, 則「平方和式」被分解為一些「二次式」: Ｓ = XA1X' + ... + XAkX', 每個二次式如同 p=1 的情形，可以預期是 Wishart 分布，並且相互獨立。二次式構成矩陣 A 可以正交對角線化 A = PDP', D 是 A  的特徵根構成的對角線矩陣，由於 A^2 = A，故其特徹根非 0 即 1, 則

XAX' = XPDP'X' = (XP)D(XP)'

若 P 只取對應 D 中非 0 元素各行構成 P*，則 PDP' = P*P*’。記得 X 各行間是 i.i.d., 各列則是 N(0,V) 的成分，意即各列是元素相互 i.i.d. 的隨機變數，故 XP* 是將諸 Xi 重新組合成 r_j 個 i.i.d. N(0,V) 隨機向量。以 X*j 代表 XP* 第 j 行，X*j = Σi p(ij)Xi，則

Cov(X*j) = Cov(Σi p(ij) Xi) = ΣiΣt p(ij)p(tj) Cov(XiXt') = Σi p(ij)^2 V = V

式中 p(ij) 是 P* 或 P 的第 j 行第 i 列元素，由於 P* 或 P 各行規格正交 (orthonormal), 故 Σi p(ij)^2 = 1。類似地，我們也可以得知諸 X*j 相互獨立。

前述 S 被分解成 k 個 XAiX' 相加，諸 Ai 是對稱、非負確定、自乘不變，且 AiAj = 0 對所有 i≠j 成立，由線性代數的結果告訴我們：A1,...,Ak 可同時對角線化。因此，存在正交矩陣 P 使

I = PP' = PIP' = Σ_i PAiP' = Σ_i Di

其中諸 Di 非 0 對角線位置正好分割了單位矩陣 I 的 r 個位置，如果 Ai 的秩是 r_i, 則 r = Σ r_i；而諸二次式 XAiX' 構成 S = XX' 的正交分解，諸成分二次式之間相互獨立。所以，在單一多變量常態群體隨機抽樣：Xi～N(μ,V), i=1,...,n, 離均差二次式分解

Σ(Xi-μ)(Xi-μ)' = Σ(Xi-Xbar)(Xi-Xbar)' + n(Xbar-μ)(Xbar-μ)'

將總變異分解成樣本內變異加樣本平均向量變異，把 W(V,p,n) 隨機矩陣分解為相互獨立的 W(V,p,n-1) 和 W(V,p,1) 兩隨機矩陣；而在 MANOVA, 多變量迴歸或一般多變量線性模型，上述變異量可做更多分解。這些構成多變量常態群體為基礎的多變量方法基石。

Hotelling T-squared, T^2, 類似 t，以相互獨立的常態隨機向量與 Wishart 隨機矩陣構建：

T^2 = rZ'W^(-1)Z, 其中 Z 與 W 獨立，Z～N(0,Ip), W～W(Ip,p,r)

式中 Ip 是 p×p 單位矩陣。此分布以 T^2(p,r) 表示，r 為其自由度。若 Z～N(0,V), W～W(V,p,r), 求 V 之對稱平方根矩陣得 V^(1/2), 則

Z* = V^(-1/2)Z ～ N(0,Ip) ;  W* = V^(-1/2)WV^(-1/2) ～ W(Ip,p,r)

而 Z*'W*^(-1)Z* = Z'W^(-1)Z, 也就是說先前 T^2 的定義取 V = Ip 並不失其一般性。

在 p = 1 時，T^2 = rZ^2/W 是 F(1,r)；對一般的 p, Hotelling T^2 也可以轉換成 F 變量:

(r-p+1)T^2/(rp) = [(r-p+1)/p] Z'W^(-1)Z ～ F(p,r-p+1)

在單群體多變量常態平均數檢定，T^2 = n(Xbar-μ0)'S^(-1)(Xbar-μ0), 其中 S 是樣本共變異矩陣

S = Σ (Xi-Xbar)(Xi-Xbar)'/(n-1) ～ W(Σ,p,n-1)/(n-1)

其中 Σ = Cov(Xi), n 為樣本大小，μ0 為虛無参數向量，即 H0: μ = μ0。在 H0 成立時 √n(Xbar-μ0) ～ N(0,Σ)。T^2 統計量欲轉換成 F 統計量其公式為

F = (n-p)T^2/((n-1)p)

在兩獨立樣本（大小分別是 m, n）多變量平均數比較，假設 X 群體為 N(μ1,Σ), Y 群體為 N(μ2,Σ), 虛無假說 H0: μ1=μ2, 則

T^2 = (Xbar-Ybar)'[S(1/m+1/n)]^(-1)(Xbar-Ybar)

 其中 S = (Σ(Xi-Xbar)(Xi-Xbar)'+Σ(Yj-Ybar)(Yj-Ybar)')/(m+n-2), 而

F = (m+n-p-1)T^2/((m+n-2)p)

此處用了一個結果：一個 W(V,p,r1) 和另一個相互獨立的 W(V,p,r2) 隨機矩陣的和是 W(V,p,r1+r2) 變量（隨機矩陣）。

對於 Wishart 隨機矩陣，其特徵值 (eigenvalue) 的分布也是被關心的話題，S 的所有特徵值 j.p.d.f. 是

C(r,p) e^(-Σλi/2) Πλi^{(r-p-1)/2} Π_{i<j}|λi-λj|

其中 C(r,p) 是一個常數，λi, i=1,...,p, 是 S 的 p 個特徵值。