「Rによるモンテカルロ法」という書籍を用いてモンテカルロ法について勉強したことをまとめていこうと思います。
本書に記載されているRのサンプルコードを写経しても意味ないので、pythonでコードを書き直しながら理解を深めていこうと思います。
古典的なモンテカルロ積分
\begin{eqnarray}
E_f[h(X)] = \int_\chi h(x)f(x) dx
\end{eqnarray}
上式の積分の計算をする際に、密度fからサンプル(X1,…..,Xm)を生成し近似として以下の経験平均を算出します。その際にfは様々な確率分布が用いられます。モンテカルロシミュレーションには正規分布や、ポアソン分布、コーシー分布が用いられることが多いようです。
\begin{eqnarray}
\overline{h}(m) = \frac{1}{m}\sum_{j=1}^m h(x_j)
\end{eqnarray}
h(m)は大数の強法則により、サンプル数を増やすとEf[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で実装すると以下のようになります。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 |
#!/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は一様乱数ですが、サンプル(X1,…….Xm)を確率分布から発生させる場合の問題もあります。
練習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種類あります。
- 被積分関数を正規分布とし、コーシー分布からサンプリングする方法
- 被積分関数をコーシー分布とし、正規分布からサンプリングする方法
試しに両方の方法で解いてみて、近似解の精度をみてみましょう。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 |
# 正規・コーシーベイズ推定量 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) |
実行結果は以下です
1 2 3 |
モンテカルロ積分(コーシー): 3.44162131364 モンテカルロ積分(正規): 3.42786016038 integ関数: 3.435061555229311 |
コーシー分布よりも、正規分布でサンプリング(正規近似)した方が精度は良いみたいです。
コーシー分布は外れ値をとる可能性が大きく、裾が広いのでその影響が出ています。
>>コーシー分布と正規分布についてはこちら
今回の被積分関数のプロット図です。青線が\(\delta(x)\)の分子、緑線が\(\delta(x)\)の分母です。
スポンサーリンク