前回以下の記事で、積分領域を囲む四角形内にランダムプロットを打っていくことによって積分領域の内外判定比率から積分値を求めることができました。今回はこの例を一般化してみます。
目次
モンテカルロ積分
一様分布に従う確率密度関数
閉区間U[i,j]において、確率変数Xが一様分布に従うとき、その確率密度関数は以下のように定義されます。
$$PDF(x)=\frac{1}{(j-i)}\quad X \sim U[i,j]$$
積分値と期待値
ここで、任意の閉区間U[i,j]での関数f(x)の積分を以下のように定義します。
$$\theta = \int_{U}f(x)dx$$
この時、f(x)について、期待値は確率密度関数を用いると
$$\begin{align}E[f(X)] &= \int_{U}PDF(x)f(x)dx\end{align}$$
と表されます。
閉区間U[i,j]について、j=i+1の場合、
$$\begin{align}E[f(X)] &= \theta \end{align}$$
が成り立ちます。
大数の弱法則と積分値
閉区間U[i,i+1]からのサンプリングの標本平均を以下のように定義します。
$$\widehat{\theta}=\frac{1}{N}\sum_{k}^n f(X_k)$$
この時、対数の弱法則から
$$\lim_{n \to \infty}P(|\widehat{\theta}-\theta |>\epsilon) = 0 (\forall\epsilon >0)$$
$$\therefore\widehat{\theta}=\frac{1}{N}\sum_{k}^n f(X_k)\simeq\theta$$
よって、閉区間U[i,i+1]における任意の関数f(x)のモンテカルロ積分は
$$ \int_{U}f(x)dx\simeq\frac{1}{N}\sum_{k}^n f(X_k)$$
任意閉区間でのモンテカルロ積分
上記の積分はj=i+1の時に確率密度関数が1になることを利用して積分近似を行なっていましたが、以下で任意閉区間U[i,j]に対するモンテカルロ積分について考えてみます。
確率密度関数を固定したアプローチ
j=i+1となる時、確率密度関数は1となることを利用すれば、
$(X \sim U[i,j])$
のもとで
$$\int_{U}f(x)dx = \sum_{k=i}^{j-1}\int_{k}^{k+1}f(x)dx\quad$$
となります。
確率密度関数が1以外の場合
上記の方法では確率密度関数の値を固定して積分区間の総和を取っていきましたが、確率密度関数の式はわかっているため、そのまま積分することができそうです。
$$\begin{align}\int_{U}f(x)dx&= \int_{U}f(x)\frac{(j-i)}{(j-i)}dx\\&=(j-i)\int_{U}f(x)PDF(x)dx\\& \simeq (j-i)\frac{1}{N}\sum_{k}^n f(X_k)\end{align}$$
Pythonで実装
以下の定積分について考えてみます。
$$\int_{0}^{5}e^{-x}sinxdx$$
Sympyで積分
まずは答えを知るためにsympyを使用して積分します。
import sympy as sp
sp.init_printing(use_unicode=True)
x = sp.symbols('x')
f=sp.sin(x)/sp.exp(x)
print(sp.integrate(f,(x,0,5)).evalf())
>>>
0.502274940083760モンテカルロ積分
与えられた定積分をモンテカルロ積分で解釈すると、以下のようになります。
$$\begin{align}\int_{0}^{5}e^{-x}sinxdx&\simeq \sum_{k=0}^{4}\int_{k}^{k+1}e^{-x}sinxdx\\&=\frac{(5)}{N}\sum_{k}^{n}e^{-x_k}sinx_k\end{align}$$
まずは、n=100の場合について実装します。
import random
n = 100
def int_monte_carlo(n,min,max):
sum = 0
for i in range(n):
x = random.uniform(min,max)
sum+= 1/np.exp(x)*np.sin(x)
return (max-min)/n*sum
print(int_monte_carlo(n,0,5))
>>>
0.5224814048691253実際の値に近いことがわかります。
ここで、
$\lim_{n \to \infty}P(|\widehat{\theta}-E(X)|>\epsilon) = 0 (\forall\epsilon >0)$
が成り立つかどうかを以下のようにして実装して確かめます。
import numpy as np
import sympy as sp
import matplotlib.pyplot as plt
import random
from tqdm import tqdm
fig =plt.figure()
sp.init_printing(use_unicode=True)
x = sp.symbols('x')
m = 10000
f=sp.sin(x)/sp.exp(x)
int_val = sp.integrate(f,(x,0,5)).evalf()
def int_monte_carlo(n,min,max):
sum = 0
for i in range(n):
x = random.uniform(min,max)
sum+= 1/np.exp(x)*np.sin(x)
return (max-min)/n*sum
max = 10000
k = 10
plt.hlines(0,0,max)
for m in tqdm(range(100,max,k)):
int_mt = int_monte_carlo(m,0,5)
plt.scatter(m,int_val-int_mt,s=1)
int_i = [int_monte_carlo(i,0,5)-int_val for i in range(100,max,k)]
i_s = [i*k for i in range(len(int_i))]
#plt.plot(i_s,int_i)
fig.savefig('result.png',dpi=500)実行結果
0に収束していく様子がわかります。
収束へのブレがわかるように、線で繋いでみます。
int_i = [int_monte_carlo(i,0,5)-int_val for i in range(100,max,k)]
i_s = [i*k for i in range(len(int_i))]
plt.plot(i_s,int_i)
最大標本抽出回数を10000回にした場合、以下のような式が成り立っていることがわかります。
$$|\lim_{n \to \infty}P(|\widehat{\theta}-E(X)|>\epsilon)| <0.5\quad (\forall \epsilon >0)$$
もう少し抽出回数を増やした場合について考えてみます。
かなり収束していることがわかります。
$$|\lim_{n \to \infty}P(|\widehat{\theta}-E(X)|>\epsilon)| <0.05\quad (\forall \epsilon >0)$$
