In [8]:
# ベータ乱数を受理・棄却法で生成
# 目標分布（ここではベータ分布）のpdfは既知とする
# 提案分布として一様分布を使用
import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize
from scipy.stats import uniform, beta
np.random.seed()
# 目標分布f
f = beta(a=2.7, b=6.3).pdf
# 提案分布g
# 提案分布から乱数生成するためにgvも保持
gv = uniform
g = gv.pdf
# 分布の上限を指定する定数Mを設定
# ベータ分布のpdfの上限値を指定すればベータ分布をすべて覆える
# 最大値を求めるためにベータ分布のpdfにマイナスをつけて
# 最小値問題に帰着させる
xopt = scipy.optimize.fmin(lambda x: -f(x), 0.0, disp=False)
M = f(xopt)[0]
print "M =", M
# 受理・棄却法
Nsim = 100000
# 提案分布gからの乱数Yを生成
Y = gv.rvs(size=Nsim)
# 一様乱数UをNsim個生成
U = uniform.rvs(size=Nsim)
# Yから受理の条件を満たすサンプルXを残して残りを棄却
X = Y[U <= f(Y) / (M * g(Y))]
print u"サンプル数: %d => %d" % (len(Y), len(X))
print u"実際の受理率 : %f" % (len(X) / float(len(Y)))
print u"理論的な受理率: %f" % (1.0 / M)
# 目標分布を描画
x = np.linspace(0.0, 1.0, 1000)
y = f(x)
plt.plot(x, y, 'r-', lw=2)
# 提案分布（一様分布）を描画
y = M * uniform.pdf(x)
plt.plot(x, y, 'g-', lw=2)
# 受理した乱数の分布を描画
plt.hist(X, bins=50, normed=True)
plt.show()
In [9]:
# 受理されたサンプルと棄却されたサンプルを描く
# 点の数が多すぎるのでNsimを小さくした
# 目標分布f
f = beta(a=2.7, b=6.3).pdf
# 提案分布g
# 提案分布から乱数生成するためにgvも保持
gv = uniform
g = gv.pdf
Nsim = 2000
# 候補密度からの乱数Yを生成
Y = uniform.rvs(size=Nsim)
# 一様乱数UをNsim個生成
U = uniform.rvs(size=Nsim)
# 受理されたサンプルと、棄却されたサンプルのインデックスを計算
acceptedIdx = U <= f(Y) / (M * g(Y))
rejectedIdx = U > f(Y) / (M * g(Y))
# 目標分布を描画
x = np.linspace(0.0, 1.0, 1000)
y = f(x)
plt.plot(x, y, 'r-', lw=2)
# 提案分布を描画
y = M * g(x)
plt.plot(x, y, 'g-', lw=2)
# 受理されたサンプルを描画
plt.scatter(Y[acceptedIdx], U[acceptedIdx] * M * g(Y[acceptedIdx]), color="red")
plt.scatter(Y[rejectedIdx], U[rejectedIdx] * M * g(Y[rejectedIdx]), color="blue")
plt.xlim((0.0, 1.0))
plt.ylim((0.0, 3.0))
plt.show()
In [13]:
# ベータ乱数を受理・棄却法で生成
# 目標分布（ここではベータ分布）のpdfは既知とする
# 提案分布として目標のベータ分布を覆うベータ分布を使用
import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize
from scipy.stats import uniform, beta
np.random.seed()
# 目標分布f
f = beta(a=2.7, b=6.3).pdf
# 提案分布g
gv = beta(a=2.0, b=6.0)
g = gv.pdf
# 分布の上限を指定する定数Mを設定
# f(x)/g(x) <= Mを満たす必要がある
# f(x)/g(x)を最大化するxoptを求める
xopt = scipy.optimize.fmin(lambda x: - f(x) / g(x), 0.0, disp=False)
# そこでの値をMとする
M = f(xopt) / g(xopt)
print "M =", M
# 受理・棄却法
Nsim = 100000
# 提案分布gからの乱数Yを生成
Y = gv.rvs(size=Nsim)
# 一様乱数UをNsim個生成
U = uniform.rvs(size=Nsim)
# Yから受理の条件を満たすサンプルXを残して残りを棄却
X = Y[U <= f(Y) / (M * g(Y))]
print u"サンプル数: %d => %d" % (len(Y), len(X))
print u"実際の受理率 : %f" % (len(X) / float(len(Y)))
print u"理論的な受理率: %f" % (1.0 / M)
# 目標分布を描画
x = np.linspace(0.0, 1.0, 1000)
y = f(x)
plt.plot(x, y, 'r-', lw=2)
# 提案分布（ベータ分布）を描画
y = M * g(x)
plt.plot(x, y, 'g-', lw=2)
# 受理した乱数の分布を描画
plt.hist(X, bins=50, normed=True)
plt.show()
In [11]:
# 受理されたサンプルと棄却されたサンプルを描く
# 点の数が多すぎるのでNsimを小さくした
# 目標分布f
f = beta(a=2.7, b=6.3).pdf
# 提案分布g
gv = beta(a=2.0, b=6.0)
g = gv.pdf
Nsim = 2000
# 候補密度からの乱数Yを生成
Y = uniform.rvs(size=Nsim)
# 一様乱数UをNsim個生成
U = uniform.rvs(size=Nsim)
# 受理されたサンプルと、棄却されたサンプルのインデックスを計算
acceptedIdx = U <= f(Y) / (M * g(Y))
rejectedIdx = U > f(Y) / (M * g(Y))
# 目標分布を描画
x = np.linspace(0.0, 1.0, 1000)
y = f(x)
plt.plot(x, y, 'r-', lw=2)
# 提案分布を描画
y = M * g(x)
plt.plot(x, y, 'g-', lw=2)
# 受理されたサンプルを描画
plt.scatter(Y[acceptedIdx], U[acceptedIdx] * M * g(Y[acceptedIdx]), color="red")
plt.scatter(Y[rejectedIdx], U[rejectedIdx] * M * g(Y[rejectedIdx]), color="blue")
plt.xlim((0.0, 1.0))
plt.ylim((0.0, 5.0))
plt.show()
In [12]:
import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize
from scipy.stats import laplace
np.random.seed()
f = laplace(loc=0, scale=4).pdf
x = np.linspace(-10, 10, 1000)
y = f(x)
plt.plot(x, y, 'r-', lw=2)
Out[12]:
In [54]:
# 正規乱数を受理・棄却法で生成
# 提案分布としてラプラス分布（二重指数分布）を使用
import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize
from scipy.stats import norm, laplace
np.random.seed()
# 目標分布f
f = norm(loc=0.0, scale=1.0).pdf
# 提案分布g
# scale = 1/alphaになる
# 練習問題2.8からもっとも効率がよいのがalpha=1なのでscale=1
gv = laplace(loc=0.0, scale=1.0)
g = gv.pdf
# 分布の上限を指定する定数Mを設定
# f(x)/g(x) <= Mを満たす必要がある
# 練習問題2.8から下式がもっとも効率のよいMとなる
alpha = 1.0
M = np.sqrt(2.0 / np.pi) * (1.0/alpha) * np.exp(alpha**2/2.0)
print "M =", M
# 受理・棄却法
Nsim = 100000
# 提案分布gからの乱数Yを生成
Y = gv.rvs(size=Nsim)
# 一様乱数UをNsim個生成
U = uniform.rvs(size=Nsim)
# Yから受理の条件を満たすサンプルXを残して残りを棄却
X = Y[U <= f(Y) / (M * g(Y))]
print u"サンプル数: %d => %d" % (len(Y), len(X))
print u"実際の受理率 : %f" % (len(X) / float(len(Y)))
print u"理論的な受理率: %f" % (1.0 / M)
# 目標分布を描画
x = np.linspace(-10, 10, 1000)
y = f(x)
plt.plot(x, y, 'r-', lw=2)
# 提案分布を描画
y = M * g(x)
plt.plot(x, y, 'g-', lw=2)
# 受理した乱数の分布を描画
plt.hist(X, bins=50, normed=True)
plt.show()
In [46]:
# 正規乱数を受理・棄却法で生成
# 任意の確率分布にしたがう乱数を生成する
import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize
import scipy.integrate
from scipy.stats import norm
np.random.seed()
# 正規化前の目標分布f
f = lambda x: np.exp(-x**2 / 2) * (np.sin(6*x)**2 + 3 * np.cos(x)**2 * np.sin(4*x)**2 + 1)
# 正規化定数
I = scipy.integrate.quad(f, -Inf, Inf)[0]
# 正規化済み（積分が1）の目標分布
f = lambda x: np.exp(-x**2 / 2) * (np.sin(6*x)**2 + 3 * np.cos(x)**2 * np.sin(4*x)**2 + 1) / I
# 提案分布g
gv = norm(loc=0.0, scale=1.0)
g = gv.pdf
# 分布の上限を指定する定数Mを設定
xopt = scipy.optimize.fmin(lambda x: - f(x) / g(x), 0.0, disp=False)[0]
# そこでの値をMとする
M = f(xopt) / g(xopt)
print "M =", M
# 受理・棄却法
Nsim = 100000
# 提案分布gからの乱数Yを生成
Y = gv.rvs(size=Nsim)
# 一様乱数UをNsim個生成
U = uniform.rvs(size=Nsim)
# Yから受理の条件を満たすサンプルXを残して残りを棄却
X = Y[U <= f(Y) / (M * g(Y))]
print u"サンプル数: %d => %d" % (len(Y), len(X))
print u"実際の受理率 : %f" % (len(X) / float(len(Y)))
print u"理論的な受理率: %f" % (1.0 / M)
# 目標分布を描画
x = np.linspace(-10, 10, 1000)
y = f(x)
plt.plot(x, y, 'r-', lw=2)
# 提案分布を描画
y = M * g(x)
plt.plot(x, y, 'g-', lw=2)
# 受理した乱数の分布を描画
plt.hist(X, bins=50, normed=True)
plt.xlim((-5, 5))
plt.show()
Out[46]: