In [27]:
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate
from scipy.stats import norm
# 単純な重点サンプリングの例
a = 5.0
# 被積分関数
f = norm.pdf
h = lambda x: x > a
y = lambda x: h(x) * f(x)
# scipy.integrateでの積分
I1 = scipy.integrate.quad(f, a, np.inf)[0]
I2 = scipy.integrate.quad(y, -np.inf, np.inf)[0]
print "scipy.integrate:", I1, I2
N = 1000
# 普通のモンテカルロ積分の場合
# サンプルxを青色の標準正規分布N(x|0,1)から生成しているため
# ほとんどのサンプルが5より小さい範囲からしか生成されない
# つまり、h(x)が0となってしまう
x = norm.rvs(size=N)
I = np.mean(h(x))
print "normal monte carlo integration:", I
# 重点サンプリングの場合
# 重点関数g(x)として平均が5にくる正規分布N(x|5,1)を使う
g = norm(loc=5, scale=1).pdf
x = norm(loc=5, scale=1).rvs(size=N)
I = np.mean(f(x) / g(x) * h(x))
print "importance sampling:", I
# グラフ描画
plt.subplot(211)
ix = np.arange(-5, 15, 0.01)
plt.plot(ix, f(ix), label="f(x)")
plt.plot(ix, g(ix), label="g(x)")
plt.plot(ix, h(ix), label="h(x)")
plt.xlim((-5, 15))
plt.ylim((0, 2))
plt.legend(loc="best")
# 被積分関数yの値がある部分をズームインして表示
plt.subplot(212)
plt.plot(ix, y(ix), label="h(x)f(x)")
plt.xlim((4.9, 7))
plt.legend(loc="best")
plt.show()
In [28]:
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate
from scipy.stats import norm
# 単純な重点サンプリングの例
a = 5.0
norms = [(0,1), (5,1), (10,1), (0,3), (5,3), (10,3)]
# 被積分関数
f = norm.pdf
h = lambda x: x > a
y = lambda x: h(x) * f(x)
# 被積分関数を描画
ix = np.arange(-5, 15, 0.01)
plt.plot(ix, f(ix), label="f(x)")
plt.plot(ix, h(ix), label="h(x)")
for m, std in norms:
g = norm(loc=m, scale=std).pdf
plt.plot(ix, g(ix), label="mean=%d std=%d" % (m, std))
plt.xlim((-5, 15))
plt.ylim((0, 2))
plt.legend(loc="best")
plt.show()
# scipy.integrateでの積分
I = scipy.integrate.quad(y, -np.inf, np.inf)[0]
print "scipy.integrate:", I
N = 2000
# 重点サンプリングの場合
# さまざまな平均・標準偏差を持つ正規分布を重点関数gとみなして
# 収束速度がどのように変わるか調査
base_loc = 231
plt.figure(figsize=(10, 6))
for m, std in norms:
plt.subplot(base_loc)
g = norm(loc=m, scale=std).pdf
x = norm(loc=m, scale=std).rvs(size=N)
x = f(x) / g(x) * h(x)
estint = np.cumsum(x) / np.arange(1, N + 1)
esterr = np.sqrt(np.cumsum((x - estint) ** 2)) / np.arange(1, N + 1)
plt.plot(estint, color='red', linewidth=2)
plt.plot(estint + 2 * esterr, color='gray')
plt.plot(estint - 2 * esterr, color='gray')
plt.axhline(y=I)
plt.ylim((0, 1e-06))
plt.title("mean:%d std:%d" % (m, std))
base_loc += 1
plt.tight_layout()
plt.show()
In [45]:
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate
from scipy.stats import uniform, norm
# 被積分関数
f = lambda x: norm(loc=0, scale=1).pdf(x)
h = lambda x: np.exp(-(x-3)**2/2) + np.exp(-(x-6)**2/2)
y = lambda x: h(x) * f(x)
# scipy.integrateでの積分
I = scipy.integrate.quad(y, -np.inf, np.inf)[0]
print "scipy.integrate:", I
N = 2000
# 通常のモンテカルロ積分
x1 = norm.rvs(size=N)
I = np.mean(h(x1))
print "normal monte carlo integrationi:", I
# 重点サンプリングによるモンテカルロ積分
# 被積分領域の範囲を覆うような一様分布を重点関数とする
# サンプリング数が少ない場合は重点サンプリングの方が精度が高い
# （注）テキストでは U(-8, -1) を使えと書いてあるけど意味がわからない
# 被積分関数の範囲からサンプル生成できないのでは？
a, b = 0, 3
u = uniform(loc=a, scale=b-a)
g = lambda x: u.pdf(x)
x2 = u.rvs(size=N)
I = np.mean(f(x2) / g(x2) * h(x2))
print "importance sampling:", I
# グラフ描画
ix = np.arange(-10, 10, 0.01)
plt.plot(ix, f(ix), label="f(x)")
plt.plot(ix, h(ix), label="h(x)")
plt.plot(ix, y(ix), label="h(x) f(x)")
plt.plot(ix, g(ix), label="g(x)")
plt.legend(loc="best")
plt.show()
# 2つの手法の収束の様子を比較
plt.subplot(211)
x1 = h(x1)
estint = np.cumsum(x1) / np.arange(1, N + 1)
esterr = np.sqrt(np.cumsum((x1 - estint) ** 2)) / np.arange(1, N + 1)
plt.plot(estint, color='red', linewidth=2)
plt.plot(estint + 2 * esterr, color='gray')
plt.plot(estint - 2 * esterr, color='gray')
plt.title("convergence (normal monte carlo integration)")
plt.ylim((0.0, 0.2))
plt.subplot(212)
x2 = f(x2) / g(x2) * h(x2)
estint = np.cumsum(x2) / np.arange(1, N + 1)
esterr = np.sqrt(np.cumsum((x2 - estint) ** 2)) / np.arange(1, N + 1)
plt.plot(estint, color='red', linewidth=2)
plt.plot(estint + 2 * esterr, color='gray')
plt.plot(estint - 2 * esterr, color='gray')
plt.title("convergence (importance sampling)")
plt.ylim((0.0, 0.2))
plt.tight_layout()
plt.show()
In [14]:
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate
from scipy.stats import norm, expon
# Rによるモンテカルロ法入門の例3.5
a = 4.5
# サンプリング数
N = 500
# 被積分関数
f = norm.pdf
h = lambda x: x > a
y = lambda x: h(x) * f(x)
# scipy.integrateでの積分
I = scipy.integrate.quad(y, -np.inf, np.inf)[0]
print "scipy.integrate:", I
# 通常のモンテカルロ積分の場合
x = norm.rvs(size=N)
I = np.mean(h(x))
print "normal monte carlo integration:", I
# (1) 重点関数として平均が4.5の正規分布を使用した重点サンプリング
g1 = norm(loc=a, scale=1).pdf
x1 = norm(loc=a, scale=1).rvs(size=N)
I = np.mean(f(x1) / g1(x1) * h(x1))
print "importance sampling (norm):", I
# (2) 重点関数として4.5で切り詰められた指数分布を使用した重点サンプリング
def g2(xlist, a):
"""aで切り詰められた指数分布
TODO: Vectorizationを使った効率的な実装は？"""
result = []
for x in xlist:
if x < a: result.append(0)
else: result.append(np.exp(-(x - a)))
return result
# g2は独自定義の関数であるため直接サンプリングできない
# g2はグラフ描画用として使う
# サンプリングはloc=0の指数分布のサンプルにaを加える（右にずらす）ことで代替する
x2 = expon(loc=0, scale=1).rvs(size=N) + a
I = np.mean(f(x2) / g2(x2, a) * h(x2))
print "importance sampling (truncated expon):", I
# グラフ描画
ix = np.arange(-10, 10, 0.01)
plt.plot(ix, f(ix), label="f(x)")
plt.plot(ix, g1(ix), label="g1(x)")
plt.plot(ix, g2(ix, a), label="g2(x)")
plt.plot(ix, expon(loc=a).pdf(ix), label="expon") # g2と同じ分布になる
plt.legend(loc="best")
plt.show()
# 収束性の評価
plt.subplot(211)
x1 = f(x1) / g1(x1) * h(x1)
estint = np.cumsum(x1) / np.arange(1, N + 1)
esterr = np.sqrt(np.cumsum((x1 - estint) ** 2)) / np.arange(1, N + 1)
plt.plot(estint, color='red', linewidth=2)
plt.plot(estint + 2 * esterr, color='gray')
plt.plot(estint - 2 * esterr, color='gray')
plt.title("convergence (g = norm)")
plt.ylim([0, 8e-6])
plt.subplot(212)
x2 = f(x2) / g2(x2, a) * h(x2)
estint = np.cumsum(x2) / np.arange(1, N + 1)
esterr = np.sqrt(np.cumsum((x2 - estint) ** 2)) / np.arange(1, N + 1)
plt.plot(estint, color='red', linewidth=2)
plt.plot(estint + 2 * esterr, color='gray')
plt.plot(estint - 2 * esterr, color='gray')
plt.title("convergence (g = truncated expon)")
plt.ylim([0, 8e-6])
plt.tight_layout()
plt.show()
In [9]:
def g(x):
if x < 4.5: return 0
return np.exp(-(x - 4.5))
ix = np.arange(-5, 15, 0.01)
iy = [g(x) for x in ix]
plt.plot(ix, iy, label="g(x)")
plt.legend(loc="best")
plt.show()
print scipy.integrate.quad(g, -np.inf, np.inf)[0]
In [8]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import t
from scipy.special import gammaln
# 図3.5の再現
# x = 0.6 としたときのベータ分布のパラメータ上の事後分布
def f(alpha, beta, lam=1, x=0.6, x0=0.5, y0=0.5):
# 対数で計算し、最後にexp()で元に戻す
temp1 = (lam+1) * (gammaln(alpha+beta) - gammaln(alpha) - gammaln(beta))
temp2 = alpha * np.log(x * x0)
temp3 = beta * np.log((1-x) * y0)
return np.exp(temp1 + temp2 + temp3)
a = np.arange(0, 150)
b = np.arange(0, 100)
A, B = np.meshgrid(a, b)
Z = f(A, B)
plt.subplot(1, 2, 1)
plt.contour(A, B, Z)
plt.xlabel("alpha")
plt.ylabel("beta")
# 多変量t分布に従う乱数を生成する
# scipy.stats.tは一変量のみ
# 例3.6の方法
# t分布に従う乱数を生成
# cholesky()はRのchol()と違って転置した値が出てくる
# http://wiki.scipy.org/NumPy_for_Matlab_Users
N = 2 * 10000
df = 3 # t分布の自由度
mu = np.array([50, 45])
E = np.array([220, 190, 190, 180]).reshape(2, 2)
x = t.rvs(df, size=N).reshape(-1, 2)
print x
y = np.transpose(np.dot(np.linalg.cholesky(E), np.transpose(x)))
# 各行に平均を加える（Rと違ってy+muでは動作しない）
for i in range(y.shape[0]):
y[i] += mu
plt.subplot(1, 2, 2)
plt.contour(A, B, Z)
plt.plot(y[:,0], y[:,1], 'r.')
plt.xlabel("alpha")
plt.ylabel("beta")
plt.xlim((0, 150))
plt.ylim((0, 100))
plt.tight_layout()
plt.show()