尤度比検定の検定統計量がχ二乗分布に従う理由を確認します。
1 ブロック行列の逆行列
尤度比検定の検定統計量を展開していくときに使うブロック行列の逆行列の求め方を確認しましょう。ほとんど四則演算です。
\(k\)行\(k\)列の正方行列を、\(r\)行と\(k-r\)行、\(r\)列と\(k-r\)列に分割します。ゼロ行列を\(O\)，恒等行列を\(1_n\)（\(n\in\{k, k-r\}\)）とおき、正方行列\(H\)を三角行列と対角行列に分解することを考えます。
\[\begin{aligned} H &= \begin{pmatrix} H_1 & H_2 \\ H_3 & H_4 \end{pmatrix} \\ &= \begin{pmatrix} 1_k & M_1 \\ O & 1_{k-r} \end{pmatrix} \begin{pmatrix} M_2 & O \\ O & M_3 \end{pmatrix} \begin{pmatrix} 1_k & O \\ M_4 & 1_{k-r} \end{pmatrix} \\ &= \begin{pmatrix} M_2 & M_2 M_4 \\ M_1 M_2 & M_1 M_2 M_4 + M_3 \end{pmatrix} \end{aligned}\]
となり、\(M_1 = H_1^{-1}H_2\), \(M_2 = H_1\), \(M_3 = H_4 - H_3 H_1^{-1} H_2\), \(M_4 = H_1^{-1}H_2\)と左右の行列の要素の対応が分かります。
逆行列を考えます。
\[\begin{aligned} H^{-1} &= \begin{pmatrix} H_1 & H_2 \\ H_3 & H_4 \end{pmatrix}^{-1} \\ &= \bigg( \begin{pmatrix} 1_k & M_1 \\ O & 1_{k-r} \end{pmatrix} \begin{pmatrix} M_2 & O \\ O & M_3 \end{pmatrix} \begin{pmatrix} 1_k & O \\ M_4 & 1_{k-r} \end{pmatrix} \bigg)^{-1} \\ &= \begin{pmatrix} 1_k & O \\ M_4 & 1_{k-r} \end{pmatrix}^{-1} \begin{pmatrix} M_2 & O \\ O & M_3 \end{pmatrix}^{-1} \begin{pmatrix} 1_k & M_1 \\ O & 1_{k-r} \end{pmatrix}^{-1} \\ &= \begin{pmatrix} 1_k & O \\ -M_4 & 1_{k-r} \end{pmatrix} \begin{pmatrix} M_2^{-1} & O \\ O & M_3^{-1} \end{pmatrix} \begin{pmatrix} 1_k & -M_1 \\ O & 1_{k-r} \end{pmatrix} \\ &= \begin{pmatrix} (H_1 - H_2 H_4^{-1} H_3)^{-1} & \cdots \\ \vdots & \ddots \end{pmatrix} \end{aligned}\]
左上のブロックだけ計算しました。他もがんばって展開して代入すれば出てくるわけですが、後の説明に使わないので省略します。
2 二次形式\({}^t\!(\xi - E[\xi]) \Sigma^{-1} (\xi - E[\xi])\)は\(\chi^2\)分布に従う
\(\xi\)は正規分布に従う要素数\(r\)のベクトル、\(\Sigma\)を\(\xi\)の分散共分散行列とします。
ここで\(z=\xi - E[\xi]\)として
\[\begin{align} {}^t\!(\xi - E[\xi]) \Sigma^{-1} (\xi - E[\xi]) &= {}^t\!z \Sigma^{-1} z \\ &= {}^t\!(\Sigma^{-\frac{1}{2}}z)(\Sigma^{-\frac{1}{2}}z) \end{align}\]
\(w=\Sigma^{-\frac{1}{2}}z\)として
\[\begin{align} {}^t\!(\xi - E[\xi]) \Sigma^{-1} (\xi - E[\xi]) &= {}^t\!w w \end{align}\]
\(E[z]=0\)なので\(E[w]=0\)になります。\(z\)の分散共分散行列も\(\Sigma\)になることに注意して、
\[\begin{align} \mbox{VAR}[w] &= (\Sigma^{-\frac{1}{2}}z) \ {}^t\!(\Sigma^{-\frac{1}{2}}z) \\ &= \Sigma^{-\frac{1}{2}} z \ {}^t\! z \ {}^t\Sigma^{-\frac{1}{2}} \\ &= \Sigma^{-\frac{1}{2}} \Sigma \ {}^t\Sigma^{-\frac{1}{2}} \\ &= 1_r \end{align}\]
ベクトル\(w\)は期待値\(0\)で分散\(1\)で共分散\(0\)になります。\({}^t\!(\xi - E[\xi]) \Sigma^{-1} (\xi - E[\xi])\)は\(w\)の\(2\)乗なので、自由度が\(r\)になる\(\chi^2\)分布に従います。
3 フィッシャー情報量の変形
\(f(x, \theta)\)を確率密度関数とします。\(x\)は観測値、\(\theta\)はパラメーターです。
\[\begin{aligned} \frac{\partial}{\partial \theta} \int f(x, \theta) dx &= \int \frac{\partial}{\partial \theta} f(x, \theta) dx \\ &= \int \bigg( \frac{\partial}{\partial \theta} \log f(x, \theta) \bigg) f(x, \theta) dx \\ &= E\bigg[ \frac{\partial}{\partial \theta} \log f(x, \theta) \bigg] \end{aligned}\]
確率密度関数の定義から常に\(\int f(x, \theta) dx = 1\)なので、\(\frac{\partial}{\partial \theta} \int f(x, \theta) dx\)は常に\(0\)になります。よって、
\[ E\bigg[ \frac{\partial}{\partial \theta} \log f(x, \theta) \bigg] = 0 \]
と常になります。上式の両辺をさらに\(\theta\)で微分すると、
\[\begin{aligned} \frac{\partial}{\partial \theta} \int \frac{\partial}{\partial \theta} \log f(x, \theta) dx &= \int \frac{\partial^2}{\partial \theta^2} \log f(x, \theta) f(x) dx \\ &\ + \int \bigg( \frac{\partial}{\partial \theta} \log f(x, \theta) \bigg)^2 f(x) dx \\ &= E\bigg[ \frac{\partial^2}{\partial \theta^2} \log f(x, \theta) \bigg] + I(\theta) \\ &= 0 \\ I(\theta) &= - E\bigg[ \frac{\partial^2}{\partial \theta^2} \log f(x, \theta) \bigg] \end{aligned}\]
となり、フィッシャー情報量を計量分析で扱いやすい形に変形できます。
4 Wilks定理
\(k\)個の要素のパラメーター・ベクトル\(\theta\)を考えます。要素\(1\)から\(r\)番目は\(0\)とした帰無仮説のパラメーター・ベクトルを\(\theta_0\)とします。要素の順番は任意に置換できるので、制約付きモデルはこの形にすることができます。
\[ \theta_0 = \begin{pmatrix} 0 \\ \theta_{0*} \end{pmatrix} \]
推定されるパラメーターを\(\hat{\theta}\)と置きます。\(\hat{\theta}\)のときの対数尤度関数\(L(X, \hat{\theta})\)のグラディエントは\(0\)です。
\[ \nabla \log L(X, \hat{\theta}) \equiv \begin{pmatrix} \nabla_A \\ \nabla_C \end{pmatrix} \log L(X, \hat{\theta}) = \begin{pmatrix} 0 \\ 0 \end{pmatrix} \]
\(X\)は観測値\(x\)の集合、データです。制約付きモデルで推定されるパラメーターを\(\hat{\theta}^r\)と置きます。対数尤度関数\(L(x, \hat{\theta}^r)\)のグラディエントは\(0\)です。
\[ \nabla_C \log L(X, \hat{\theta}^r) = 0 \]
フィッシャー情報行列\(I(\theta_0)\)は以下のように書けます。
\[\begin{aligned} I(\theta_0) &= -\frac{1}{n} E \big[ \nabla^{\otimes 2} \log L(X, \theta_0)\big] \\ &= \frac{1}{n} E \big[ \nabla \log L(X, \theta_0)^{\otimes 2}\big] \\ &= \begin{pmatrix} A & B \\ ^t\!B & C \end{pmatrix} \end{aligned}\]
\(n\)で割っているのは、尤度関数は\(n\)個の観測値による確率密度関数の積になるためです。
\(A\)は\(r \times r\)，\(B\)は\(r \times (k-r)\)、\(C\)は\((k-r) \times (k-r)\)行列です。
スコア統計量\(S\)は以下のように書けます。
\[ \frac{1}{\sqrt{n}} \equiv \frac{1}{\sqrt{n}} \nabla \log L(X, \theta_0) \equiv \begin{pmatrix} \xi \\ \eta \end{pmatrix} \]
\(\nabla \log L(X, \hat{\theta})\)を\(\theta_0\)のまわりで一階のテーラー展開をすると、
\[ \nabla \log L(X, \hat{\theta}) = \nabla \log L(X, \theta_0) + \nabla^2 \log L(X, \theta_0)(\hat{\theta} - \theta_0) + \cdots \]
\(\nabla \log L(X, \hat{\theta}) = 0\)なので、
\[ 0 = \nabla \log L(X, \theta_0) + \nabla^2 \log L(X, \theta_0)(\hat{\theta} - \theta_0) + \cdots \]
両辺を\(\sqrt{n}\)で割り、フィッシャー情報行列を使って書き換えると、
\[ 0 \approx \frac{1}{\sqrt{n}} \nabla \log L(X, \theta_0) - I(\theta_0) \sqrt{n} (\hat{\theta} - \theta_0) \]
よって
\[ \begin{pmatrix} \xi \\ \eta \end{pmatrix} \approx I(\theta_0) \sqrt{n} (\hat{\theta} - \theta_0) \]
\[ \eta \approx C \sqrt{n} (\hat{\theta}_{*} - \theta_{0*}) \]
尤度比検定の統計量は以下のように定義されています。
\[ -2 \log \Lambda \equiv \log \frac{T_1}{T_2}, \quad T_1 \equiv \frac{L(X, \hat{\theta})}{L(X, \theta_0)}, \quad T_2 \equiv \frac{L(X, \hat{\theta}^r)}{L(X, \theta_0)} \]
\(\log T_1\)をテイラー展開し、\(\nabla \log L(X, \hat{\theta}) = 0\)に注意して整理します。
\[\begin{aligned} \log T_1 &= \log L(X, \hat{\theta}) - \log L(X, \theta_0) \\ &= \log L(X, \theta_0) \\ &\ + \nabla \log L(X, \hat{\theta})(\hat{\theta} - \theta_0) \\ &\ + \frac{1}{2}\ \!^t\!(\hat{\theta} - \theta_0) \nabla^{\otimes 2} \log L(X, \theta_0) (\hat{\theta} - \theta_0) \\ &\ + \cdots - \log L(X, \theta_0) \\ &\approx \frac{1}{2}\ {}^t\!(\hat{\theta} - \theta_0) \nabla^{\otimes 2} \log L(X, \theta_0) (\hat{\theta} - \theta_0) \end{aligned}\]
\(\sqrt{n}\frac{1}{n}\sqrt{n}=1\)に注意すると、\(2 \log T_1\)は以下となります。
\[\begin{aligned} 2 \log T_1 &\approx \sqrt{n}\ ^t\!(\hat{\theta} - \theta_0) \frac{1}{n} \nabla^{\otimes 2} \log L(X, \theta_0) (\hat{\theta} - \theta_0) \sqrt{n} \\ &= \begin{matrix} {} \\ {} \end{matrix}^t\!\!\!\begin{pmatrix} \xi \\ \eta \end{pmatrix} I^{-1}(\theta_0) \begin{pmatrix} \xi \\ \eta \end{pmatrix} \end{aligned}\]
同様に、
\[ 2 \log T_2 \approx {}^t\!\eta \ C^{-1} \eta \]
尤度比検定の統計量は、以下のように近似できます。
\[ -2 \log \Lambda \approx \begin{matrix} {} \\ {} \end{matrix}^t\!\!\! \begin{pmatrix} \xi \\ \eta \end{pmatrix} I^{-1}(\theta_0) \begin{pmatrix} \xi \\ \eta \end{pmatrix} - {}^t\eta C^{-1}\eta \]
この式を展開するための補題を整理します。
\[\begin{align} \frac{1}{\sqrt{n}} \nabla_A \log L(X, \hat{\theta}^r) &\approx \nabla_A \frac{1}{\sqrt{n}} \log L(X, \theta_0) \\ &\ + \frac{1}{n} \nabla_A {}^t \nabla \log L(X, \theta_0) \begin{pmatrix} 0 \\ \sqrt{n} (\hat{\theta}_* - \theta_{0*}) \end{pmatrix} \end{align}\]
\(\log L(X, \hat{\theta}^r)\)に\(\theta_0\)まわりでテイラー展開を用いました。
\[ \begin{pmatrix} \nabla_A \\ O \end{pmatrix} \begin{pmatrix} {}^t\nabla_A & {}^t\nabla_C \end{pmatrix} = \begin{pmatrix} \nabla_A {}^t\nabla_A & \nabla_A {}^t\nabla_C \\ O {}^t\nabla_A & O {}^t\nabla_C \end{pmatrix} \]
と、\(n^{-1}\nabla^{\otimes 2}\log L(X, \theta_0) \approx - I(\theta_0)\)と、\(\eta \approx C \sqrt{n} (\hat{\theta}_{*} - \theta_{0*})\)に注意して、
\[ \frac{1}{\sqrt{n}} \nabla_A \log L(X, \hat{\theta}^r) \approx \xi - B\sqrt{n}(\hat{\theta}_* - \theta_{0*}) = \xi - BC^{-1}\eta \]
\({}^t\!({}^t\!\xi\ {}^t\!\eta)\)を分割します。
\[ \begin{pmatrix} \xi \\ \eta \end{pmatrix} = v_A + v_C \]
\[\begin{align} v_A &\equiv \begin{pmatrix} \xi - BC^{-1}\eta \\ 0 \end{pmatrix} \\ v_C &\equiv \begin{pmatrix} BC^{-1}\eta \\ \eta \end{pmatrix} \end{align}\]
以下に注意して、
\[\begin{align} \big( {}^t\!B \quad C \big) I^{-1}(\theta_0) &= \big( {}^t\!B \quad C \big) \begin{pmatrix} A & B \\ {}^t\!B & C \end{pmatrix}^{-1} = ( 0 \quad 1_{k-r} ) \end{align}\]
\(v_CI^{-1}(\theta)\)は次のように書き換えられます。
\[\begin{align} v_CI^{-1}(\theta) &= \begin{matrix} {} \\ {} \end{matrix}^t\!\!\!\begin{pmatrix} BC^{-1}\eta \\ \eta \end{pmatrix} I^{-1}(\theta_0) \\ &= {}^t\big( C^{-1} \eta \big) \big( {}^t\!B \quad C \big) I^{-1}(\theta_0) \\ &= {}^t\big( C^{-1}\eta\big) ( 0 \quad 1_{k-r} ) \end{align}\]
\(v_A\)は\(r\)番目より後は\(0\)なので、ブロック行列にして左上のブロック以外は\(0\)になることに注意して、ブロック行列の逆行列の公式をあてはめます。
\[\begin{align} {}^t\!v_A I^{-1}(\theta_0) v_A &= {}^t (\xi - BC^{-1}\eta) (A - BC^{-1}{}^t\!B)^{-1} (\xi - BC^{-1}\eta) \end{align}\]
\[\begin{align} {}^t\!v_C I^{-1}(\theta_0) v_A &= {}^t\big( C^{-1}\eta\big) ( 0 \quad 1_{k-r} ) \begin{pmatrix} \xi - BC^{-1}\eta \\ 0 \end{pmatrix} \\ &= 0 \end{align}\]
対称行列の逆行列も対称行列なことに注意します。
\[\begin{align} {}^t\!v_A I^{-1}(\theta_0) v_C &= {}^t ({}^t\!v_C \ {}^t\!I^{-1}(\theta_0) v_A) \\ &= {}^t ({}^t\!v_C I^{-1}(\theta_0) v_A) = 0 \\ \end{align}\]
\[\begin{align} \\ {}^t\!v_C I^{-1}(\theta_0) v_C &= {}^t\big( C^{-1}\eta\big) ( 0 \quad 1_{k-r} ) \begin{pmatrix} BC^{-1}\eta \\ \eta \end{pmatrix} \\ &= {}^t ({}^t \eta C^{-1} \eta) \\ &= {}^t \eta C^{-1} \eta \end{align}\]
\[ \begin{matrix} {} \\ {} \end{matrix}^t\!\!\! \begin{pmatrix} \xi \\ \eta \end{pmatrix} I^{-1}(\theta_0) \begin{pmatrix} \xi \\ \eta \end{pmatrix} = {}^t\!(v_A + v_C) I^{-1}(\theta_0) (v_A + v_C) \]
\[ = (\xi - BC^{-1}\eta) (A - BC^{-1}{}^t\!B)^{-1} (\xi - BC^{-1}\eta) \\ + {}^t \eta C^{-1} \eta \]
よって、
\[ -2 \log \Lambda \approx {}^t\!(\xi - BC^{-1}\eta) (A - BC^{-1}{}^t\!B)^{-1} (\xi - BC^{-1}\eta) \]
となることが分かります。
\({}^t\!({}^t\!\xi\ {}^t\!\eta)\)は\(\mathit{N}(0, I^{-1}(\theta_0))\)に漸近するため、\(\xi - BC^{-1}\eta\)も正規分布に漸近します。分散は
\[ E\big[ (\xi - BC^{-1}\eta) {}^t\!(\xi - BC^{-1}\eta) \big] \approx (A - BC^{-1}{}^t\!B)^{-1} \]
となるので1、\(-2 \log \Lambda\)が自由度\(r\)の\(\chi^2\)分布に従うと看做せることが分かります。
\(\xi\)がテーラー展開の1次の項、\(-BC^{-1}\eta\)が2次の項になるため、\(-BC^{-1}\eta\)の交差項は無視できます。\(\xi{}^t\!\xi\)を考えればよくなりますが、これは\(I(\theta)\)の左上のブロックになるため、その分散共分散行列は\(I^{-1}(\theta)\)の左上のブロックになります。↩︎