Help us understand the problem. What is going on with this article?

非心カイ二乗分布に関する出題

More than 1 year has passed since last update.

はじめに

  • 統計検定一級の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$となる。

おわりに

  • こんな問題を出題されても、解答を記すには余白が狭すぎる。
fred55
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away