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

最近は東京大学出版会の統計学入門を用いて、統計学に再入門しています。標本分散の定義に $n-1$ が出てくるのはなぜだっけ？というところからよくわかっていなかったので、何回かに分けてその理由をメモしていきます。
推定とはなにかというところから出発します。若干迂遠に感じられるかもしれませんが、間違い探しだと思ってお付き合いください。

- - -

確率的に変動する要素がある場合、わたしたちはその変動をできるだけ正確に見積もりたいのに対して、厳密に調べつくすことは現実的ではないという問題に遭遇します。
例えばさいころを振ることを考えても、実際にはわたしたちは有限回しかさいころを振るという行動はできません。
したがって、「理想的なさいころであればそれぞれの目が出る確率は $\frac16$ になるはずだ」ということに対して、手元にあるさいころの背後にある確率分布を精密に調べることは不可能です。
そうであれば、実験、調査、測定などによって得られたデータ（これを標本と呼ぶ）のどのような性質を見れば、確率分布の持つ平均や分散などといった特徴（これを母数と呼ぶ）を調べられるでしょうか？また、性質の見方には指標があるでしょうか？ということを考える必要が生じます。
この記事群では一致推定量、不偏推定量という観点から、母平均と母分散を"推定"することを考えていきます。*1

とりあえずは同じ形のさいころを200個振ってみます。有限個とは言え疲れそうなので、計算機にやらせましょう。

> N = 200
> x = as.integer(runif(N,1,7))
> x
  [1] 6 2 5 6 1 1 1 6 2 4 6 1 3 2 3 5 4 6 5 5 4 6 3 5 6
 [26] 5 6 6 5 3 6 4 6 5 5 4 5 3 6 4 4 2 6 2 3 1 6 1 4 5
 [51] 2 5 6 2 3 5 2 1 1 1 5 3 3 3 5 1 3 5 2 2 2 4 2 3 5
 [76] 5 3 3 4 6 6 1 2 6 4 2 3 6 1 3 2 1 3 3 3 4 5 5 2 6
[101] 3 5 5 6 4 4 4 1 3 1 5 2 6 2 1 3 1 3 5 3 4 1 1 6 1
[126] 2 1 3 2 2 4 5 3 3 1 1 2 1 1 5 3 5 5 4 4 3 1 2 2 5
[151] 1 6 1 4 4 6 3 3 4 6 3 2 2 3 1 3 6 3 3 2 3 5 1 5 3
[176] 6 6 6 1 1 4 2 5 6 4 4 2 6 3 4 3 2 1 2 4 2 2 5 5 4

Rの擬似乱数を使って1,2,3,4,5,6 の数字をランダムに200個表示させています。
つまり、 $P(X=k)= \frac16,\quad (k =1,2,\cdots ,6)$ なる確率分布を持つ確率変数 $X$ の実現値を並べてみました。
ただし、統計学とはあくまでこれを推定するという立場ですから、これらの数字がどのように生成されているかはわかっておらず、200個の数字の羅列 $x = (x_n)_{n=1}^{200}$ だけがわかっています。ですので母平均が $\mu =\frac72$ だということはわかっていないということに注意しておきます。

得られたデータから、さいころの平均を推定しよう！

さて、ここから確率変数 $X$ の平均 $E[X]$ を推定しましょう。200個のデータがあるので、標本平均
$\dfrac{x_1 + x_2 +\cdots +x_{200}}{200}$
は良さそうです。標本の大きさに対して標本平均がどう推移していくかを調べたいので、1個目からn個目までのデータの平均値
$y_n = \dfrac{x_1 + x_2 +\cdots +x_{n}}{n}$
の列を考え、その推移を考えていくことにしましょう。

> y = cumsum(x)/c(1:N)
> y
  [1] 6.000000 4.000000 4.333333 4.750000 4.000000
  [6] 3.500000 3.142857 3.500000 3.333333 3.400000
 [11] 3.636364 3.416667 3.384615 3.285714 3.266667
 [16] 3.375000 3.411765 3.555556 3.631579 3.700000
 [21] 3.714286 3.818182 3.782609 3.833333 3.920000
 [26] 3.961538 4.037037 4.107143 4.137931 4.100000
 [31] 4.161290 4.156250 4.212121 4.235294 4.257143
 [36] 4.250000 4.270270 4.236842 4.282051 4.275000
 [41] 4.268293 4.214286 4.255814 4.204545 4.177778
 [46] 4.108696 4.148936 4.083333 4.081633 4.100000
 [51] 4.058824 4.076923 4.113208 4.074074 4.054545
 [56] 4.071429 4.035088 3.982759 3.932203 3.883333
 [61] 3.901639 3.887097 3.873016 3.859375 3.876923
 [66] 3.833333 3.820896 3.838235 3.811594 3.785714
 [71] 3.760563 3.763889 3.739726 3.729730 3.746667
 [76] 3.763158 3.753247 3.743590 3.746835 3.775000
 [81] 3.802469 3.768293 3.746988 3.773810 3.776471
 [86] 3.755814 3.747126 3.772727 3.741573 3.733333
 [91] 3.714286 3.684783 3.677419 3.670213 3.663158
 [96] 3.666667 3.680412 3.693878 3.676768 3.700000
[101] 3.693069 3.705882 3.718447 3.740385 3.742857
[106] 3.745283 3.747664 3.722222 3.715596 3.690909
[111] 3.702703 3.687500 3.707965 3.692982 3.669565
[116] 3.663793 3.641026 3.635593 3.647059 3.641667
[121] 3.644628 3.622951 3.601626 3.620968 3.600000
[126] 3.587302 3.566929 3.562500 3.550388 3.538462
[131] 3.541985 3.553030 3.548872 3.544776 3.525926
[136] 3.507353 3.496350 3.478261 3.460432 3.471429
[141] 3.468085 3.478873 3.489510 3.493056 3.496552
[146] 3.493151 3.476190 3.466216 3.456376 3.466667
[151] 3.450331 3.467105 3.450980 3.454545 3.458065
[156] 3.474359 3.471338 3.468354 3.471698 3.487500
[161] 3.484472 3.475309 3.466258 3.463415 3.448485
[166] 3.445783 3.461078 3.458333 3.455621 3.447059
[171] 3.444444 3.453488 3.439306 3.448276 3.445714
[176] 3.460227 3.474576 3.488764 3.474860 3.461111
[181] 3.464088 3.456044 3.464481 3.478261 3.481081
[186] 3.483871 3.475936 3.489362 3.486772 3.489474
[191] 3.486911 3.479167 3.466321 3.458763 3.461538
[196] 3.454082 3.446701 3.454545 3.462312 3.465000

標本平均 $y = (y_n)_{n=1}^{200}$ の推移を調べるために、横軸に $n$ 、縦軸に $y_n$ を取ると次のグラフが得られます。

> plot(y, type = "l")

f:id:wakabame:20180421154917p:plain

データの平均から、さいころの平均は推定できるかな？

上のグラフを見ると、標本平均 $y_n$ はデータの量が増えるにつれて母平均（真の平均）である $\frac72$ に近づいているような気がします。
それでは、さいころの話題から離れて、母平均が $\mu$ であるような確率変数 $X$ *2があるとします。自然数 $n$ に対して、 $n$ 個の独立な確率変数の列 $X_1, X_2 ,\cdots ,X_n$ を $X$ と同じ分布を持つものとし、それらによる平均を
\begin{equation}
\overline{X} = \dfrac{X_1+X_2+\cdots+X_n}{n}
\end{equation}
とします。さいころの例を用いると、この実現値が $y_n$ に対応します。
$\mu$ を推定するのにあたって、 $n$ を大きくすれば $\overline{X}$ が $\mu$ に近づくことは重要な性質です。しかしながら、 $\overline{X}$ もまた確率的に変動しうるので、 $\overline{X}$ が $\mu$ に近い確率は $n$ を大きくすると1に近づくか？という風に言い換えなくてはなりません。これを数式に直すと、

任意の $\varepsilon >0$ に対して、 $n\rightarrow \infty$ のとき $P\left(|\overline{X}-\mu |< \varepsilon \right) \rightarrow 1$ が成立すること

と表せます。逆に、これを満たすとき、 $\overline{X}$ は $\mu$ の一致推定量である、と言います。

標本平均は母平均の一致推定量になっていることを示そう！

これは確率変数が連続型であれ、離散型であれ、緩い条件のもとで示すことができます。
まずは連続型であるときを考えます。 $X$ の確率分布が連続型であるとき、確率密度を表す非負の関数 $f$ が存在して、
$\begin{align*} &P( a < X \le b ) = \int_a^b f(x)\,dx \\ &\mu = E[X] = \int_\mathbb{R} x f(x)\,dx\\ &\sigma^2 =V[X] =\int_\mathbb{R} (x-\mu)^2f(x)\,dx \end{align*}$
が成立するのでした。 $\overline{X}$ の確率密度は $f(x_1)f(x_2)\cdots f(x_n)$ であることから、平均と分散の加法性、
$\begin{align*} &\mu_n := E[\overline{X}] =\frac{1}{n}\int_{\mathbb{R}^n}(x_1+x_2+\cdots + x_n) f(x_1)f(x_2)\cdots f(x_n)\, dx_1dx_2\cdots dx_n =\mu \\ &\sigma_n^2 :=V[\overline{X}] =\dfrac{\sigma^2}{n} \end{align*}$
が成立します。 $\bar{x} =\frac{x_1+x_2+\cdots + x_n}{n}$ 、 $x =(x_1,x_2,\cdots x_n)$ 、 $F(x)=f(x_1)f(x_2)\cdots f(x_n)$ と書くことにすると最後の等式から、
$\begin{align*} 1&= \frac{1}{ \sigma_n^2} \int_{\mathbb{R}^n} (\bar{x}-\mu)^2F(x)\,dx\\ &\ge \frac{n^2}{ \sigma^2} \int_{\{ |\bar{x}-\mu| \ge \varepsilon\}} (\bar{x}-\mu)^2 F(x)\,dx\\ &\ge \frac{n^2 \varepsilon^2}{\sigma^2} \int_{\{|\bar{x}-\mu| \ge \varepsilon\} } F(x)\,dx\\ &=\frac{n^2 \varepsilon^2}{\sigma^2} P(|\bar{X}-\mu| \ge \varepsilon) \end{align*}$
が成立します。ここで１行目から２行目への不等式は積分範囲を制限していることから、２行目から３行目への不等式は被積分部分において $(\bar{x}-\mu)^2 \ge \varepsilon^2$ であることから従います。両辺を $\frac{n^2 \varepsilon^2}{\sigma^2}$ で割ると、
$P(|\bar{X}-\mu| \ge \varepsilon) \le \frac{\sigma^2}{\varepsilon^2n^2}\rightarrow 0,\quad$ as $n\rightarrow \infty$
が得られ、ここから一致性が従います。
また、ここでの証明が為されるには母平均 $\mu$ と母分散 $\sigma^2$ はそれぞれ定義されなくてはならないことを念のために注意しておきます。つまりコーシー分布のような裾の広い分布では平均や分散を考えられないので、あらかじめ確率分布には可積分性についてのある程度の仮定が必要になります。

続いて離散型の場合を考えます。 $\bar{X}$ の実現値が $M$ 個あるとし、その値を $(\bar{x}_i)_{i=1}^M$ とします。それぞれの確率を $P(\bar{X} = \bar{x_i} )= p_i$ とし、添字集合 $I$ を $I=\{i = 1,2,\cdots ,M| \, |\bar{x}_i-\mu| \ge \varepsilon\}$ と定めれば、先ほどと同様の手順で一致性を示すことができます。