python 「 - [Rによるモンテカルロ法入门」の练习问题3.4

Posted 2021-05-10

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()