書記が数学やるだけ#786 モンテカルロ法

積分の近似法であるモンテカルロ法について。

問題

説明

モンテカルロ法は乱数を用いて積分を近似する方法である。

派生の一つとして重点サンプリング法がある。

解答

最も基本的なモンテカルロ法において，期待値と分散は以下の通りで，サンプルサイズの平方根レートで近似できることがわかる。

重点サンプリング法において，理論上最も効率の良いQは以下の式で表せる。

ではPythonで実装していく。モンテカルロ法の実装自体は至ってシンプルで，以下の複雑な関数の積分においてscipy.integrateと高い一致性を示している。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import uniform, norm
import scipy.integrate
import matplotlib.patches as pat

a, b = 0, 1
f = lambda x: (np.cos(50 * x) + np.sin(20 * x)) ** 2

#被積分関数
x = np.linspace(-0.1, 1.1, 1000)
y = f(x)
plt.plot(x, y)

#作図
ix = np.arange(a, b, 0.001)
iy = f(ix)
verts = [(a, 0)] + list(zip(ix, iy)) + [(b,0)]
poly = pat.Polygon(verts, facecolor='0.8', edgecolor='k')
plt.gca().add_patch(poly)

#scipy.integrateで積分を計算
I = scipy.integrate.quad(f, a, b)[0]
print("scipy.integrate:", I)

#モンテカルロ積分
N = 100000
x = uniform(loc=a, scale=b-a).rvs(size=N)
I = (b - a) * np.mean(f(x))
print("モンテカルロ積分:", I)

#出力
#scipy.integrate: 0.9652009360501453
#モンテカルロ積分: 0.9641977525223131

次に(2)の積分について。標準正規分布の5以上の部分の積分は非常に小さい値であり，基本的モンテカルロ法では0と出力されてしまう。

#被積分関数
f = norm.pdf
h = lambda x: x > 5
y = lambda x: h(x) * f(x)

#scipy.integrateでの積分
I1 = scipy.integrate.quad(f, 5, np.inf)[0]
I2 = scipy.integrate.quad(y, -np.inf, np.inf)[0]
print("scipy.integrate:", I1, I2)

#基本的モンテカルロ積分
N = 1000
x = norm.rvs(size=N)
I = np.mean(h(x))
print("normal monte carlo integration:", I)

#作図
ix = np.arange(-5, 15, 0.01)
plt.plot(ix, f(ix), label="f(x)")
plt.plot(ix, h(ix), label="h(x)")
plt.xlim((-5, 15))
plt.ylim((0, 2))
plt.legend(loc="best")

#出力
#scipy.integrate: 2.8665157035203983e-07 2.866527562360824e-07
#normal monte carlo integration: 0.0

ここでg(x)=N(x|5,1)として重点サンプリング法を用いることで，うまく計算することができるようになる。

#重点サンプリング，g(x)=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)

#グラフ描画
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")

#出力
#importance sampling: 2.660535762457889e-07

本記事のもくじはこちら：