「Rによるモンテカルロ法」という書籍を用いてモンテカルロ法について勉強したことをまとめていこうと思います。

本書に記載されているRのサンプルコードを写経しても意味ないので、pythonでコードを書き直しながら理解を深めていこうと思います。

古典的なモンテカルロ積分

\begin{eqnarray}
E_f[h(X)] = \int_\chi h(x)f(x) dx
\end{eqnarray}

上式の積分の計算をする際に、密度fからサンプル（X₁,…..,X_m）を生成し近似として以下の経験平均を算出します。その際にfは様々な確率分布が用いられます。モンテカルロシミュレーションには正規分布や、ポアソン分布、コーシー分布が用いられることが多いようです。

\begin{eqnarray}
\overline{h}(m) = \frac{1}{m}\sum_{j=1}^m h(x_j)
\end{eqnarray}

h(m)は大数の強法則により、サンプル数を増やすとE_f[h(x)]に収束します。

試しに例3.3を解いてみます。

例3.3

次の関数を一様乱数を用いて、独立同分布の乱数を発生させて、\(\sum h(Ui)/n\)によって\(\int h(x) dx\)を近似することができます。

\begin{eqnarray}
h(x) = \{\cos(50x)+ \sin(20x)\}^2
\end{eqnarray}

pythonで実装すると以下のようになります。

#!/usr/bin/env python
# -*- coding: utf-8
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad
import random


integ = []       #任意の積分範囲の積分値を格納する配列
std_list = []   #各回数（i）における標準誤差を格納する配列
a = []          #各回数における積分値の平均値を格納する配列
error1 = []     #95％信頼区間の値(+)
error2 = []     #95％信頼区間の値(-)
for i in range (1,2001):
    j = random.random()
    integ.append((np.cos(50*j)+np.sin(20*j))**2)
    std  = np.std(integ)/np.sqrt(i)
    std_list.append(std)
    a.append(np.average(integ))
    error1.append(a[i-1] + 2 * std_list[i-1])
    error2.append(a[i-1] - 2 * std_list[i-1])

#integ関数による積分
y ,error = quad(lambda x : (np.cos(50*x)+np.sin(20*x))**2, 0, 1)
print(y)

#グラフ描画
plt.figure(figsize=(7, 4))
plt.xlabel("number")
plt.ylabel("value")
plt.plot(error1 ,color = "r", label = "+0.05CI")
plt.plot(error2 , color = "g", label = "-0.05CI")
plt.plot(a, color = "b", label = "average")
plt.show()

#!/usr/bin/env python

# -*- coding: utf-8

import numpy as np

import matplotlib.pyplot as plt

from scipy.integrate import quad

import random

integ = [] #任意の積分範囲の積分値を格納する配列

std_list = [] #各回数（i）における標準誤差を格納する配列

a = [] #各回数における積分値の平均値を格納する配列

error1 = [] #95％信頼区間の値(+)

error2 = [] #95％信頼区間の値(-)

for i in range (1,2001):

j = random.random()

integ.append((np.cos(50*j)+np.sin(20*j))**2)

std = np.std(integ)/np.sqrt(i)

std_list.append(std)

a.append(np.average(integ))

error1.append(a[i-1] + 2 * std_list[i-1])

error2.append(a[i-1] - 2 * std_list[i-1])

#integ関数による積分

y ,error = quad(lambda x : (np.cos(50*x)+np.sin(20*x))**2, 0, 1)

print(y)

#グラフ描画

plt.figure(figsize=(7, 4))

plt.xlabel("number")

plt.ylabel("value")

plt.plot(error1 ,color = "r", label = "+0.05CI")

plt.plot(error2 , color = "g", label = "-0.05CI")

plt.plot(a, color = "b", label = "average")

plt.show()

＞＞描画結果と詳しい解説は以下の記事

例題3.3は一様乱数ですが、サンプル(X₁,…….X_m)を確率分布から発生させる場合の問題もあります。

練習3.1

以下のコーシーベイズ推定量についてx = 4のときを想定して解いてみます。

\begin{eqnarray}
\delta(x) = \frac{ \int_{-\infty}^{\infty} \frac{\theta}{1 + \theta^2} e^{-(x-\theta)^2/2} d\theta }{\int_{-\infty}^{\infty} \frac{1}{1 + \theta^2} e^{-(x-\theta)^2/2} d\theta}
\end{eqnarray}

このままでは確率分布が使えないので、式変形をします。

\begin{eqnarray}
\delta(x) &=& \frac{ \int_{-\infty}^{\infty} \frac{\theta}{1 + \theta^2} e^{-(x-\theta)^2/2} d\theta }{\int_{-\infty}^{\infty} \frac{1}{1 + \theta^2} e^{-(x-\theta)^2/2} d\theta}\\
&=& \frac{\int_{-\infty}^{\infty} \theta \frac{1}{1+\theta^2} \frac{1}{\sqrt{2\pi}} \frac{\theta}{1 + \theta^2} e^{-(\theta – x)^2 /2} d\theta }{\int_{-\infty}^{\infty} \frac{1}{1+\theta^2} \frac{1}{\sqrt{2\pi}} \frac{1}{1 + \theta^2} e^{-(\theta – x)^2 /2} d\theta}
\end{eqnarray}

こうすると、正規分布とコーシー分布が出てきます。

解き方には2種類あります。

被積分関数を正規分布とし、コーシー分布からサンプリングする方法
被積分関数をコーシー分布とし、正規分布からサンプリングする方法

試しに両方の方法で解いてみて、近似解の精度をみてみましょう。

# 正規・コーシーベイズ推定量
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad
from scipy.stats import cauchy,norm

x=4
N = 100000

#被積分関数のプロット
x_ra = np.arange(-1, 9, 0.01)
y_1 = x_ra*np.exp(-((x - x_ra )**2)/2)/(1+x_ra**2)
y_2 = np.exp(-((x - x_ra )**2)/2)/(1+x_ra**2)
#グラフ描画
plt.figure(figsize=(7, 4))
plt.xlabel("x")
plt.ylabel("P(x)")
plt.plot(x_ra, y_1, color = "b", label = "average")
plt.plot(x_ra, y_2, color = "g", label = "average")
plt.show()


#モンテカルロ積分(コーシー分布からサンプリングする)

#分子をモンテカルロ積分
hoge = cauchy.rvs(size=N)
I_nume = np.mean(hoge * norm(loc=x).pdf(hoge))

#分母をモンテカルロ積分
hoge = cauchy.rvs(size=N)
I_denom = np.mean(norm(loc=x).pdf(hoge))

I = I_nume/I_denom
print(u"モンテカルロ積分(コーシー)：",I)


#モンテカルロ積分(正規分布からサンプリングする)

#分子をモンテカルロ積分
hoge = norm(loc=x).rvs(size=N)
I_nume = np.mean(hoge * cauchy.pdf(hoge))

#分母をモンテカルロ積分
hoge = norm(loc=x).rvs(size=N)
I_denom = np.mean(cauchy.pdf(hoge))

I = I_nume/I_denom
print(u"モンテカルロ積分(正規)：",I)


#通常のinteg関数による積分
f_1 = lambda t: t * norm(loc=x).pdf(t) * cauchy.pdf(t)
f_2 = lambda t: norm(loc=x).pdf(t) * cauchy.pdf(t)
I_bunsi = quad(f_1, -np.inf, np.inf)[0]
I_bunbo = quad(f_2, -np.inf, np.inf)[0]
print(u"integ関数：",I_bunsi/I_bunbo)

# 正規・コーシーベイズ推定量

import numpy as np

import matplotlib.pyplot as plt

from scipy.integrate import quad

from scipy.stats import cauchy,norm

x=4

N = 100000

#被積分関数のプロット

x_ra = np.arange(-1, 9, 0.01)

y_1 = x_ra*np.exp(-((x - x_ra )**2)/2)/(1+x_ra**2)

y_2 = np.exp(-((x - x_ra )**2)/2)/(1+x_ra**2)

#グラフ描画

plt.figure(figsize=(7, 4))

plt.xlabel("x")

plt.ylabel("P(x)")

plt.plot(x_ra, y_1, color = "b", label = "average")

plt.plot(x_ra, y_2, color = "g", label = "average")

plt.show()

#モンテカルロ積分(コーシー分布からサンプリングする)

#分子をモンテカルロ積分

hoge = cauchy.rvs(size=N)

I_nume = np.mean(hoge * norm(loc=x).pdf(hoge))

#分母をモンテカルロ積分

hoge = cauchy.rvs(size=N)

I_denom = np.mean(norm(loc=x).pdf(hoge))

I = I_nume/I_denom

print(u"モンテカルロ積分(コーシー)：",I)

#モンテカルロ積分(正規分布からサンプリングする)

#分子をモンテカルロ積分

hoge = norm(loc=x).rvs(size=N)

I_nume = np.mean(hoge * cauchy.pdf(hoge))

#分母をモンテカルロ積分

hoge = norm(loc=x).rvs(size=N)

I_denom = np.mean(cauchy.pdf(hoge))

I = I_nume/I_denom

print(u"モンテカルロ積分(正規)：",I)

#通常のinteg関数による積分

f_1 = lambda t: t * norm(loc=x).pdf(t) * cauchy.pdf(t)

f_2 = lambda t: norm(loc=x).pdf(t) * cauchy.pdf(t)

I_bunsi = quad(f_1, -np.inf, np.inf)[0]

I_bunbo = quad(f_2, -np.inf, np.inf)[0]

print(u"integ関数：",I_bunsi/I_bunbo)

実行結果は以下です

モンテカルロ積分(コーシー)： 3.44162131364
モンテカルロ積分(正規)： 3.42786016038
integ関数： 3.435061555229311

モンテカルロ積分(コーシー)： 3.44162131364

モンテカルロ積分(正規)： 3.42786016038

integ関数： 3.435061555229311

コーシー分布よりも、正規分布でサンプリング（正規近似）した方が精度は良いみたいです。

コーシー分布は外れ値をとる可能性が大きく、裾が広いのでその影響が出ています。

＞＞コーシー分布と正規分布についてはこちら

今回の被積分関数のプロット図です。青線が\(\delta(x)\)の分子、緑線が\(\delta(x)\)の分母です。

記事が役に立ったらシェア！

スポンサードリンク

モンテカルロ法の勉強の一貫として、円周率を導出するプログラムを書いたので載せます。

今回はこちらの本を用いて勉強しました。数値計算の基本的な手法について網羅されております。

モンテカルロ法とは？

乱数を用いてサンプル分布の情報を得る手法の総称です。乱数を用いていれば、その乱数の分布がなんであろうが、全てモンテカルロ法になるわけです。

かなり広義な意味です。さらにその乱数の確率分布に正規分布とかベルヌーイ分布、コーシー分布、指数分布とかがあるわけです。

円周率の計算

使用したコードです。

#!/usr/bin/env python
# -*- coding: utf-8
import math
import random
import numpy as np
import matplotlib.pyplot as plt
import pylab

result=[]　
result_diff=[]
for j in range(500):#円周率計算を500回繰り返す
    circle=[]
    hoge=[]
    for i in range(10000):　#10000個のデータ点作成
        x_rand = random.random()
        y_rand = random.random()
        hoge.append([x_rand,y_rand])
        if 0 <= (x_rand)**2+(y_rand)**2 <= 1:#半径1の円内にある時は別でカウントする
            circle.append([x_rand,y_rand])

    count = len(circle)
    world = len(hoge)
    ratio = 4*count/world
    difference = ratio - math.pi
    result.append(ratio)
    result_diff.append(difference)
print(sum(result)/500)#500回の円周率の平均を求める
pylab.plot(result)#円周率の500回の推移を描画する
pylab.show()

#円周率と計算値の差のヒストグラム
x = np.arange(-0.6,0.6,0.01)
plt.xlabel("pi-result")
plt.ylabel("N")
plt.hist(result_diff,bins=40,range=(-0.055,0.055))
plt.show()

#!/usr/bin/env python

# -*- coding: utf-8

import math

import random

import numpy as np

import matplotlib.pyplot as plt

import pylab

result=[]　

result_diff=[]

for j in range(500):#円周率計算を500回繰り返す

circle=[]

hoge=[]

for i in range(10000):　#10000個のデータ点作成

x_rand = random.random()

y_rand = random.random()

hoge.append([x_rand,y_rand])

if 0 <= (x_rand)**2+(y_rand)**2 <= 1:#半径1の円内にある時は別でカウントする

circle.append([x_rand,y_rand])

count = len(circle)

world = len(hoge)

ratio = 4*count/world

difference = ratio - math.pi

result.append(ratio)

result_diff.append(difference)

print(sum(result)/500)#500回の円周率の平均を求める

pylab.plot(result)#円周率の500回の推移を描画する

pylab.show()

#円周率と計算値の差のヒストグラム

x = np.arange(-0.6,0.6,0.01)

plt.xlabel("pi-result")

plt.ylabel("N")

plt.hist(result_diff,bins=40,range=(-0.055,0.055))

plt.show()

擬似乱数を用いていますので、プログラムを動かすたびに違う結果になります。

ちなみに500回円周率を計算した時の平均値は3.141804799999997になりました。グラフから見てもあってそうです。

さすがにデータ点10000点×500回なので、ほぼ真の値に収束しています。

計算回数を増やせば増やすほど、値が収束する現象は「大数の法則」として知られています。

各回数での円周率の計算値がどれくらい真の値からずれているのか、ずれに頻度はどれくらいかを視覚化するためヒストグラムを作成しました。

横軸が（円周率-計算で求めた円周率）、縦軸が各範囲における頻度です。今回はビン数を40（x軸の範囲内で何等分するか）にしました。

スポンサードリンク

乱数プロットの画像描画

モンテカルロ法のwiki風な画像を作成しました。試しにデータ点が10000と1000の時で計算しました。左側がデータ点10000、右側がデータ点1000の時です。半径1の円内にあるものが赤、それ以外が青です。

上の図を作成するために使用したコードを載せます。

#!/usr/bin/env python
# -*- coding: utf-8
import math
import random
import numpy as np
import matplotlib.pyplot as plt
import pylab


x_list=[]
y_list=[]
x_list2=[]
y_list2=[]
for i in range(10000):
    x_rand = random.random()#データ点のx座標作成
    y_rand = random.random()#データ点のy座標作成
    if 0 <= (x_rand)**2+(y_rand)**2 <= 1:　#半径1の円内にある時
        x_list2.append(x_rand)
        y_list2.append(y_rand)
    else:#それ以外
        x_list.append(x_rand)
        y_list.append(y_rand)

#グラフ作成
pylab.figure(figsize=(7, 7))#グラフを正方形にする
plt.title("n=10000")
plt.scatter(x_list2,y_list2,c="red",marker=".")
plt.scatter(x_list,y_list,c="blue",marker=".")
plt.show()

#!/usr/bin/env python

# -*- coding: utf-8

import math

import random

import numpy as np

import matplotlib.pyplot as plt

import pylab

x_list=[]

y_list=[]

x_list2=[]

y_list2=[]

for i in range(10000):

x_rand = random.random()#データ点のx座標作成

y_rand = random.random()#データ点のy座標作成

if 0 <= (x_rand)**2+(y_rand)**2 <= 1:　#半径1の円内にある時

x_list2.append(x_rand)

y_list2.append(y_rand)

else:#それ以外

x_list.append(x_rand)

y_list.append(y_rand)

#グラフ作成

pylab.figure(figsize=(7, 7))#グラフを正方形にする

plt.title("n=10000")

plt.scatter(x_list2,y_list2,c="red",marker=".")

plt.scatter(x_list,y_list,c="blue",marker=".")

plt.show()

本当はwikiみたいにデータ点数を変えた時のプロット描画を見たかったので、for文を回して各データ点数ごとにグラフを作成し、それを動画にするプログラムも書きたかったのですが、疲れたので今回はここまでにします。

wiki風動画作成は別途記事にしようと思います。

記事が役に立ったらシェア！

Tag: モンテカルロ法

pythonでモンテカルロう！（3章）

モンテカルロ法で円周率を計算する。（python）

モンテカルロ法とは？

円周率の計算

乱数プロットの画像描画