中心極限定理を実感する｜二項分布でシミュレートしてみた

確率論の重要な定理として中心極限定理があります．

大雑把に言えば，中心極限定理とは「同じ分布に従う試行を何度も繰り返すと，トータルで見れば正規分布っぽくなる」という定理です．

本記事の目的は「中心極限定理がどういうものか実感しようという」というもので，独立なベルヌーイ分布の確率変数列$\{X_n\}$に対して中心極限定理が成り立つ様子をプログラミングでシミュレーションします．

なお，本記事ではJuliaというプログラミング言語を扱っていますが，本記事の主題は中心極限定理のイメージを理解することなので，Juliaのコードが分からなくても問題ないように話を進めます．

準備
1. ベルヌーイ分布
2. 二項分布
中心極限定理

準備

まずは準備として

ベルヌーイ分布
二項分布

について確認しておきましょう．

ベルヌーイ分布

表が出る確率が$p$のいびつなコインを投げて

表が出れば$1$点
裏が出れば$0$点

という「ゲーム$X$」を考えます．

このゲーム$X$のように確率的に事象が変化する事柄の結果（事象）に応じて，値を返す関数を確率変数といいます．

ベルヌーイ分布の厳密に定義を述べると以下のようになります（分からなければ飛ばしても問題ありません）．

確率空間$(\Omega,\mathcal{F},\mathbb{P})$を

$\Omega=\{0,1\}$
$\mathcal{F}=2^{\Omega}$（$\Omega$の冪集合）
$\mathbb{P}(\{1\})=p$, $\mathbb{P}(\{0\})=1-p$

で定める．また，可測空間$(S,\mathcal{S})$を$S=\{0,1\}$, $\mathcal{S}=2^{S}$とするとき，

$\begin{align*}X(0)=0,\quad X(1)=1\end{align*}$

で定まる可測写像$X:\Omega\to S$はベルヌーイ分布 (Bernulli distribution)に従うといい，$X\sim B(1,p)$と表す．

上のゲーム$X$において，

「表」はベルヌーイ分布の$1\in\Omega$
「裏」はベルヌーイ分布の$0\in\Omega$

に相当し，

$1$点はベルヌーイ分布の$1\in S$
$0$点はベルヌーイ分布の$0\in S$

に相当します．つまり，上のゲーム$X$はベルヌーイ分布に従うということができます．

二項分布

ベルヌーイ分布でも考えた表が出る確率が$p$のいびつなコインを$n$回投げて，表が出た回数を得点とする「ゲーム$Y$」を考えます．

このゲーム$Y$はゲーム$X$を繰り返し行って点数を加算していくゲームと言えるので，ゲーム$Y$の$k$回目のコイン投げをゲーム$X_k$と呼ぶことにすれば

$\begin{align*}Y=X_1+X_2+\dots+X_n\end{align*}$

と表すことができますね．

二項分布の厳密に定義を述べると以下のようになります（こちらも分からなければ飛ばしても問題ありません）．

$(\Omega, \mathcal{F},\mathbb{P})$を上のベルヌーイ分布の定義での確率空間とし，確率空間$(\Omega’,\mathcal{F}’,\mathbb{P}’)$を

$\Omega’=\Omega^n$
$\mathcal{F}’=2^{\Omega’}$
$\mathbb{P}'(\{(x_1,x_2,\dots,x_n)\})=\prod_{k=1}^{n}\mathbb{P}(\{x_k\})$

で定める．また，可測空間$(S,\mathcal{S})$を$S=\{0,1,\dots,n\}$, $\mathcal{S}=2^{S}$とするとき，

$\begin{align*}Y((x_1,x_2,\dots,x_n))=\sum_{k=1}^{n}x_k\end{align*}$

で定まる可測写像$X:\Omega’\to S$は二項分布 (binomial distribution)に従うといい，$Y\sim B(n,p)$と表す．

上のゲーム$Y$において

「表が$1$回も出ない」は二項分布の$0\in\Omega$
「表がちょうど$1$回出る」は二項分布の$1\in\Omega$
「表がちょうど$2$回出る」は二項分布の$2\in\Omega$
……
「表がちょうど$n$回出る」は二項分布の$n\in\Omega$

に相当し，

$0$点は二項分布の$0\in S$
$1$点は二項分布の$1\in S$
$2$点は二項分布の$2\in S$
……
$n$点は二項分布の$n\in S$

に相当します．つまり，上のゲーム$Y$は二項分布$B(n,p)$に従うということができます．

中心極限定理

それでは，中心極限定理のイメージの説明に移りますが，そのために二項分布をシミュレートしていきます．

二項分布のシミュレート

ここでは$p=0.3$の二項分布$B(n,p)$を考えます．

つまり，「表が30%の確率で出る歪んだコインを$n$回投げたときに，合計で何回表が出るか」を考えます．

$n=10$のとき

$n=10$の場合，つまり$B(10,0.3)$を考えましょう．

このとき，「表が$30\%$の確率で出る歪んだコインを$10$回投げたときに，合計で何回表が出るか」を考えることになるわけですが，表が$3$回出ることもあるでしょうし，$1$回しか出ないことも，$7$回出ることもあるでしょう．

しかし，さすがに$10$回投げて$1$回も表が出なかったり，$10$回表が出るということはあまりなさそうに思えますね．

ということで，「表が$30\%$の確率で出る歪んだコインを$10$回投げて，表が出る回数を記録する」という試行を$100$回やってみましょう．

Juliaのコードを表示

using Plots;gr()

p = 0.3  #コイン投げの表の確率
n = 10  #コイン投げをする回数
N = 100  #試行回数
face = zeros(N)  #N回それぞれの試行での表の回数を全て0にセットする

for I = 1:N  #N回中I回目を考える
    coin = rand(n)  #n個の乱数（0〜1）を生成する(n個のサイコロをふる)
    for i = 1:n
        if coin[i] < p
            face[I] += 1  #i個目の乱数が0.3未満なら表とみなして，I回目の表の表の回数に+1
        end
    end
end

plot(face, label = "", color = :red, marker = :auto)  #各試行での表の回数

using Plots;gr()

p = 0.3 #コイン投げの表の確率

n = 10 #コイン投げをする回数

N = 100 #試行回数

face = zeros(N) #N回それぞれの試行での表の回数を全て0にセットする

for I = 1:N #N回中I回目を考える

coin = rand(n) #n個の乱数（0〜1）を生成する(n個のサイコロをふる)

for i = 1:n

if coin[i] < p

face[I] += 1 #i個目の乱数が0.3未満なら表とみなして，I回目の表の表の回数に+1

end

plot(face, label = "", color = :red, marker = :auto) #各試行での表の回数

結果は以下の図になりました．

1回目は表が$1$回も出なかったようで，17回目と63回目と79回目に表が$6$回出ていてこれが最高の回数ですね．

この図を見ると，$3$回表が出ている試行が最も多いように見えますね．

そこで，表が出た回数をヒストグラムに直してみましょう．

Juliaのコードを表示

plot(face, label = "", color = :red, st = :histogram)  #ヒストグラム

1	plot(face, label = "", color = :red, st = :histogram) #ヒストグラム

確かに，$3$回表が出た試行が最も多く$30$回となっていますね．

$n=30$のとき

$n=30$の場合，つまり$B(30,0.3)$を考えましょう．

つまり，「$30$回コインを投げて表の回数を記録する」というのを1回の試行として，この試行を$10000$回行ったときのヒストグラムを出力すると以下のようになりました．

Juliaのコードを表示

using Plots;gr()

p = 0.3  #コイン投げの表の確率
n = 30  #コイン投げをする回数
N = 10000  #試行回数
face = zeros(N)  #N回それぞれの試行での表の回数を全て0にセットする

for I = 1:N  #N回中I回目を考える
    coin = rand(n)  #n個の乱数（0〜1）を生成する(n個のサイコロをふる)
    for i = 1:n
        if coin[i] < p
            face[I] += 1  #i個目の乱数が0.3未満なら表とみなして，I回目の表の表の回数に+1
        end
    end
end

plot(face, label = "", color = :red, st = histogram)  #各試行での表の回数

using Plots;gr()

p = 0.3 #コイン投げの表の確率

n = 30 #コイン投げをする回数

N = 10000 #試行回数

face = zeros(N) #N回それぞれの試行での表の回数を全て0にセットする

for I = 1:N #N回中I回目を考える

coin = rand(n) #n個の乱数（0〜1）を生成する(n個のサイコロをふる)

for i = 1:n

if coin[i] < p

face[I] += 1 #i個目の乱数が0.3未満なら表とみなして，I回目の表の表の回数に+1

end

plot(face, label = "", color = :red, st = histogram) #各試行での表の回数

先ほどより，ガタガタではなく少し滑らかに見えてきました．

そこで，もっと$n$を大きくしてみましょう．

$n=100$のとき

$n=100$の場合，つまり$B(100,0.3)$を考えましょう．

試行回数$1000000$回でシミュレートすると，以下のようになりました(コードは省略)．

とても綺麗な釣鐘型になりましたね！

釣鐘型の確率密度関数として有名なものといえば正規分布ですね．

このように，二項分布$B(n,p)$は$n$を大きくしていくと，正規分布のような雰囲気を醸し出すことが分かりました．

中心極限定理

二項分布$B(n,p)$に従う確率変数$Y$は，ベルヌーイ分布$B(1,p)$に従う独立な確率変数$X_1,\dots,X_n$の和として表せるのでした：$Y=X_1+\dots+X_n$．

この和$Y$が$n$を大きくすると正規分布の確率密度関数のような形状に近付くことは上でシミュレートした通りですが，実は$X_1,\dots,X_n$がベルヌーイ分布でなくても，独立同分布の確率変数$X_1,\dots,X_n$の和でも同じことが起こります．

このような同一の確率変数の和について成り立つ次の定理を中心極限定理といいます．

厳密に書けば以下のようになります．

平均$\mu\in\R$，分散$\sigma^2\in(0,\infty)$の独立同分布に従う確率変数列$X_1,X_2,\dots$に対して

$\begin{align*}Z_n:=\frac{1}{\sqrt{n\sigma^2}}\sum_{k=1}^{n}(X_k-\mu)\end{align*}$

で定まる確率変数列$Z_1,Z_2,\dots$は，標準正規分布に従う確率変数$Z$に法則収束する：

$\begin{align*}Z_n\xrightarrow[]{\mathcal{L}}Z\sim\mathcal{N}(0,1).\end{align*}$

細かい言い回しなどは，この記事ではさほど重要ではありませんので，ここでは「$n$が十分大きければ確率変数

$\begin{align*}Z_n=\frac{1}{\sqrt{n\sigma^2}}\sum_{k=1}^{n}(X_k-\mu)\end{align*}$

は$n$が大きいときだいたい標準正規分布に従う」という程度の理解で問題ありません．

この式を変形すると

$\begin{align*}X_1+\dots+X_n=\sqrt{n\sigma^2}Z_n+n\mu\end{align*}$

となります．

中心極限定理より，$n$が十分大きければ$Z_n$は標準正規分布に従う確率変数$Z$に近いので，確率変数$X_1+\dots+X_n$は確率変数$\sqrt{n\sigma^2}Z+n\mu$に近いと言えますね．

確率変数に数をかけても縮尺が変わるだけですし，数を足しても平行移動するだけなので，結果として$X_1+\dots+X_n$は正規分布と同じ釣鐘型に近くなるわけですね．

シミュレートして実感する

先ほどシミュレートした$n=100$の場合のヒストグラムは$1000000$回のシミュレートなので，ヒストグラムの度数を$1000000$で割ると$B(100,0.3)$の確率関数がシミュレートされますね．

Juliaのコードを表示

using Plots;gr()
using LaTeXStrings

p = 0.3  #コイン投げの表の確率
n = 100  #コイン投げをする回数
N = 1000000  #試行回数
face = zeros(N)  #N回それぞれの試行での表の回数を全て0にセットする

for I = 1:N  #N回中I回目を考える
    coin = rand(n)  #n個の乱数（0〜1）を生成する(n個のサイコロをふる)
    for i = 1:n
        if coin[i] < p
            face[I] += 1  #i個目の乱数が0.3未満なら表とみなして，I回目の表の回数に+1
        end
    end
end

distribution = zeros(n)  #n回表が出る確率を全て0にセットする

for i = 1:n  #i回表が出る場合を考える
    for I = 1:N
        if face[I] == i
            distribution[i] += 1/N  #I回目の試行で表がi回出ていれば，i回の表が出る確率に+1/N
        end
    end
end

plot(distribution, label = "", color = :red, linewidth=2)  #B(100,0.3)の分布関数

using Plots;gr()

using LaTeXStrings

p = 0.3 #コイン投げの表の確率

n = 100 #コイン投げをする回数

N = 1000000 #試行回数

face = zeros(N) #N回それぞれの試行での表の回数を全て0にセットする

for I = 1:N #N回中I回目を考える

coin = rand(n) #n個の乱数（0〜1）を生成する(n個のサイコロをふる)

for i = 1:n

if coin[i] < p

face[I] += 1 #i個目の乱数が0.3未満なら表とみなして，I回目の表の回数に+1

end

distribution = zeros(n) #n回表が出る確率を全て0にセットする

for i = 1:n #i回表が出る場合を考える

for I = 1:N

if face[I] == i

distribution[i] += 1/N #I回目の試行で表がi回出ていれば，i回の表が出る確率に+1/N

end

plot(distribution, label = "", color = :red, linewidth=2) #B(100,0.3)の分布関数

一般に，ベルヌーイ分布$B(1,p)$に従う確率変数$X$は

平均は$p$
分散は$p(1-p)$

であることが知られています．

よって，中心極限定理より，二項分布$B(100,0.3)$に従う確率変数$X_1+\dots+X_{100}$ ($X_1,\dots,X_n\sim B(1,0.3)$は，確率変数

$\begin{align*} &\sqrt{100\times0.3\cdot(1-0.3)}Z+100\cdot0.3\cdot(1-0.3) \\=&\sqrt{21}Z+30 \end{align*}$

に十分近いはずです．この確率変数は

平均は$30$
分散は$21$

の正規分布に従うので，この確率密度関数を上でシミュレートした$B(100,0.3)$の確率関数と重ねて表示させると

Juliaのコードを表示

mu = n*p
Sigma = n*p*(1-p)
f(x) = exp(-(x-mu)^2/(2*Sigma))/sqrt(2*pi*Sigma)  #平均30，分散21の正規分布の確率密度関数

x = range(0,100; length = 2001)
plot!(x,f.(x); label = L"\mathcal{N}(30,21)", color = :blue, linewidth=1, style = :dash)  #中心極限定理