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 と同じ分布を持つものとし、それらによる平均を
 \overline{X} = \dfrac{X_1+X_2+\cdots+X_n}{n}
とします。さいころの例を用いると、この実現値が  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}
が成立します。ここで1行目から2行目への不等式は積分範囲を制限していることから、2行目から3行目への不等式は被積分部分において   (\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\} と定めれば、先ほどと同様の手順で一致性を示すことができます。

また、証明の中で積分範囲を制限して作った不等式をチェビシェフの不等式と呼び、今回示した標本平均の一致性を大数の法則と呼びます。

参考文献

この記事での記法はなるべく次の参考文献のそれに準拠させました:

読み方が間違っているなどのご指摘があればコメント下されば幸いです。
次回は不偏性について考えていきます。

*1:実は個人的に計算でつまづいたところがあったので、自分用メモの側面もあります。

*2:ある確率分布に従って数字がランダムに出てくるもの