はじめに
現在この記事は加筆中です。ここでは多変量解析で使われる多次元正規分布及びその周辺の定理について纏めていこうと思います。
多次元正規分布
p次元の確率変数(確率ベクトル)Xが平均\(\mu_{p \times 1}\), 分散共分散行列 \(\Sigma_{p \times p}\) を持つことを以下の記号で表すとします。
$$X \sim \mathcal{N}_p(\mu_{p \times 1},\Sigma_{p \times p})$$
そもそもp次元の確率変数Xが多次元正規分布に従うことの定義は、任意のp次元のベクトル\(a \in \mathbb{R}^p \)に対して\( a^T X \)が１次元の正規分布に従う(または退化する)ことです。よく考えれば、この事は確率変数Xの特性関数\(\phi_X(a)=E[\exp(ia^TX)]\)が\(\exp(a^T\mu – \frac{1}{2}a^T \Sigma a)\)となることと同じだと言えます。
また\(X \sim \mathcal{N}_p(\mu_{p \times 1},\Sigma_{p \times p})\) であるときXの確率密度関数\(f(x) \)は次の形で表されます。
$$ \frac{1}{(2\pi)^{p/2}} \frac{1}{|\Sigma|^{1/2}} \exp(-\frac{1}{2} (x-\mu)^T \Sigma^{-1} (x-\mu)) $$
多次元正規分布の条件付き分布
p次元正規分布に従う確率変数Xを考えます。またXを\((X_1, X_2)\)の2つに分解します。(q次元とp–q次元に分解)
$$ X = \left( \begin{array}{c} X_1 \\ X_2 \end{array} \right) \sim \mathcal{N}_p(\mu, \Sigma)$$
さらに\( \mu, \Sigma \)を次のように分解しておきます。
$$\mu = \left( \begin{array}{c} \mu_1 \\ \mu_2 \end{array} \right), \Sigma =
\left( \begin{array}{cc} \Sigma_{11} & \Sigma_{12} \\ \Sigma_{21} & \Sigma_{22} \\ \end{array} \right) $$
ここでは \(X_2 | X_1 \) の条件付き分布の導出を書きます。ポイントは以下のような線形変換を考えることです。
$$ \left(\begin{array}{cc} I_q & O \\ -\Sigma_{21}\Sigma_{11}^{-1} & I_{p-q} \\ \end{array} \right) \left( \begin{array}{c} X_1 \\ X_2 \end{array} \right) $$
以下の関係式から簡単にこの分布を求められます。(以下の関係式は特性関数を用いればすぐに導出可能です)
$$X \sim \mathcal{N}_p(\mu, \Sigma) \rightarrow AX+b \sim \mathcal{N}_p(A\mu+b, A\Sigma A^T)$$
よって線形変換後の分布は次の式になります。
$$ \left(\begin{array}{cc} I_q & O \\ -\Sigma_{21}\Sigma_{11}^{-1} & I_{p-q} \\ \end{array} \right) \left( \begin{array}{c} X_1 \\ X_2 \end{array} \right) \sim \mathcal{N}_p(\left( \begin{array}{c} \mu_1 \\ \mu_2 – \Sigma_{21}\Sigma_{11}^{-1} \mu_1 \end{array} \right), \left(\begin{array}{cc} \Sigma_{11} & O \\ O & \Sigma_{22}-\Sigma_{21}\Sigma_{11}^{-1}\Sigma_{12} \\ \end{array} \right))$$
ここで嬉しいことに共分散が\( O \)になりました。多次元正規分布では共分散が\( O \)である(or無相関)ことと、独立であることは同値になります。この事も特性関数から簡単に証明できます。2つの確率変数\(W,Z\)について、\(W,Z\)が独立であれば\(W\)の周辺分布と、\(W|Z\)の条件付き分布は等しくなります。この事を用いれば先程の分布について、\( X_2-\Sigma_{21}\Sigma_{11}^{-1}X_1 \)の周辺分布と、\( X_2-\Sigma_{21}\Sigma_{11}^{-1}X_1 | X_1 \)の条件付き分布は等しいことがわかります。つまり、
$$ X_2-\Sigma_{21}\Sigma_{11}^{-1}X_1 | X_1=x_1 \sim \mathcal{N}_{p-q}(\mu_2-\Sigma_{21}\Sigma_{11}^{-1}\mu_1, \Sigma_{22} – \Sigma_{21} \Sigma_{11}^{-1} \Sigma_{12}) $$
となります。ここで\(X_1=x_1\)として条件を与えた瞬間に、\(x_1\)は定数と見なして構いません。そのため両辺に\(\Sigma_{21} \Sigma_{11}^{-1} x_1\)を加えることで、\(X_2 | X_1\)の分布は、
$$ X_2 | X_1=x_1 \sim \mathcal{N}_{p-q}(\mu_2+\Sigma_{21}\Sigma_{11}^{-1}(x_1-\mu_1), \Sigma_{22} – \Sigma_{21} \Sigma_{11}^{-1} \Sigma_{12}) $$
となります。
ウィシャート分布の定義
ウィシャート分布は英語ではWishart Distributionと書きます。今、\(X_1,X_2, …, X_n \sim^{iid} \mathcal{N}_{p}(0,\Sigma) \)である時、\(\sum_{j=1}^n X_j X_j^T \sim Wishart(\Sigma,n)\)であると定義します。見てもわかるようにウィシャート分布は\(p \times p\)の対称行列となります。\(p^2\)個の確率変数を一つにまとめた確率行列と考えれば良いでしょう。
また \(X= \left( \begin{array}{c} X_1^T \\ X_2^T \\ … \\ X_n^T \end{array} \right) \)のような行列のことをデータ行列(data matrix)と呼びます。データ行列を用いると、\(\sum_{j=1}^n X_j X_j^T = X^TX \)と書くことが出来ます。
ウィシャート分布に関する定理1
いま\(X_1,X_2, …, X_n \sim^{iid} \mathcal{N}_{p}(0,\Sigma) \)とし、\(X= \left( \begin{array}{c} X_1^T \\ X_2^T \\ … \\ X_n^T \end{array} \right) \)とする。さらに\(C_{n \times n}\)を対称な冪等行列(idempotent)で(rankC = r\)とする。ただし\(C\)は定数の行列または\(X_j (j=1,2 .. n)\)とは独立であるとする。このとき以下が成立する。
$$ X^T C X \sim Wishart(r,\Sigma)$$
(証明) まずは線形代数の基本的な事項を整理しておきます。(実数)対称行列は直交行列\(U_{n \times n}\)によって対角化が可能である。\(U\)が直交行列であるとは\(U^TU=I_n\)が成り立つことを言います。次に\(C\)が冪等行列とは\(C^2=C\)のような関係が成り立つことを言います。さらに冪等行列の固有値は1と0しか持ちません。(証明は容易) さらにこの事を使えば\(rank C = tr C = r\)と分かります。
いま以上の事から、ある直交行列\(U\)が存在して、
$$ C = U \Lambda U^T, \Lambda = \left(\begin{array}{cc} I_r & O \\ O & O_{n-r} \\ \end{array} \right) $$
と書くことが出来ます。
よって
$$X^TCX = X^T U \Lambda U^T X$$
が成立する。
ここで\(Y\)を次のように定義します。
$$Y= \left( \begin{array}{c} Y_1^T \\ Y_2^T \\ … \\ Y_n^T \end{array} \right) = U^T X $$
すると
$$X^TCX = Y^T \Lambda Y = \Sigma_{j=1}^r Y_jY_j^T$$
となること分かります。
そこでもし\(Y_1, Y_2, …, Y_n \sim^{iid} \mathcal{N}_p(0,\Sigma)\)であれば上記の命題が成立することが分かりますので、この事を証明します。いま、
$$Y^T = [Y_1 Y_2 … Y_n] = [X_1 X_2 … X_n] U$$
でありさらにUを成分毎に表示すると
$$[Y_1 Y_2 … Y_n] = [X_1 X_2 … X_n] \left(\begin{array}{cccc} u_{11} & u_{12} & … & u_{1n} \\ u_{21} & u_{22} & … & u_{2n} \\ … & … & … & … \\ u_{n1} & u_{n2} & … & u_{nn} \\ \end{array} \right)$$
ゆえに
$$Y_i = \Sigma_{k=1}^n u_{ki} X_k, Y_j = \Sigma_{l=1}^n u_{lj} X_l,$$
となります。
そこで\(Cov[Y_i,Y_j]\)を考えると、
$$Cov[Y_i,Y_j] = \Sigma_{k=1}^n \Sigma_{l=1}^n u_{ki} u_{lj} Cov[X_k,X_l]$$
今\(X_1, X_2, …, X_n\)は独立なので\(k=l\)以外の共分散は\(O\)になるはずです。そこで
$$ = \Sigma_{k=1}^n u_{ki} u_{kj} Cov[X_k,X_k] = \Sigma \cdot (\Sigma_{k=1}^n u_{ki} u_{kj})$$
となります。
また\(U = [u_1 u_2 … u_n]\)というように、今度は\(U\)を列ベクトル\( u_1, …, u_n \)で表現してやると、
$$ Cov[Y_i,Y_j]= \Sigma \cdot (\Sigma_{k=1}^n u_{ki} u_{kj}) = \Sigma \cdot (u_i^T u_j) = \Sigma ({\rm if}: i=j) , O ({\rm if}: i \neq j)=\Sigma \cdot \delta_{i,j}$$
となります。なぜならば\(u_1, … u_n\)は正規直交ベクトルだからです。
次に\(E[Y_i] = 0\)であることは容易に分かるでしょう。また\(X^{*} = \left( \begin{array}{c} X_1 \\ X_2 \\ … \\ X_n \end{array} \right) \)という\(np\)次元の確率変数は\(np\)次元正規分布に従うことも容易に分かるでしょう。
つまり
$$ X^{*} \sim \mathcal{N}_{np}(\left( \begin{array}{c} 0 \\ 0 \\ … \\ 0 \end{array} \right), \left( \begin{array}{cccc} \Sigma & O & … & O \\ O & \Sigma & … & O \\ … & … & … & … \\ O & O & … & \Sigma \end{array} \right)$$
いま\(Y^{*} = \left( \begin{array}{c} Y_1 \\ Y_2 \\ … \\ Y_n \end{array} \right) \)は\(X^{*}\)の線形変換によって得られるので、\(Y^{*}\)もまた\(np\)次元の正規分布に従うことが分かります。さらに先程の結果\( Cov[Y_i,Y_j] = \Sigma \cdot \delta_{i,j} , E[Y_i]=0 \)をあわせると、\(Y^{*}\)も同様に、
$$Y^{*} \sim \mathcal{N}_{np}(\left( \begin{array}{c} 0 \\ 0 \\ … \\ 0 \end{array} \right), \left( \begin{array}{cccc} \Sigma & O & … & O \\ O & \Sigma & … & O \\ … & … & … & … \\ O & O & … & \Sigma \end{array} \right) $$
であることが分かります。また多次元正規分布では、無相関と独立は同値であるので\(Y_1, Y_2, …, Y_n\)が独立であり、さらにそれぞれが\(N_{p}(0,\Sigma)\)に従うことが分かりますので、題意が示されました。
ウィシャート分布に関する定理2
\(M \sim Wishart(n,\Sigma_{p \times p})\)であるとする。
さらに\( \left( \begin{array}{cc} M_{11} & M_{12} \\ M_{21} & M_{22} \end{array} \right) \)と分解をしておきます。ただし\(M_{11}\)は\(a \times a\)の正方行列であり、\(M_{22}\)は\((p-a)\times(p-a)\)の正方行列とします。
このとき以下の性質が成立します。
$$ M_{22.1} := M_{22}-M_{21}M_{11}^{-1}M_{12} \sim Wishart(\Sigma_{22.1},n-a) , {\rm where} \Sigma_{22.1}=\Sigma_{22}-\Sigma_{21}\Sigma_{11}^{-1}\Sigma_{12} $$
(証明) まず\(M \sim Wishart(\Sigma,n)\)であることから\(M = X^TX\)と書き表すことができ、しかも\(X= \left( \begin{array}{c} X_1^T \\ X_2^T \\ … \\ X_n^T \end{array} \right), X_j \sim^{iid} \mathcal{N}_p(0,\Sigma) \)とすることができます。
次に各\(X_j\)に関して前から\(a\)個と後ろ\(p-a\)個に分けておきます。つまり\(X_j = \left( \begin{array}{c} X_{j,1} \\X_{j,2} \end{array} \right) \)とします。
また\(X\)を縦に\(X=[X_{(1)} X_{(2)}]\)というように区切っておきます。つまり\(X_{(1)}=\left( \begin{array}{c} X_{1,1}^T \\X_{2,1}^T \\ … \\ X_{n,1}^T \end{array} \right)\)かつ\(X_{(2)}=\left( \begin{array}{c} X_{1,2}^T \\X_{2,2}^T \\ … \\ X_{n,2}^T \end{array} \right)\)です。
すると\(M=\left( \begin{array}{c} X_{(1)}^TX_{(1)} & X_{(1)}^TX_{(2)} \\ X_{(2)}^TX_{(1)} & X_{(2)}^TX_{(2)} \end{array} \right)\)となります。
そのことから、\(M_{22.1}=X_{(2)}^TX_{(2)}-X_{(2)}^TX_{(1)} (X_{(1)}^TX_{(1)})^{-1}X_{(1)}^T X_{(2)} = X_{(2)}^T(I_{n}-X_{(1)} (X_{(1)}^TX_{(1)})^{-1}X_{(1)}^T)X_{(2)}\)となります。
また簡単のため\(H_{(1)} = X_{(1)} (X_{(1)}^TX_{(1)})^{-1}X_{(1)}^T\)と置きます。よって\(M_{22.1}=X_{(2)}^T(I_n-H_{(1)})X_{(2)}\)となります。
ここで\(H_{(1)}\)はハット行列(Hat Matrix)と言われる行列で、\(H_{(1)}^2=H_{(1)}, H_{(1)}^T=H_{(1)}\)となり、さらに\(H_{(1)}X_{(1)}=X_{(1)}\)となることに注意しましょう。(証明は容易) つまり\(H_{(1)}\)は対称な冪等行列であることに注意しておきます。さらに\(I_n-H_{(1)}\)も同様に対称な冪等行列になっています。(証明は容易)
いま\(H_{(1)}X_{(1)}=X_{(1)}\)に注意すれば、\(M_{22.1} = (X_{(2)}-X_{(1)}\Sigma_{11}^{-1}\Sigma_{12})^T (I-H_{(1)})(X_{(2)}-X_{(1)}\Sigma_{11}^{-1}\Sigma_{12}) \)となることが分かります。
(※さらに言えば、\((I_n-H_{(1)})X_{(1)}=O_{n}\)になるので結局、加えた\(X_{(1)}\Sigma_{11}^{-1}\Sigma_{12}\)は打ち消されることとなる。)
ここで\(X_{(2.1)} := X_{(2)}-X_{(1)}\Sigma_{11}^{-1}\Sigma_{12}\)と置くと、\(M_{22.1} = X_{(2.1)}^T(I-H_{(1)})X_{(2.1)}\)となります。
次に少し計算をすると \(X_{(2.1)} = \left( \begin{array}{c} … \\ (X_{j,2} – \Sigma_{21}\Sigma_{11}^{-1}X_{j,1})^T \\ … \end{array} \right)_{n \times (p-a)}\)の形で書けることが分かります。
また条件付き分布の定理で出てきたように、\( X_{j,2} – \Sigma_{21}\Sigma_{11}^{-1}X_{j,1} \sim^{iid} \mathcal{N}_{p}(0,\Sigma_{22}-\Sigma_{21}\Sigma_{11}^{-1}\Sigma_{12}) = \mathcal{N}_{p}(0,\Sigma_{22.1}) \)が成り立ちます。
加えて、各\(j\)について\(X_{j,2} – \Sigma_{21}\Sigma_{11}^{-1}X_{j,1}\)と\(X_{j,1}\)は独立であったことに注意します。すると、\(X_{(2.1)}\)と\(X_{(1)}\)は独立であることが分かります。よって\(X_{(2.1)}\)と\(H_{(1)}\)も独立となります。(\(X_{(2.1)}\)と\(I-H_{(1)}\)も独立)
よって先程の「ウィシャート分布に関する定理1」で\(C=I_n-H_{(1)}\)と見立てて、\( rank(I_n-H_{(1)})=tr(I_n-H_{(1)})=n-a\)となります。
(※\( trAB=trBA, tr(A+B)=trA+trB \)の関係を使います。→ \( trH_{(1)} = tr X_{(1)}(X_{(1)}^T X_{(1)})^{-1} X_{(1)}^T \)であり、行列の積が定義できるならば順序を入れ替えてもtraceは変わらないので、\(= tr (X_{(1)}^T X_{(1)})^{-1} X_{(1)}^T X_{(1)} = tr I_{a} = a \) であることに注意すれば分かるでしょう。)
よって以上のことから証明が完了しました。
ホテリングのTスクエアの定義
英語では\( {\rm Hotelling’s} \hspace{4pt} T^2\)などと書かれるようです。
\(d \sim \mathcal{N}_{p}(0,I_{p}) , M \sim Wishart(n,I_{p})\)とし\(d,M\)は独立とする。
このとき\(nd^TM^{-1}d\)の分布をホテリングのTスクエア分布と呼び、\(nd^TM^{-1}d \sim T^2(p;n)\)などと表記する。
次に証明するがHotelling’s T-squared 分布はF分布の定数倍によって表される。
ホテリングのTスクエアとF分布の関係についての証明
ホテリングのTスクエアとF分布の間には、\( T^2(p;n) = \frac{np}{n-p+1} F_{p,n-p+1}\)という関係が成り立つ。
(証明 STEP1)
まず\(d \in \mathbb{R}^p\)を0でないベクトルとする。次に \(D = [D_1 d]_{p \times p}\)という正則行列を考える。さらに\(d^TD_1 = 0\)を満たすように\(D_1\)を取っているものとする。
次に\(M = X^TX \sim Wishart(I_{p},n)\)とする。ここで\(X\)はデータ行列とし、各行を構成する\(X_j^T\)について\(X_j \sim^{iid} \mathcal{N}_p(0,I_p)\)であることに注意しておく。今\( D^{-1}X_j \sim^{iid} \mathcal{N}_p(0,(D^TD)^{-1}) \)となるので、\( D^{-1}M{D^{-1}}^T = D^{-1}X^TX{D^{-1}}^T \sim Wishart((D^TD)^{-1},n)\)となる。
いま\(N := D^{-1}M{D^{-1}}^T \)と書くこととします。さらに\(N = \left( \begin{array}{cc} N_{11} & N_{12} \\ N_{21} & N_{22} \end{array} \right)\)と分解しておきます。ただし\(N_{11}\)の大きさは\((p-1)\times(p-1)\)で\(N_{22}\)は\(1 \times 1\)の大きさ(つまりスカラー)とします。また\(N\)と同じレイアウトで\((D^TD)^{-1}\)も分解しておきます。ただしはじめに与えた条件に注意しておくと\((D^TD)^{-1} = \left( \begin{array}{cc} (D_{1}^TD_{1})^{-1} & O \\ O & (d^Td)^{-1} \end{array} \right) \)となります。
次に「ウィシャート分布に関する定理2」を用いて\(N_{22.1}\)の分布を求めると、\(N_{22.1} \sim Wishart((D^TD)^{-1}_{22.1}, n-(p-1))\) となります。さらに書き換えると\(Wishart((d^Td)^{-1}, n-p+1))\)となります。
またウィシャート分布において行列パラメータ\(\Sigma\)がスカラー(例えば\(\sigma^2\))になるときはカイ二乗分布の定数倍となります。(証明は容易) そして \(d^Td N_{22.1} \sim Wishart(1,n-p+1) = \chi^2_{(n-p+1)}\)となります。
最後にブロック行列の逆行列の公式から、\((N^{-1})_{22} = N_{22.1}^{-1} \)の関係が成り立ちます。いま\((N^{-1})_{22} = (D^TM^{-1}D)_{22} = d^TM^{-1}d\)だったことに注意すると\(N_{22.1} = \frac{1}{d^TM^{-1}d}\)となります。よって\(\frac{d^Td}{d^TM^{-1}d} \sim \chi^2_{n-p+1}\)ということがわかりました。
(証明 STEP2)
次に\(d\)を\(d \sim \mathcal{N}_p)(0,I_p\)であると仮定します。
さらに\(d\)とは独立なウィシャート分布に従う確率行列\(M = X^TX \sim Wishart(I_{p},n)\)を考えます。
いま連続型の確率変数であることから\(d \neq 0_p (a.s)\)が成り立つことに注意します。またSTEP1 の証明では、\(0_p\)ではないどのようなベクトル\(d\)を与えても\(d^T M^{-1} d \sim \chi^2_{(n-p+1)}\)であったことに注意します。よって\(d \neq 0_p (a.s)\)であることから、\(d \sim \mathcal{N}_p(0,I_p)\)の条件下でも、\(\frac{d^Td}{d^T M^{-1} d} \sim \chi^2_{n-p+1}\)であることに注意します。
また\(\frac{d^Td}{d^T M^{-1} d} = 1/(\frac{d}{|d|} M^{-1} \frac{d}{|d|}) \)と書き換えると、\(\frac{d^Td}{d^T M^{-1} d}\)は\(|d|=(d^Td)^{1/2}\)とは無関係に決まるとわかります。(\(\frac{d}{|d|}\)は単位化されたベクトルなので既に大きさに関する情報は保持していないからです。) つまり\(\frac{d^Td}{d^T M^{-1} d}\)と\(d^Td\)は独立となります。いま\(d^Td \sim \chi^2_{(p)}\)であることに注意しつつ、F分布の定義に立ち返ると、\(\frac{d^Td / p}{\frac{d^Td}{d^T M^{-1} d} / (n-p+1)} \sim F_{p,n-p+1}\)となります。
あとは式を整理することで命題の主張が得られます。
ウィルクスのラムダの定義
本記事の最後にウィルクスのラムダについて述べたいと思います。
\(\)