pythonでモンテカルロう!(3章)




 

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

>>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で実装すると以下のようになります。

 

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

例題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種類あります。

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

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

実行結果は以下です

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

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

 

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

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

 

 

 


スポンサーリンク

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

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



スポンサードリンク

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

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

 

 

モンテカルロ法とは?

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

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

円周率の計算

使用したコードです。

 

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

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

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

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

 


 

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

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

 



スポンサードリンク

乱数プロットの画像描画

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

 

  

 

 

 

 

 

 

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

 

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

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

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