Rにたくさんさいころを振らせて、さいころの分散を求めてみる

前回記事Rにさいころをたくさん振らせて、さいころの平均を求めてみる - わかばめにっきでは、Rの擬似乱数を使って標本平均が母平均の一致推定量であることを示しました。今回は後で述べる不偏標本分散という量が母分散の不偏推定量であることを調べていきます。

ところで母平均と母分散の定義ってなんだっけ？

確率変数 $X$ が離散型、つまり実現値の個数が $M$ 個であって、それぞれの確率が $P(X = x_i) = p_i, \ i=1,\cdots,M$ のときを考えます。このとき、母平均 $\mu$ 、母分散 $\sigma^2$ はそれぞれ

$\begin{align} \mu &= \sum_{i=1}^M x_i p_i\\ \sigma^2 &= \sum_{i=1}^M (x_i-\mu)^2p_i \end{align}$

により与えられます。一方で確率変数 $X$ の連続型のときは、確率密度を $f(x)$ とすると、母平均 $\mu$ 、母分散 $\sigma^2$ はそれぞれ

$\begin{align} \mu &= \int_{\mathbb{R}} xf(x)\,dx\\ \sigma^2 &= \int_{\mathbb{R}} (x-\mu)^2f(x)\,dx \end{align}$

により与えられます。母平均は確率変数の期待値を、母分散はどれくらい散らばりやすい確率変数かを表しています。

理想的なさいころの母平均と母分散はどうなる？

確率変数 $X$ が理想的な確率分布 $P(X = i) = 1/6, \ i=1,\cdots,M$ に従うとき、そのまま計算すれば、

$\begin{align} \mu &= \sum_{i=1}^6 i \times\frac16\\ &= \frac72\\ \sigma^2 &= \sum_{i=1}^6 \left(i-\frac72\right)^2\times\frac16\\ &= \frac16\left(\frac{25}{4}+\frac{9}{4}+\frac{1}{4}+\frac{1}{4}+\frac{9}{4}+\frac{25}{4}\right)\\ &= \frac{35}{12} \end{align}$

となり母平均は 3.5、母分散はおおむね 2.917 くらいということがわかります。わたしたちはこれから、この数値を推定していくにあたり、不偏推定量という概念を導入します。

推定量が不偏ってどういうことさ？

不偏性 | 統計用語集 | 統計WEBによると、不偏性とは「その推定量が平均的に過大にも過小にも母数を推定しておらず、推定量の期待値が母数に等しいこと」です。ですが初めてこれを読んでも意味を理解するのは難しいのではないでしょうか。
そこで、「さいころを５個しか振れない」という条件のもとで、母分散を調べることを試みましょう。まずは前回と同様にさいころを計算機に振らせます。

> x = as.integer(runif(n,1,7))
> x
[1] 2 1 5 1 2

わたしたちは5個の数字の羅列 $x_1, \cdots ,x_5$ しか知らないので、母分散の定義に含まれる母平均すら現時点では推定の対象です。ひとまず標本の平均値 $\bar{x} = \frac{x_1+\cdots+x_5}{5}$ とおき、二乗誤差の平均値

$\displaystyle\dfrac15\sum_{i=1}^5\left( x_i -\bar{x} \right)^2$

を用いて分散を推定していきましょう。

> sum((x-mean(x))^2)/5
[1] 2.16

このステップを何度か繰り返してみましょう。

> x = as.integer(runif(n,1,7))
> sum((x-mean(x))^2)/5
[1] 3.44
> x = as.integer(runif(n,1,7))
> sum((x-mean(x))^2)/5
[1] 1.76
> x = as.integer(runif(n,1,7))
> sum((x-mean(x))^2)/5
[1] 3.76
> x = as.integer(runif(n,1,7))
> sum((x-mean(x))^2)/5
[1] 1.76
> x = as.integer(runif(n,1,7))
> sum((x-mean(x))^2)/5
[1] 1.84
> x = as.integer(runif(n,1,7))
> sum((x-mean(x))^2)/5
[1] 2.16

当然ですがブレまくりますね。それに、こころなしか小さめに出ているような気がします。そこで、数学を使ってこの推定の筋がいいかを考えましょう。つまり、「5個の標本から計算した標本分散」を確率変数とみなして、確率変数が持つ性質の1つである平均を求め、それがもとの母分散に一致しているかを調べます。これが一致していたらそんなに悪い指標ではないと言えそうですよね。

また、話を一般化するために標本の大きさを5でなく $n$ とします。すなわち $X$ と同分布な確率変数 $X_i,\, i=1,\cdots ,n$ を各々独立であるとし、標本平均と標本誤差平均をそれぞれ

$\begin{align} \bar{X} &= \frac{1}{n}\sum_{i=1}^n X_i\\ S^2 &= \frac{1}{n}\sum_{i=1}^n (X_i- \bar{X})^2 \end{align}$

と定めます。それでは、母分散を推定するために定めた標本誤差平均の平均値である $E[S^2]$ を求めましょう。ここでの計算は $X$ は連続型、すなわち確率密度 $f(x)$ を持つものとします。このとき独立性の仮定から、 $i \neq j$ に対して

$\begin{equation} E[X_i^2]=\mu_2,\quad E[X_i X_j]=\mu_1^2 \end{equation}$

が成立します。ただしここで $\mu_1$ 、 $\mu_2$ はそれぞれ確率変数 $X$ の1次モーメント、2次モーメント：

$\begin{align} \mu_1&= \int_{\mathbb{R}} xf(x)\,dx\ (=\mu)\\ \mu_2&= \int_{\mathbb{R}} x^2f(x)\,dx\ (=\sigma^2 +\mu^2) \end{align}$

を表すものとします。実際、 $g(X_1,\cdots,X_n)$ の確率密度は $f(x_1)f(x_2)\cdots f(x_n)$ なので、例えば

$\begin{align} E[X_1^2]&=\int_\mathbb{R}\int_\mathbb{R}\cdots\int_\mathbb{R} x_1^2 f(x_1)f(x_2)\cdots f(x_n)\, dx_1 dx_2 \cdots dx_n \\ &= \int_\mathbb{R} x_1^2f(x_1)dx_1 \int_\mathbb{R}f(x_2)dx_2 \cdots\int_\mathbb{R} f(x_n)dx_n\\ &= \mu_2 \times 1\times \cdots \times 1\\ E[X_1 X_2]&=\int_\mathbb{R}\int_\mathbb{R}\cdots\int_\mathbb{R} x_1 x_2 f(x_1)f(x_2)\cdots f(x_n)\, dx_1 dx_2 \cdots dx_n \\ &= \int_\mathbb{R} x_1f(x_1)dx_1 \int_\mathbb{R} x_2f(x_2)dx_2 \int_\mathbb{R}f(x_3)dx_3 \cdots\int_\mathbb{R} f(x_n)dx_n\\ &= \mu_1 \times \mu_1 \times 1\times \cdots \times 1\\ \end{align}$

となることから理解できます。また、 $\displaystyle \sum_{i=1}^n X_i = n\bar{ X}$ を用いると、 $i =1,\cdots ,n$ に対して

$\displaystyle \sum_{i=1}^n( X_i -\bar{ X})^2= \displaystyle \sum_{i=1}^n( X_i^2 -2\bar{ X} X_i +\bar{ X }^2) = \displaystyle \sum_{i=1}^n X_i ^2- n\bar{ X}^2$

に注意することにより、

$E[ S^2]=\dfrac{1}{n}\displaystyle \sum_{i=1}^n E[ X_i ^2]- E[ \bar{ X}^2]$

を得ます。第1項は独立の仮定から $\mu_2$ に一致します。第2項は相異なる $i$ 、 $j$ の組み合わせは $n(n-1)$ 個あることにより

$\begin{align}E[ \bar{ X}^2] &= \dfrac{1}{n^2}\sum_{i=1}^n E\left[ X_i\sum_{j=1}^n X_j\right]\\ &= \dfrac{1}{n^2}\sum_{i=1}^n E\left[ X_i\left(X_i+\sum_{j\neq i} X_j\right)\right]\\ &= \dfrac{1}{n^2}\sum_{i=1}^n E\left[ X_i^2\right] + \dfrac{1}{n^2} \sum_{i\neq j} E\left[ X_iX_j\right] \\ &= \dfrac{n}{n^2} \mu_2 +\frac{n(n-1)}{n^2} \mu_1^2 \end{align}$

となり、分散の定義から $\mu_2 - \mu_1^2 = \sigma^2$ により

$\begin{align} E[ S^2]&=\dfrac{1}{n}\displaystyle \sum_{i=1}^n E[ X_i ^2]- E[ \bar{ X}^2]\\ &= \mu_2 -\left( \frac{1}{n} \mu_2 + \frac{n-1}{n} \mu_1^2\right)\\ &= \frac{n-1}{n} \sigma^2 \end{align}$
を得ます。したがって $S^2$ は平均して母分散の値より $\frac{n-1}{n}$ 倍だけ小さく現れてしまいます*1。

実際にデータで確かめてみよう。

確率変数 $S^2$ の平均値を求めるために、 $S^2$ の実現値を10000回求めてみます。

ss = c()

M = 10000
n = 5
for (i in 1:M){
x = as.integer(runif(n,1,7))
m = mean(x)
y = (x - m)^2
S2 = sum(y)/5
ss <- append(ss, S2)
}
hist(ss)

f:id:wakabame:20180424213547p:plain
$X$ は離散値なので、でこぼこなヒストグラムが現れました*2。この平均を求めてみましょう。

> mean(ss)
[1] 2.322912

たしかに真の平均とくらべて 4/5 倍くらいの大きさになっていそうです。

したがって、 $n$ 個の数値データを持つ標本から母分散を求めるには、
$\begin{equation} s^2 =\dfrac{1}{n-1}\displaystyle\sum_{i=1}^n (X_i- \bar{X})^2 \end{equation}$
なる確率変数を推定量とすればよいことがわかりました。このように定めた $s^2$ を不偏標本分散と呼びます。
さらに、 $E[s^2] = \sigma^2$ のように期待値を取れば母数に一致するような推定量を、その母数の不偏推定量であると呼びます。
例えば標本平均 $\bar{X}$ は母平均 $\mu$ の不偏推定量です。試しにこれを示してみてください。