はじめに
- 統計検定一級の2014年統計数理の問4では、非心カイ二乗分布について出題されている。
- 公式の解答例を読んでも全く理解できず、検索しても全く手掛かりが得られない。
- 自分が納得できる解答を考えてみた。
非心カイ二乗分布の定義
$X\sim N(\mu,1)$の場合、$Y\sim X^2\ $を自由度1、非心度$\lambda=\mu^2$の非心カイ二乗分布と定義する。
確率密度関数と積率母関数は以下のようになる。
\begin{align*}
f(x)&=\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{(x-\mu)^2}{2}\right),\quad y=x^2,x=\sqrt{y}\\
f(y)&=\frac{1}{2\sqrt{y}}\frac{1}{\sqrt{2\pi}}\left[\exp\left(-\frac{(y-2\sqrt{y}\mu+\mu^2)}{2}\right)
+\exp\left(-\frac{(y+2\sqrt{y}\mu+\mu^2)}{2}\right)\right]\\
&=\frac{1}{2\sqrt{y}}\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{y+\mu^2}{2}\right)
\left[\exp(\sqrt{y}\mu) + \exp(-\sqrt{y}\mu)\right]\\
&=\frac{1}{\sqrt{2\pi y}}\exp\left(-\frac{y+\lambda}{2}\right)
\cosh\left(\sqrt{y\lambda}\right)\quad(\lambda=\mu^2)\\
M_Y(t)&=\int_0^\infty \frac{1}{\sqrt{2\pi y}}\exp\left(-\frac{y(1-2t)+\lambda}{2}\right)
\cosh(\sqrt{y\lambda}) dy\\
&=\int_0^\infty \frac{1}{\sqrt{2\pi z(1-2t)}}\exp\left(-\frac{z+\lambda}{2}\right)
\cosh\left(\sqrt{z\frac{\lambda}{1-2t}}\right) dz\\
&\lambda'=\frac{\lambda}{1-2t},\lambda'-\lambda=\frac{2\lambda t}{1-2t}\\
&=\int_0^\infty \frac{1}{\sqrt{2\pi z(1-2t)}}
\exp\left(-\frac{z+\lambda'}{2}\right)
\exp\left(-\frac{\lambda-\lambda'}{2}\right)
\cosh\left(\sqrt{z\lambda'}\right) dz\\
&=\frac{\exp\left(\frac{\lambda t}{1-2t}\right)}{\sqrt{1-2t}}\\
\end{align*}
積率母関数の形式から、パラメータ$\lambda$が加法的となることが分かる。
自由度k、非心度$\lambda$の非心カイ二乗分布の確率密度関数と積率母関数は以下の通り。
(https://en.wikipedia.org/wiki/Noncentral_chi-squared_distribution)
\begin{align*}
Z&=\sum_{i=1}^kX_i^2,\quad X_i\sim N(\mu_i,1),
\quad \lambda=\sum_{i=1}^k \mu_i^2\\
f(z|k,\lambda)&=\frac{1}{2}\exp\left(-\frac{z+\lambda}{2}\right)
\left(\frac{z}{2}\right)^{\frac{k}{2}-1}
\sum_{n=0}^\infty\frac{1}{n!\Gamma\left(\frac{k}{2}+n\right)}\left(\frac{\lambda z}{4}\right)^n\\
M_Z(t)&=\frac{\exp\left(\frac{\lambda t}{1-2t}\right)}{(1-2t)^{\frac{k}{2}}}
\end{align*}
問題の抜粋
ばね秤りを用いて5個の物体の重さを10回量る。次の2通りの測定法を考える。
(a) 1つずつそれぞれ2回量る。
(b) 2つずつ10通り$({}_5C_2=10)$の異なる組合せについて各1回量る。
測定誤差は互いに独立に平均0、分散$\sigma^2\ $の正規分布に従うとして、
それぞれの測定法につき、以下の各問に答えよ。
(1) 物体の重さの真値を$\theta=(\theta_1,\dots,\theta_5)'$、測定値を$x=(x_1,\dots,x_{10})'$、
計画行列を$X$、誤差を$\epsilon=(\epsilon_1,\dots,\epsilon_{10})'$として、線形モデル$x=X\theta+\epsilon$を仮定する。
$\theta_i$の最小二乗推定量を$\hat{\theta}_i(i=1,\dots,5)$として、それを導く正規方程式を書け。
(2) $\sigma^2$は未知として、帰無仮説「$H_0:\theta_1=\cdots=\theta_5$」を
対立仮説「$H_1:\theta_1=\cdots=\theta_5$でない」に対して検定するためのF統計量を求めよ。
(3) 対立仮説$H_1$の下でのF統計量の非心度$\lambda$を求めよ。
ただし、$Y_1$および$Y_2$が互いに独立にそれぞれ非心度$\lambda$で自由度$f_1$の非心カイ二乗分布、
および自由度$f_2$のカイ二乗分布に従うとき、
F統計量$F=(Y_1/f_1)/(Y_2/f_2)$の非心度は$\lambda$であると定義される。
解答例
問題では測定値を表す記号として$x$を使用しているが、ここでは代わりに$Y$を使用する。
(1) 最小二乗推定量の正規方程式
$H_1$と$H_0$の下での線形モデルと正規方程式は以下の通り。$n=10,r=5$とする。
すべての要素が1で長さが$r$の縦ベクトルを$L$とし、$X_0=XL$とする。$\theta_0$は長さ1のベクトル。
\begin{align*}
H_1:Y&=X\theta+\epsilon,\hat{\theta}=(X'X)^{-1}X'Y,\dim(X)=(n,r)\\
H_0:Y&=X_0\theta_0+\epsilon,\hat{\theta}_0=(X_0'X_0)^{-1}X_0'Y,\dim(X_0)=(n,1)
\end{align*}
(2) F統計量の定義
尤度比検定を適用する。線形モデルの最大尤度は$(2\pi\hat{\sigma}^2)^{-\frac{n}{2}}\exp(-\frac{n}{2})$となる。
ここで、$\hat{\sigma}^2$は測定誤差の分散の最尤推定量であり、$H_1:\frac{R}{n}$および$H_0:\frac{S}{n}$となる。
従って尤度比は、それぞれのモデルによる残差平方和の比$\frac{S}{R}$で表現できる。
$H_0$のモデルは$H_1$のモデルに制約を追加したものなので、常に$S\ge R$となる。
そこで$\frac{S}{R}$の代わりに$\frac{S-R}{R}=\frac{Q}{R}$を使用する。
\begin{align*}
R&=\|Y-X\hat{\theta}\|^2\\
S&=\|Y-X_0\hat{\theta}_0\|^2\\
Q&=S-R=\|X\hat{\theta}-X_0\hat{\theta}_0\|^2
\end{align*}
さらに既知のF分布に従うように、自由度を考慮してF統計量を以下のように定義する。
\begin{align*}
F&=\frac{Q/(r-1)}{R/(n-r)}\\
R&=\|Y-X\hat{\theta}\|^2=Y'(I-X(X'X)^{-1}X')Y=Y'AY\\
Q&=Y'(X(X'X)^{-1}X'-X_0(X_0'X_0)^{-1}X_0')Y=Y'A_0Y
\end{align*}
ここで$A,A_0$を以下のように定義した。
\begin{align*}
A&=I-X(X'X)^{-1}X'\\
A_0&=X(X'X)^{-1}X'-X_0(X_0'X_0)^{-1}X_0'
\end{align*}
$A,A_0$はべき等行列$(A^2=A,A_0^2=A_0)$であり、固有値は1または0となる。
(3) $A,A_0$の固有値分解
$A,A_0$の固有ベクトルを求めるために、$R^n$を三つの部分空間に直交直和分解する。
\begin{align*}
R^n&=V_1 \oplus V_2 \oplus V_3\\
V_1&=\{x|X \perp x\},\dim(V_1)=n-r,V_1=\langle u_{r+1},\dots,u_n \rangle\\
V_2&=\{x|X \parallel x, X_0 \perp x\},\dim(V_2)=r-1,V_2=\langle u_2,\dots,u_r \rangle\\
V_3&=\{x|X_0 \parallel x\},\dim(V_3)=1,V_3=\langle u_1 \rangle,
\quad u_i'u_j=\delta_{ij}
\end{align*}
ここで、$X \perp x\ $は、$x$が$X$を構成する全ての列ベクトルと直交することを表し、
$X \parallel x$は、$x$が$X$を構成する列ベクトルの一次結合で表現できることを表す。
$u_i(i=1,\dots,n)$は、それぞれの部分空間の正規直交基底であり、全体では$R^n$の正規直交基底を成す。
$X_0 \parallel x$なら、$X \parallel x\ $となることに注意。
$R^n=V_1\oplus V_1^c$と分解して、さらに$V_1^c=V_2\oplus V_3$と分解している。
$x \in V_1 \rightarrow Ax=x$であり、$AX=0$であることから、$x \in V_2\oplus V_3 \rightarrow Ax=0\ $である。
同様に、$x \in V_2 \rightarrow A_0x=x\ $であり、$x \in V_1\oplus V_3 \rightarrow A_0x=0\ $である。
$V_1$を生成する$u_{r+1},\dots,u_n$は、$A$の固有値1に対応する固有ベクトルを構成しており、
$V_2$を生成する$u_2,\dots,u_r$は、$A_0$の固有値1に対応する固有ベクトルを構成している。
従って対称行列の固有値分解により、以下の等式が成立する。
\begin{align*}
A=\sum_{i=r+1}^n u_iu_i',\quad
A_0=\sum_{i=2}^r u_iu_i'
\end{align*}
(4) 帰無仮説$H_0$の下でのF統計量の分布
帰無仮説$H_0$の下での$R,Q,F$の分布は、$Y=X_0\theta_0+\epsilon$と表現されることから、以下のようになる。
パラメータの真の値$\theta_0$は消去されることに注意。
\begin{align*}
H_0:R&=(\theta_0'X_0'+\epsilon')(I-X(X'X)^{-1}X')(X_0\theta_0+\epsilon)\\
&=\epsilon'(I-X(X'X)^{-1}X')\epsilon=\epsilon' A \epsilon\\
&=\textstyle\sum_{i=r+1}^n \epsilon_i^2\sim\sigma^2\chi^2_{n-r}\\
H_0:Q&=(\theta_0'X_0'+\epsilon')(X(X'X)^{-1}X'-X_0(X_0'X_0)^{-1}X_0')(X_0\theta_0+\epsilon)\\
&=\epsilon'(X(X'X)^{-1}X'-X_0(X_0'X_0)^{-1}X_0')\epsilon=\epsilon' A_0 \epsilon\\
&=\textstyle\sum_{i=2}^r \epsilon_i^2\sim\sigma^2\chi^2_{r-1}\\
H_0:F&=\frac{Q/(r-1)}{R/(n-r)}\sim F(r-1,n-r)
\end{align*}
ここで、$\epsilon_i=u_i'\epsilon$と定義した。
定義から、$\epsilon_i$は正規分布に従い、$E[\epsilon_i]=0,E[\epsilon_i\epsilon_j]=0(i\neq j),Var(\epsilon_i)=\sigma^2\ $である。
$H_0:R$は、自由度$n-r$の$\chi^2$分布に従い、
$H_0:Q$は、自由度$r-1$の$\chi^2$分布に従う。
$H_0:F$は、自由度$(r-1,n-r)$のF分布に従う。
(5) 対立仮説$H_1$の下でのF統計量の分布
同様に、対立仮説$H_1$の下での$R,Q,F$の分布、および非心度$\lambda$は以下のようになる。
パラメータの真の値$\theta$が一部残ることに注意。
\begin{align*}
H_1:R&=(\theta'X'+\epsilon')(I-X(X'X)^{-1}X')(X\theta+\epsilon)\\
&=\epsilon'(I-X(X'X)^{-1}X')\epsilon=\epsilon' A \epsilon\\
&=\textstyle\sum_{i=r+1}^n \epsilon_i^2\sim\sigma^2\chi^2_{n-r}\\
H_1:Q&=(\theta'X'+\epsilon')(X(X'X)^{-1}X'-X_0(X_0'X_0)^{-1}X_0')(X\theta+\epsilon)\\
&=\textstyle\sum_{i=2}^r (u_i'X\theta+\epsilon_i)^2
\sim\sigma^2\chi^2_{r-1}(\lambda=\sum_{i=2}^r (u_i'X\theta)^2/\sigma^2)\\
H_1:F&=\frac{Q/(r-1)}{R/(n-r)}\sim F(r-1,n-r,\lambda)\\
\lambda&=\textstyle\sum_{i=2}^r (u_i'X\theta)^2/\sigma^2
=\theta'X'A_0X\theta/\sigma^2\\
&=\theta'X'(X(X'X)^{-1}X'-X_0(X_0'X_0)^{-1}X_0')X\theta/\sigma^2\\
&=\theta'(X'X-X'X_0(X_0'X_0)^{-1}X_0'X)\theta/\sigma^2
\end{align*}
$H_1:R$は、自由度$n-r$の$\chi^2$分布に従い、
$H_1:Q$は、自由度$r-1$、非心度$\lambda$の非心$\chi^2$分布に従う。
$H_1:F$は、自由度$(r-1,n-r)$、非心度$\lambda$の非心F分布に従う。
非心度の計算
非心度$\lambda=\theta'(X'X-X'X_0(X_0'X_0)^{-1}X_0'X)\theta/\sigma^2$を計算する。
ここで、$(X'X-X'X_0(X_0'X_0)^{-1}X_0'X)L=0$となるため、
$\lambda=(\theta-\bar{\theta})'(X'X-X'X_0(X_0'X_0)^{-1}X_0'X)(\theta-\bar{\theta})/\sigma^2$として良い。
ただし、$\bar{\theta}=\frac{1}{r}\sum_{i=1}^r\theta_i$とし、$\theta-\bar{\theta}L$を$\theta-\bar{\theta}$と略記する。
測定法(a)の場合
\begin{align}
X&=\begin{pmatrix}
1 & 0 & 0 & 0 & 0\\
1 & 0 & 0 & 0 & 0\\
0 & 1 & 0 & 0 & 0\\
0 & 1 & 0 & 0 & 0\\
0 & 0 & 1 & 0 & 0\\
0 & 0 & 1 & 0 & 0\\
0 & 0 & 0 & 1 & 0\\
0 & 0 & 0 & 1 & 0\\
0 & 0 & 0 & 0 & 1\\
0 & 0 & 0 & 0 & 1\end{pmatrix},
X'X=\begin{pmatrix}
2 & 0 & 0 & 0 & 0\\
0 & 2 & 0 & 0 & 0\\
0 & 0 & 2 & 0 & 0\\
0 & 0 & 0 & 2 & 0\\
0 & 0 & 0 & 0 & 2
\end{pmatrix},
X_0=1_{10},X'X_0=2_5
\end{align}
ここで、$1_{10}$はすべての要素が1で長さが10の縦ベクトルを表す。$2_{5}$も同様。
従って、$(\theta-\bar{\theta})'X'X_0=0$となるので、$\lambda=2\sum_{i=1}^5(\theta_i-\bar{\theta})^2/\sigma^2$となる。
測定法(b)の場合
\begin{align}
X&=\begin{pmatrix}
1 & 1 & 0 & 0 & 0\\
1 & 0 & 1 & 0 & 0\\
1 & 0 & 0 & 1 & 0\\
1 & 0 & 0 & 0 & 1\\
0 & 1 & 1 & 0 & 0\\
0 & 1 & 0 & 1 & 0\\
0 & 1 & 0 & 0 & 1\\
0 & 0 & 1 & 1 & 0\\
0 & 0 & 1 & 0 & 1\\
0 & 0 & 0 & 1 & 1\end{pmatrix},
X'X=\begin{pmatrix}
4 & 1 & 1 & 1 & 1\\
1 & 4 & 1 & 1 & 1\\
1 & 1 & 4 & 1 & 1\\
1 & 1 & 1 & 4 & 1\\
1 & 1 & 1 & 1 & 4
\end{pmatrix},
X_0=2_{10},X'X_0=8_5
\end{align}
同様に、$\lambda=3\sum_{i=1}^5(\theta_i-\bar{\theta})^2/\sigma^2$となる。
おわりに
- こんな問題を出題されても、解答を記すには余白が狭すぎる。