In [4]:

    

%pylab inline

# Box-Mullerアルゴリズムを用いて標準正規分布N(0,1)に従う乱数を生成
import numpy as np
from scipy.stats import uniform, norm
import matplotlib.pyplot as plt

# 独立した一様分布からそれぞれ一様乱数を生成
np.random.seed()
N = 100000
rv1 = uniform(loc=0.0, scale=1.0)
rv2 = uniform(loc=0.0, scale=1.0)
U1 = rv1.rvs(N)
U2 = rv2.rvs(N)

# Box-Mullerアルゴリズムで正規分布に従う乱数に変換
# 2つの一様分布から2つの標準正規分布が得られる
X1 = np.sqrt(-2 * np.log(U1)) * np.cos(2 * np.pi * U2)
X2 = np.sqrt(-2 * np.log(U1)) * np.sin(2 * np.pi * U2)

# 変換した乱数の分布と標準正規分布の真のpdfを描画
rv = norm(loc=0, scale=1)

# X1の分布
plt.subplot(2, 1, 1)
nbins = 50
plt.hist(X1, nbins, normed=True)
x = np.linspace(rv.ppf(0.01), rv.ppf(0.99), 1000)
y = rv.pdf(x)
plt.plot(x, y, 'r-', lw=2)
plt.xlim((rv.ppf(0.01), rv.ppf(0.99)))

# X2の分布
plt.subplot(2, 1, 2)
plt.hist(X2, nbins, normed=True)
x = np.linspace(rv.ppf(0.01), rv.ppf(0.99), 1000)
y = rv.pdf(x)
plt.plot(x, y, 'r-', lw=2)
plt.xlim((rv.ppf(0.01), rv.ppf(0.99)))

plt.show()


    

    




Populating the interactive namespace from numpy and matplotlib






    

In [1]:

    

%pylab inline

# 中心極限定理（CLT）による正規乱数生成器
import numpy as np
from scipy.stats import uniform, norm
import matplotlib.pyplot as plt

# 独立した一様分布からそれぞれ一様乱数を生成
np.random.seed()
N = 100000

# 12個の一様乱数を使用
K = 12

# (-0.5, 0.5)の一様乱数を生成
rv = uniform(loc=-0.5, scale=1.0)
U = rv.rvs(K * N).reshape((K, -1))

# K個の一様分布の和から正規分布を生成
Z = np.sum(U, axis=0)

# 変換した乱数の分布と標準正規分布の真のpdfを描画
plt.figure(1)
nbins = 50
plt.hist(Z, nbins, normed=True)

rv = norm(loc=0, scale=1)
x = np.linspace(rv.ppf(0.01), rv.ppf(0.99), 1000)
y = rv.pdf(x)
plt.plot(x, y, 'r-', lw=2)
plt.xlim((rv.ppf(0.01), rv.ppf(0.99)))

plt.show()


    

    




Populating the interactive namespace from numpy and matplotlib






    

In [37]:

    

# Box-MullerアルゴリズムとCLTの比較
import numpy as np
from scipy.stats import uniform, norm
import matplotlib.pyplot as plt

np.random.seed()
N = 100000

# Box-Mullerアルゴリズムによる正規乱数
rv1 = uniform(loc=0.0, scale=1.0)
rv2 = uniform(loc=0.0, scale=1.0)
U1 = rv1.rvs(N)
U2 = rv2.rvs(N)
X1 = np.sqrt(-2 * np.log(U1)) * np.cos(2 * np.pi * U2)
X2 = np.sqrt(-2 * np.log(U1)) * np.sin(2 * np.pi * U2)

# CLTによる正規乱数
K = 12
rv = uniform(loc=-0.5, scale=1.0)
U = rv.rvs(K * N).reshape((K, -1))
Z = np.sum(U, axis=0)

# 重ねて描画
plt.figure(1)
nbins = 50
plt.hist(X1, nbins, normed=True)
plt.hist(Z, nbins, normed=True)

rv = norm(loc=0, scale=1)
x = np.linspace(-4, 4, 1000)
y = rv.pdf(x)
plt.plot(x, y, 'r-', lw=2)
plt.xlim((-4, 4))


    

    
Out[37]:





(-4, 4)