正規分布のシミュレーション
中心極限定理
同じ分布に従う独立な確率変数
の和は、正規分布に近づくという中心極限定理のシミュレーションを行う。まず、確率変数が連続変数の場合について行い、次に離散変数の場合について行う。
連続な確率変数の場合
確率変数が一様分布に従うものとする。
まず、2つの確率変数の和
の分布について調べる。1,000,000回のシミュレーションを行い、度数分布を描くと図１のようになる。
図１
正規分布のグラフは、赤色の曲線で描かれている。和のヒストグラムが3角形の形で正規分布に沿っていることがわかる。
シミュレーションのスクリプトは、リスト１に示すものである。
リスト１
import numpy
as np
import scipy.stats
as ss
import matplotlib.pyplot
as plt
n_simu = 1_000_000
sample_size = 2
"""
サンプルサイズ２のサンプリングを1,000,000回行う
"""
samples = ss.uniform.rvs(loc=0.0,
scale=1.0, size=(n_simu,sample_size))
samples = (np.sum(samples,
axis=1) - (0.5 * sample_size)) / \
((sample_size *(1/12.0))**0.5)
xcoord = np.linspace(-6,
6, 1000)
ycoord = ss.norm.pdf(xcoord,loc=0,
scale=1)
plt.hist(samples, density=True, bins=np.linspace(-6, 6, 100))
plt.plot(xcoord, ycoord)
plt.title(f'Uniform,
Sample Size = {sample_size}', fontsize=16)
plt.show()
和をとる確率変数の個数を100個に増やしてシミュレーションを行ってみた。図２のヒストグラムを得た。
図２
和
の分布は、滑らかに正規分布（赤色の曲線）に沿っていることがわかる。
図２のグラフは、リスト２のスクリプトで描いた。
リスト２
import numpy
as np
import scipy.stats
as ss
import matplotlib.pyplot
as plt
n_simu = 1_000_000
sample_size = 100
"""
サンプルサイズ100のサンプリングを1,000,000回行う
"""
samples = ss.uniform.rvs(loc=0.0,
scale=1.0, size=(n_simu,sample_size))
samples = (np.sum(samples,
axis=1) - (0.5 * sample_size)) / \
((sample_size
*(1/12.0))**0.5)
xcoord = np.linspace(-6,
6, 1000)
ycoord = ss.norm.pdf(xcoord,loc=0,
scale=1)
plt.hist(samples, density=True, bins=np.linspace(-6, 6, 100))
plt.plot(xcoord, ycoord)
plt.title(f'Uniform,
Sample Size = {sample_size}', fontsize=16)
plt.show()
各確率変数が確率1/2のベルヌーイ分布に従うものとする。
と置くとき、
のサンプリングを1,000,000回行いヒストグラムを描くと、図３となった。
図３
３つの値をとる棒グラフであり、赤色で描かれた正規分布のグラフがなぞるように重なっている。
スクリプトをリスト３に示す。
リスト３
import numpy
as np
import scipy.stats
as ss
import matplotlib.pyplot
as plt
n_simu = 1_000_000
sample_size = 2
"""
サンプルサイズ２のサンプリングを1,000,000回行う
"""
samples = ss.bernoulli.rvs(0.5,
size=(n_simu, sample_size))
samples = (np.sum(samples,
axis=1) - 0.5 * sample_size) / \
((sample_size * (0.5**2))**0.5)
xcoord = np.linspace(-6,
6, 1000)
ycoord = ss.norm.pdf(xcoord,loc=0,
scale=1)
plt.hist(samples, density=True, bins=np.linspace(-6, 6, 12), rwidth=0.1)
plt.plot(xcoord, ycoord)
plt.title(f'Bernoulli, Sample size = {sample_size}',
fontsize=16)
plt.show()
足す乱数の個数を100個とし、
のサンプリングを1,000,000回行うと、図４のようになった。
図４
ヒストグラムは、正規分布によく合うようになっている。
スクリプトをリスト４に示す。
リスト４
import numpy
as np
import scipy.stats
as ss
import matplotlib.pyplot
as plt
n_simu = 1_000_000
sample_size = 100
"""
サンプルサイズ100のサンプリングを1,000,000回行う
"""
samples = ss.bernoulli.rvs(0.5,
size=(n_simu,sample_size))
samples = (np.sum(samples,
axis=1) - 0.5*sample_size) / \
((0.25*sample_size)**0.5)
xcoord = np.linspace(-6,
6, 1000)
ycoord = ss.norm.pdf(xcoord,loc=0,
scale=1)
plt.hist(samples, density=True, bins=np.linspace(-6, 6.0, 60))#,rwidth=1)
plt.plot(xcoord, ycoord)
plt.title(f'Bernoulli, Sample size = {sample_size}',
fontsize=16)
plt.show()
参考文献
Hogg, R. V., McKean, J. W., & Craig, A. T. (2005) Introduction to Mathematical Statistics, sixth edition. Pearson Education, Inc.
Johnson, N. L., Kotz, S., & Balakrishnan, N. (1994) Continuous Univariate Distributions, volume 1, second edition. John Wiley & Sons, Inc.
岡本安晴（2019）いまさら聞けないPythonでデータ分析．丸善出版