WWWを検索
サイト内を検索

Rによる適合度のχ二乗検定

Check Tweet

何かと使われるχ²検定の理屈を抑えておきましょう。

適合度のχ二乗検定についての数理的な理屈を確認した上で、Rによる実践方法を確認しましょう。

1 χ二乗分布の導出

自由度と同じ数の標準正規分布に従う確率変数の二乗和の分布がχ二乗分布です。

1.1 自由度\(1\)

標準正規分布に従う確率変数\(x\)が\((-\infty, \sqrt{z})\)に入る確率は、

\[ \int_{-\infty}^{\sqrt{z}} f(x) dx = \int_{-\infty}^{\sqrt{z}} \frac{1}{\sqrt{2\pi}}\exp\big(-\frac{x^2}{2}\big) dx \]

となります。

\(y = x^2\)として、\(y\)が非負なので2点\(-x\)と\(x\)が1点\(y\)に写されることに注意し、\(x\)と\(y\)を変数変換をします。

\[ \frac{dy}{dx} = 2x = 2\sqrt{y} \]

\[ \int_0^{z} f(y) dy = \int_0^{z} 2\frac{1}{\sqrt{2\pi}}\exp\big(-\frac{y}{2}\big)\frac{1}{2\sqrt{y}} dy \]

\(\Gamma(1/2) = \sqrt{\pi}\)に注意して整理すると、

\[ \int_0^{z} f(y) dy = \int_0^{z} \frac{1}{2^\frac{1}{2} \Gamma(1/2)} y^{-\frac{1}{2}} \exp\big(-\frac{y}{2}\big) dy \]

と自由度\(1\)のχ²分布に従う確率変数が、\([0, z)\)の間をとる確率になります。

1.2 自由度\(n\)

自由度\(1\)のχ²分布は\(k=1/2\), \(\theta=2\)のガンマ関数で、後述するガンマ関数の再生性から直ちに自由度\(n\)のχ²分布に従う確率変数の形状を言えますが、自由度\(2\)以上のχ二乗分布について畳み込み計算で示します。

\[ y_{n-1} = x_1^2 + \cdots + x_{n-1}^2 \]

が、χ二乗分布に従うとします。\(x_i\ (i\in[1, n-1])\)は、標準正規分布に従う確率変数です。このとき\(y_{n-1} + y_1\)の確率密度関数は、

\[\begin{align} f_n(y) &= \int_0^y f_{n-1}(t) f_1(x-t) dt \\ &= \frac{\exp\big( -\frac{y}{2} \big)}{2^\frac{n}{2} \Gamma\big(\frac{n-1}{2} \big)\sqrt{\pi}} \int^y_0 t^\frac{n-3}{2} (y-t)^{-\frac{1}{2}}dt \end{align}\]

となり、\(u=t/x\)で変数変換をすると

\[\begin{align} &= \frac{\exp\big( -\frac{y}{2} \big) y^{\frac{n-3}{2} - \frac{1}{2} + 1} }{2^\frac{n}{2} \Gamma\big(\frac{n-1}{2} \big) \Gamma\big(\frac{1}{2} \big)} \int^1_0 u^\frac{n-3}{2}(1-u)^{-\frac{1}{2}}du \\ &= \frac{\exp\big( -\frac{y}{2} \big) y^{\frac{n-3}{2} - \frac{1}{2} + 1} }{2^\frac{n}{2} \Gamma\big(\frac{n-1}{2} \big) \Gamma\big(\frac{1}{2} \big)} B \bigg( \frac{n-1}{2}, \frac{1}{2} \bigg) \\ &= \frac{\exp\big( -\frac{y}{2} \big) y^{\frac{n-3}{2} - \frac{1}{2} + 1} }{2^\frac{n}{2} \Gamma\big(\frac{n-1}{2} \big) \Gamma\big(\frac{1}{2} \big)} \frac{\Gamma\bigg(\frac{n-1}{2}\bigg)\Gamma\bigg(\frac{1}{2}\bigg)}{\Gamma\bigg(\frac{n}{2}\bigg)} \\ &= \frac{1}{2^\frac{n}{2} \Gamma\big( \frac{n}{2} \big)} \exp\big( -\frac{y}{2} \big) y^{\frac{n}{2} - 1} \end{align}\]

\(k=n/2\), \(\theta=2\)のガンマ分布になりました。

1.3 ガンマ分布のモーメント母関数と再生性

中心極限定理を使う場合、モーメント母関数もしくは特性関数の性質を認めているわけで、そちらから計算した方が楽です。

\[\begin{align} M_X(t) &= \int \exp(tx) \frac{x^{k-1} \exp\big(-\frac{x}{\theta}\big)}{\Gamma(k)\theta^k} dx \\ &= (1-\theta t)^{-k} \int \frac{x^{k-1}\exp\big(-\frac{1-\theta t}{\theta} \big)}{\Gamma(k)\big(\frac{\theta}{1-\theta t}\big)^k} dx \end{align}\]

積分の内側の\(\theta/(1-\theta t)\)を\(z\)と置くと、ガンマ関数の確率密度関数になり、確率密度関数の定義から積分の値が\(1\)になることが分かります。

よってガンマ分布のモーメント母関数は、

\[ M_X(t) = (1-\theta t)^{-k} \]

となります。

よってパラメーターが\(\theta\)と\(k=k_1\)、パラメーターが\(\theta\)と\(k=k_2\)のガンマ分布に従う2つの確率変数の和が従う確率分布のモーメント母関数は、

\[ M_X^{k=k_1}(t)M_X^{k=k_2}(t) = (1-\theta t)^{-k_1-k_2} = M_X^{k=k_1+k_2} \]

となり、パラメーターが\(\theta\)と\(k1+k2\)のガンマ分布のモーメント母関数になるので、モーメント母関数の性質から\(\theta\)が同一であれば再生性があることが分かります。

2 行列式の補題

恒等行列\(I\)と、正則行列\(A\in\mathbb{R}^{n \times n}\)と、列ベクトル\(\mu, \mu', \nu \in \mathbb{R}^{n}\)を考えます。

\[\begin{align} &\det(I + \mu'{}^t\!\nu) \\ &= \det\begin{pmatrix} I + \mu'{}^t\!\nu & \mu' \\ 0 & 1 \end{pmatrix} \\ &= \det\begin{pmatrix} I & 0 \\ {}^t\!\nu & 1 \end{pmatrix}\det\begin{pmatrix} I + \mu'{}^t\!\nu & \mu' \\ 0 & 1 \end{pmatrix}\det\begin{pmatrix} I & 0 \\ -{}^t\!\nu & 1 \end{pmatrix} \\ &= \det\begin{pmatrix} I & u \\ 0 & 1+{}^t\!\nu\mu' \end{pmatrix} \\ &= (1 + {}^t\!\nu\mu') \det(I) \end{align}\]

さらに、

\[\begin{align} &\det(A + \mu'{}^t\!\nu) \\ &= \det(A(I + A^{-1}\mu'{}^t\!\nu)) \\ &= \det(A)\det(I + A^{-1}\mu'{}^t\!\nu) \end{align}\]

\(\mu' = A\mu\)として補題を適応すると

\[ \det(A + \mu{}^t\!\nu) = \det(A)(1 + {}^t\!\nu\mu) \]

3 シャーマン・モリソンの公式

初見だと使い道が無さそうな気もするのですが、BFGS法の導出に使われていたりと意外に有能なのがSherman–Morrison formulaです。具体的には行列\(A\in\mathbb{R}^{n \times n}\)と、列ベクトル\(\mu, \nu \in \mathbb{R}^{n}\)に対して、

1. \(A+\mu{}^t\!\nu\)が正則 ⇔ \(1+{}^t\!\nu A^{-1}\mu \ne 0\)
2. (1)のとき、\((A + \mu{}^t\!\nu)^{-1} = A^{-1} - \frac{A^{-1}\mu{}^t\!\nu A^{-1}}{1+{}^t\!\nu A^{-1}\mu}\)

の二つを主張しています。

(1)の\(1 + {}^t\!\nu A^{-1}\mu \ne 0\) ⇒ \(A + \mu{}^t\!\nu\)が正則は、(2)の右辺に\(A + \mu{}^t\!\nu\)を乗じると、

\[\begin{align} &(A + \mu{}^t\!\nu)(A^{-1} - \frac{A^{-1}\mu{}^t\!\nu A^{-1}}{1+{}^t\!\nu A^{-1}\mu}) \\ &= (A^{-1} - \frac{A^{-1}\mu{}^t\!\nu A^{-1}}{1+{}^t\!\nu A^{-1}\mu})(A + \mu{}^t\!\nu) \\ &= AA^{-1} + \mu{}^t\!\nu A^{-1} - \frac{AA^{-1}\mu{}^t\!\nu A^{-1} + \mu{}^t\!\nu A^{-1}\mu{}^t\!\nu A^{-1}}{1 + {}^t\!\nu A^{-1}\mu} \\ &= I + \mu{}^t\!\nu A^{-1} - \frac{\mu(1 + {}^t\!\nu A^{-1}\mu){}^t\!\nu A^{-1}}{1 + {}^t\!\nu A^{-1}\mu} \\ &= I + \mu{}^t\!\nu A^{-1} - \mu{}^t\!\nu A^{-1} \\ &= I \end{align}\]

と、\(A^{-1} - \frac{A^{-1}\mu{}^t\!\nu A^{-1}}{1+{}^t\!\nu A^{-1}\mu}\)が\(A + \mu{}^t\!\nu\)の逆行列となることから分かります。逆行列があれば正則です。同時に(2)が示されます。

\(1 + {}^t\!\nu A^{-1}\mu = 0\)だとすると、

\[\begin{align} &\det(A + \mu{}^t\!\nu) \\ &= \det(A(I + (A^{-1}\mu){}^t\!\nu)) \\ &= \det(A)(1 + {}^t\!\nu A^{-1}\mu) \\ &= 0 \end{align}\]

となり、\(A + \mu{}^t\!\nu\)は逆行列が無いので非正則となります。

\(k\)種類の因子にそれぞれ観測値がある適合度の検定を考えます。\(X_i\)を多項分布に従う観測値のベクトルで、ある要素だけが\(1\)で、他は\(0\)となります。集計表に記入された生起数は\(\sum_{i=1}^n X_i\)です。

4 適合度検定の統計量がχ二乗分布に従うことの証明

証明方法は多くある¹のですが、シャーマン・モリソンの公式を使ったものを確認しましょう。

帰無仮説が示すその期待値\(E[X_i]\)を\(p\)とおきます。\(p\)はベクトルで、その要素を\(p_1, p_2, \cdots, p_k\)と書きます。中心極限定理から

\[ \sqrt{n}(\bar{X}_n - p) \xrightarrow[n \to \infty]{\mathrm{d}} \mathcal{N}_k(0, \Sigma) \]

となり、k次元の正規分布で近似できることが分かります。ここで\(\bar{X}_n\)は観測数\(n\)の観測値の平均ベクトルであり、\(\Sigma\)は分散共分散行列です。多項分布ですから、分散が\(p_i(1-p_i)\)で、共分散は\(-p_ip_j\ (i \ne j)\)となります。また、\(\sum_{i=1}^k p_i=1\)なので、\(\Sigma\)の各行と各列は線形従属しており、逆行列をつくれません。

観測値から\(k\)番目の要素を除外したベクトルを\(X_i^*\)，期待値から\(k\)番目の要素を除外したベクトルを\(p^*\)と置きます。

\[ X_i^* = \begin{pmatrix} X_{i1} \\ X_{i2} \\ \vdots \\ X_{ik-1} \end{pmatrix},\qquad p^* = \begin{pmatrix} p_1 \\ p_2 \\ \vdots \\ p_{k-1} \end{pmatrix} \]

\(X_i^*\)に対応する分散共分散行列は、

\[\begin{align} \Sigma^* &= \begin{pmatrix} p_1(1-p_1) & -p_1p_2 & \cdots & -p_1p_{k-1} \\ -p_1p_2 & p_2(1-p_2) & \cdots & -p_2p_{k-1} \\ \vdots & \vdots & \ddots & \vdots \\ -p_1p_{k-1} & -p_2p_{k-1} & \cdots & p_{k-1}(1-p_{k-1}) \end{pmatrix} \\ &= \begin{pmatrix} p_1 & 0 & \cdots & 0 \\ 0 & p_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & p_{k-1} \end{pmatrix} - \begin{pmatrix} p_1p_1 & p_1p_2 & \cdots & p_1p_{k-1} \\ p_1p_2 & p_2p_2 & \cdots & p_2p_{k-1} \\ \vdots & \vdots & \ddots & \vdots \\ p_1p_{k-1} & p_2p_{k-1} & \cdots & p_{k-1}p_{k-1} \end{pmatrix}\\ &= \begin{pmatrix} p_1 & 0 & \cdots & 0 \\ 0 & p_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & p_{k-1} \end{pmatrix} - \begin{pmatrix} p_1 \\ p_2 \\ \vdots \\ p_{k-1} \end{pmatrix}\begin{pmatrix} p_1 & p_2 & \cdots & p_{k-1} \end{pmatrix} \\ &= \mathrm{diag}(p) - p^*{}^t\!p^* \end{align}\]

となります。

\[ 1 - {}^t\!p^* \mathrm{diag}(p^*)^{-1} p = \sum_{i=1}^k p_i - \sum_{i=1}^{k-1} p_i = p_k \]

\[ \mathrm{diag}(p^*)^{-1} p = \begin{pmatrix} 1 \\ 1 \\ \vdots \\ 1 \end{pmatrix} \]

上の二つの式に注意して、\((\mathrm{diag}(p) - p^*{}^t\!p^*)^{-1}\)に、\(A=\mathrm{diag}(p^*)\), \(\mu=\nu=p^*\)とおいて、前節で導出したSherman–Morrison formulaを用いると、

\[ (\Sigma^*)^{-1} = \begin{pmatrix} 1/p_1 & 0 & \cdots & 0 \\ 0 & 1/p_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1/p_{k-1} \end{pmatrix} + \frac{1}{p_k}\begin{pmatrix} 1 & 1 & \cdots & 1 \\ 1 & 1 & \cdots & 1 \\ \vdots & \vdots & \ddots & \vdots \\ 1 & 1 & \cdots & 1 \end{pmatrix} \]

が得られます。

ここからχ二乗統計量が何を意味するのか調べていきましょう。集計表に記入される因子\(i\)の観測値は\(n\bar{X}_i\)、その期待値は\(np_i\)になりますが、そのχ二乗統計量を変形します。

\[\begin{align} \chi^2 &= \sum_{i=1}^k \frac{(n\bar{X}_i - np_i)^2}{np_i} \\ &= n \sum_{i=1}^k \frac{(\bar{X}_i - p_i)^2}{p_i} \\ &= n \bigg[ \sum_{i=1}^{k-1} \frac{(\bar{X}_i - p_i)^2}{p_i} + \frac{(\bar{X}_k - p_k)^2}{p_k} \bigg] \\ &= n \bigg[ \sum_{i=1}^{k-1} \frac{(\bar{X}_i - p_i)^2}{p_i} + \frac{(\sum_{i=1}^{k-1}(\bar{X}_i - p_i))^2}{p_k} \bigg] \end{align}\]

最後の変形は\(\sum_{i=1}^{k}\bar{X}_i = 1\), \(\sum_{i=1}^{k}p_i = 1\)から、\(\sum_{i=1}^{k}(\bar{X}_i - p_i)=0\)となることを使いました。

行列を使って書き直します。

\[\begin{align} \chi^2 &= n\ {}^t\!(\bar{X}^* - p^*) (\Sigma^*)^{-1} (\bar{X}^* - p^*) \end{align}\]

二次形式になりました。

中心極限定理から、上の式の平方根である\(\sqrt{n} (\Sigma^*)^{-1/2} (\bar{X}^* - p^*)\)は\(\mathcal{N}_{k-1}(0, I_{k-1})\)に分布収束します。よって、

\[\begin{align} \chi^2 \xrightarrow[n \to \infty]{\mathrm{d}} {}^t\!\mathcal{N}_{k-1}(0, I_{k-1}) \mathcal{N}_{k-1}(0, I_{k-1}) \end{align}\]

となり、χ二乗統計量は\(k-1\)次元の標準正規分布の二乗、つまり自由度\(k-1\)のχ二乗分布に収束することが分かります。これで適合度のχ二乗検定を行う前提を確認できました。

5 Rによる適合度検定

適合度検定（goodness-of-fit）を試します。

まずは教科書的に統計量から計算してみます。

# 練習用にデータセットを乱数から生成
set.seed(522)
n <- 25
(X <- rmultinom(n, 1, prob = c(0.1, 0.5, 0.4)))

     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
[1,]    0    0    0    0    0    1    0    0    0     0     0     0
[2,]    0    0    0    0    0    0    0    1    0     1     1     0
[3,]    1    1    1    1    1    0    1    0    1     0     0     1
     [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22]
[1,]     0     1     0     0     0     0     0     0     0     0
[2,]     1     0     1     0     1     1     1     1     0     1
[3,]     0     0     0     1     0     0     0     0     1     0
     [,23] [,24] [,25]
[1,]     0     0     0
[2,]     0     1     1
[3,]     1     0     0

# 合計する
(sX <- apply(X, 1, sum))

[1]  2 12 11

# 自由度を計算
k <- length(sX)
# 帰無仮説を設定
p <- 1/k
# χ二乗統計量を計算
q <- sum((sX - n*p)^2/(n*p))
# p値を計算
pchisq(q, k-1, lower.tail = FALSE)

[1] 0.02625234

次により実践的に、ビルトイン関数を使います。

# χ二乗検定をかける(帰無仮説はp = c(1/3, 1/3, 1/3))
chisq.test(sX)

    Chi-squared test for given probabilities

data:  sX
X-squared = 7.28, df = 2, p-value = 0.02625

# χ二乗検定をかける(帰無仮説はp = c(0.1, 0.5, 0.4))
chisq.test(sX, p = c(0.1, 0.5, 0.4))

Warning in chisq.test(sX, p = c(0.1, 0.5, 0.4)):
カイ自乗近似は不正確かもしれません

    Chi-squared test for given probabilities

data:  sX
X-squared = 0.22, df = 2, p-value = 0.8958

簡単ですね。

---

1. Benhamou and Melot (2018)に7種類紹介されています。↩︎