モンテカルロ積分のはなし
モンテカルロ積分って？
モンテカルロ積分とは解析的に解くのが難しい積分を、数値計算によって近似値を求める手法の一つです。
まず、ある積分
を考えます。
f(x) を何らかの確率密度関数 p(x) を使ってg(x)p(x) の形に書き換えられるとき、
\int f(x) dx = \int g(x)p(x)dx = \mathbb{E}_p[g(x)], \quad X\sim p
と期待値の形に変換できます。
あとは、分布 p(x) から乱数 X_1, X_2, \dots, X_n をサンプリングし、g(X_i)の標本平均をとると、
\hat{I}_n = \frac{1}{n}\sum_{i=1}^n g(X_i)
大数の法則より、n が大きいとき
\hat{I}_n \rightarrow \mathbb{E}_p[g(X)]
となり、目的の積分の近似値が得られます。
具体的な計算例
積分
I = \int_0^{\infty} \sin{(x)}e^{-x^2}dx
について、一様乱数 U\sim \text{Uniform}(0, 1) を用いてモンテカルロ積分で近似値を求めることを考えます。
近似1: 一様乱数
X = U/(1-U) と変換する。ヤコビアンは
\frac{dX}{dU} = \frac{1}{(1-U)^2}
※ dX/dU\gt0より、乱数の範囲(0, 1) と積分範囲 (0, \infty) が一対一で対応
積分値は、
\begin{align*}
I &= \int_0^{\infty} \sin{(x)}e^{-x^2}dx\\
&= \int_0^1 \sin{\left(\frac{u}{1-u}\right)}e^{-\left(\frac{u}{1-u}\right)^2}\cdot\frac{1}{(1-u)^2}du\\
&= \mathbb{E}_{U}\left[\sin{\left(\frac{u}{1-u}\right)}e^{-\left(\frac{u}{1-u}\right)^2}\cdot\frac{1}{(1-u)^2}\right],\quad U\sim \text{Uniform}(0, 1)
\end{align*}
よって、
g_1(u) = \sin{\left(\frac{u}{1-u}\right)}e^{-\left(\frac{u}{1-u}\right)^2}\cdot\frac{1}{(1-u)^2}
と置くと、
\hat{I}_n = \frac{1}{n}\sum_{i=1}^n g_1(U_i)
で、I の近似値を計算できます。
近似2: 指数分布に変換
モンテカルロ積分では、積分値を様々な分布の期待値の形式に書き換えることができます。
今、Y = X^2 と置くと、X = \sqrt{Y}, dx/dy = 1 / 2\sqrt{Y} です。
よって、
\begin{align*}
I &= \int_0^{\infty} \sin{(x)}e^{-x^2}dx\\
&= \int_0^{\infty}\sin{(\sqrt{y})}e^{-y}\cdot\frac{1}{2\sqrt{y}}dy\\
&= \int_0^{\infty}\frac{\sin{(\sqrt{y})}}{2\sqrt{y}}\cdot e^{-y}dy\\
&= \mathbb{E}_Y\left[\frac{\sin{(\sqrt{y})}}{2\sqrt{y}}\right], \quad Y\sim \text{Exp}(1)
\end{align*}
これは、\text{Exp}(1) に関する期待値の形になっています。
一様分布 U\sim \text{Uniform}(0, 1) に従う乱数について、Y = -\log{U} と変換すると、
指数分布 \text{Exp}(1) に従う乱数が得られることから、
\begin{align*}
\hat{I}_n &= \frac{1}{n}\sum_{i=1}^ng_2(U_i)
\end{align*}
ただし、
g_2(U) = \frac{\sin{(\sqrt{-\log{U}})}}{2\sqrt{-\log{U}}}
近似1, 2で一様乱数から I を推定するための関数g_1, g_2 が得られました。
しかし、実は u\rightarrow 1 で g_1 は (1-u)^{-2}が発散してしまい、
g_2 は理論的に有限値に収束しますが、-\log{U} が極端に小さいときに丸め誤差が大きくなるという欠点があります。
そこで数値的に安定な別の変換を考えます。
近似3: レイリー分布に変換
元の積分を部分積分を使って変形します。
\begin{align*}
I &= \int_0^{\infty} \sin{(x)}e^{-x^2}dx\\
&= \left[-\cos{(x)}e^{-x^2}\right]_0^{\infty} - \int_0^{\infty}2x\cos{(x)}e^{-x^2}dx\\
&= 1 - \int_0^{\infty}\cos{(x)}\cdot 2xe^{-x^2}dx
\end{align*}
ここで、第二項に出てきた 2xe^{-x^2} は、レイリー分布という確率分布の密度関数になっています。
レイリー分布は一つのパラメータ\sigma をもち、確率密度関数は
\text{Rayleigh}(\sigma) = \frac{x}{\sigma^2}\exp{-\frac{x^2}{2\sigma^2}}
したがって、
\begin{align*}
I = 1 - \mathbb{E}[\cos{W}], \quad W\sim \text{Rayleigh}(1/\sqrt{2})
\end{align*}
また、レイリー分布は一様乱数 U \sim \text{Uniform}(0, 1) を使って
で生成できることから、
g_3(u) = \cos{\sqrt{-\log{u}}}
として、
\hat{I}_n = \frac{1}{n}\sum_{i=1}^ng_3(U_i)
で推定できます。
\cos{x} は常に[-1, 1] の範囲に収まるため、数値的に安定になります。
数値シミュレーション
最後に実際に乱数を生成して、I = \int_0^{\infty} \sin{(x)}e^{-x^2}dx を推定したときの結果がどのようになるかを確認してみます。
import numpy as np
from matplotlib import pyplot as plt
import japanize_matplotlib
n = 1000
u = np.random.default_rng(42).uniform(0, 1, n)
eps = 1e-12
def g1(u, eps=eps):
u = np.clip(u, 0.0, 1.0 - eps)
x = u / (1.0 - u)
return np.sin(x) * np.exp(-x * x) / (1.0 - u) ** 2
def g2(u, eps=eps):
y = -np.log( np.clip(u, eps, 1.0) )
s = np.sqrt(y)
val = np.sin(s) / (2.0 * s)
return val
def g3(u, eps=eps):
x = np.sqrt(-np.log(np.clip(u, eps, 1.0)))
return 1.0 - np.cos(x)
def running_mean(u, g):
vals = g(u)
csum = np.cumsum(vals)
n = np.arange(1, len(vals)+1)
return csum / n
rm1 = running_mean(u, g1)
rm2 = running_mean(u, g2)
rm3 = running_mean(u, g3)
plt.figure(figsize=(7,4))
plt.plot(rm1, label="g1 Direct")
plt.plot(rm2, label="g2 Exponential")
plt.plot(rm3, label="g3 Rayleigh")
plt.axhline(I_true, linestyle="--")
plt.xlabel("n")
plt.ylabel("推定値")
plt.legend()
plt.show()
いずれも、解析値に近い値に収束しています。
しかし、一様分布を直接使った g_1 は n=1000 でもかなり値が暴れているし、指数分布変換を行った g_2 も まだ収束しきっていないように見えます。
一方で、レイリー分布に変換した g_3 はかなり値が安定しており、理論値への収束も早いことが確認できます。
おわりに
モンテカルロ積分について、様々な求め方と数値シミュレーションを実施しました。
間違いなどありましたら、コメントでご指摘いただけると幸いです。
参考
Discussion