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

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

かなり大雑把に言えば，中心極限定理とは

「同じ分布に従う試行を何度も繰り返すと，トータルで見れば正規分布っぽい分布に近付く」

という定理です．

もう少し数学の言葉を用いて説明するならば，「独立同分布の確率変数列$\{X_n\}$の和$\sum_{k=1}^{n}X_k$は，$n$が十分大きければ正規分布に従う確率変数に近い」という定理です．

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

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

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

準備

まずは準備として

ベルヌーイ分布
二項分布

を復習します．

ベルヌーイ分布

最初に説明するベルヌーイ分布は「コイン投げの表と裏」のような，２つの事象が一定の確率で起こるような試行に関する確率分布です．

いびつなコインを考えて，このコインを投げたときに表が出る確率を$p$とし，このコインを投げて

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

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

$X(\text{表})=1$
$X(\text{裏})=0$

と表すことにしましょう．

雑な言い方ですが，このゲーム$X$はベルヌーイ分布$B(1,p)$に従うといい，$X\sim B(1,p)$と表します．

このように確率的に事象が変化する事柄(いまの場合はコイン投げ)に対して，結果に応じて値(いまの場合は$1$点と$0$点)を返す関数を確率変数といいますね．

つまり，上のゲーム$X$は「ベルヌーイ分布に従う確率変数」ということができます．

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

$\Omega=\{0,1\}$，$\mathcal{F}=2^{\Omega}$($\Omega$の冪集合)とし，関数$\mathbb{P}:\mathcal{F}\to[0,1]$を

$\begin{align*} \mathbb{P}(\{1\})=p,\quad \mathbb{P}(\{0\})=1-p \end{align*}$

で定めると，$(\Omega,\mathcal{F},\mathbb{P})$は確率空間となる．

また，$S=\{0,1\}$，$\mathcal{S}=2^{S}$とすると$(S,\mathcal{S})$は可測空間で，写像$X:\Omega\to S$を

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

で定めると，$X$は$(\Omega,\mathcal{F})$から$(S,\mathcal{S})$への可測写像となる．

このとき，$X$はベルヌーイ分布 (Bernulli distribution)に従うといい，$X\sim B(1,p)$と表す．

このベルヌーイ分布の定義をゲーム$X$に当てはめると

$1\in\Omega$が「表」
$0\in\Omega$が「裏」

に相当し，

$1\in S$が$1$点
$0\in S$が$0$点

に相当します．

$\Omega$と$S$は同じく$0$と$1$からなる集合ですが，意味が違うので注意して下さい．

二項分布

先程のベルヌーイ分布で考えたゲーム$X$を$n$回行うことを考え，このゲームを「ゲーム$Y$」としましょう．

つまり，コインを$n$回投げて，表が出た回数を得点とするのがゲーム$Y$ですね．

ゲーム$X$を繰り返し行うので，何回目に行われたゲームなのかを区別するために，$k$回目に行われたゲーム$X$を$X_k$と表すことにしましょう．

このゲーム$Y$は$X_1,X_2,\dots,X_n$の得点を足し合わせていくので

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

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

このとき，ゲーム$Y$もやはり確率変数で，このゲーム$Y$は二項分布$B(n,p)$に従うといい，$Y\sim B(n,p)$と表します．

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

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

$\Omega’=\Omega^n$，$\mathcal{F}’=2^{\Omega}$とし，測度$\mathbb{P}’:\mathcal{F}\to[0,1]$を

$\begin{align*} \mathbb{P}'(\{(k_1,k_2,\dots,k_n)\}) =&\prod_{i=1}^{n}\mathbb{P}(\{k_i\}) \\\biggl(=&\mathbb{P}(\{k_1\})\times\mathbb{P}(\{k_2\})\times\dots\times\mathbb{P}(\{k_n\})\biggr) \end{align*}$

で定めると，$(\Omega’,\mathcal{F}’,\mathbb{P}’)$は確率空間となる．

また，$S=\{0,1,\dots,n\}$，$\mathcal{S}=2^{S}$とすると$(S,\mathcal{S})$は可測空間で，写像$Y:\Omega\to S$を

$\begin{align*} Y((k_1,k_2,\dots,k_n))=\sum_{i=1}^{n}k_i \end{align*}$

で定めると，$Y$は$(\Omega’,\mathcal{F}’)$から$(S,\mathcal{S})$への可測写像となる．

このとき，$Y$は二項分布 (binomial distribution)に従うといい，$Y\sim B(n,p)$と表す．

$k=k_1+k_2+\dots+k_n$ ($k_i\in\Omega$)なら，$\mathbb{P}(\{(k_1,k_2,\dots,k_n)\})$は$n$回コインを投げて$k$回表が出る確率がなので，反復試行の考え方から

$\begin{align*} \mathbb{P}((k_1,k_2,\dots,k_n)) =&\binom{n}{k}p^k(1-p)^{n-k} \\=&\frac{n!}{k!(n-k)!}p^k(1-p)^{n-k} \end{align*}$

となりますね．

この二項分布の定義をゲーム$Y$に当てはめると

$0\in\Omega$が「表が$1$回も出ない」
$1\in\Omega$が「表がちょうど$1$回出る」
$2\in\Omega$が「表がちょうど$2$回出る」
……
$n\in\Omega$が「表がちょうど$n$回出る」

に相当し，

$0\in S$が$0$点
$1\in S$が$1$点
$2\in S$が$2$点
……
$n\in S$が$n$点

に相当します．

中心極限定理

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

二項分布のシミュレート

ここでは$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

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

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

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

using Plots;gr()

#コイン投げの表の確率

p = 0.3

#コイン投げをする回数

n = 10

#試行回数

N = 100

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

face = zeros(N)

#N回中I回目を考える

for I = 1:N

#n個の乱数を生成する(n個のサイコロをふる)

coin = rand(n)

for i = 1:n

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

if coin[i] < p

face[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 2	#ヒストグラム 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

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

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

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

using Plots;gr()

#コイン投げの表の確率

p = 0.3

#コイン投げをする回数

n = 30

#試行回数

N = 10000

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

face = zeros(N)

#N回中I回目を考える

for I = 1:N

#n個の乱数を生成する(n個のサイコロをふる)

coin = rand(n)

for i = 1:n

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

if coin[i] < p

face[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*}$

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

この式を変形すると

$\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 = <span class="w-impact">1000000</span>

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

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

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

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

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

using Plots;gr()

using LaTeXStrings

#コイン投げの表の確率

p = 0.3

#コイン投げをする回数

n = 100

#試行回数

N = <span class="w-impact">1000000</span>

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

face = zeros(N)

#N回中I回目を考える

for I = 1:N

#n個の乱数を生成する(n個のサイコロをふる)

coin = rand(n)

for i = 1:n

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

if coin[i] < p

face[I] += 1

end

#n回表が出る確率を全て0にセットする

distribution = zeros(n)

#i回表が出る場合を考える

for i = 1:n

for I = 1:N

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

if face[I] == i

distribution[i] += 1/N

end

#B(100,0.3)の分布関数

plot(distribution, label = "", color = :red, linewidth=2)

一般に，ベルヌーイ分布$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のコードを表示

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

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

#平均30，分散21の正規分布の確率密度関数

mu = n*p

Sigma = n*p*(1-p)

f(x) = exp(-(x-mu)^2/(2*Sigma))/sqrt(2*pi*Sigma)

#中心極限定理

x = range(0,100; length = 2001)

plot!(x,f.(x); label = L"\mathcal{N}(30,21)", color = :blue, linewidth=1, style = :dash)

となり，確かに近いことが見てとれますね！

確かにシミュレーションから中心極限定理が成り立っていそうなことが分かりましたね．