モンテカルロ積分 (2)

モンテカルロ積分（2014/7/28）のつづき。今回は、Rによるモンテカルロ法入門の例3.3と練習問題3.1のモンテカルロ積分の例を検証した。

例3.3

例3.3は簡単な積分の問題。

これも前回と同じように確率分布 p(x) が明示的にないため一様分布を追加する。積分区間が [0, 1] なので b - a = 1 となり無視できる。

実行結果は

scipy.integrate: 0.96520093605
モンテカルロ積分: 0.966150133367

練習問題3.1 正規・コーシー-ベイズ推定量の積分

練習問題3.1はもっと複雑な積分だ。どうやら正規・コーシー-ベイズ推定量と呼ばれるものらしい。何に使うのかはよくわからないけど・・・

これがなぜ正規・コーシー-ベイズ推定量と呼ばれるのか謎だったのだが、上の式を下の式のように変形すると正規分布とコーシー分布が出てくることがわかった。が突然現れるけど分子と分母でクリアされるので追加しても問題ない。

分子と分母に出てくる積分をモンテカルロ積分で計算してみる。分子と分母はほとんど同じ形なので分子だけ説明。この積分は二通りの方法でモンテカルロ積分が適用できることがわかる。

(1) コーシー分布をpとみなしてコーシー分布からサンプリングする場合

(2) 正規分布をpとみなして正規分布からサンプリングする場合

初めて一様分布以外からサンプリングするモンテカルロ積分の例が出てきた。両方のパターンで積分を計算し、結果を比較してみよう。ついでにscipy.integrateの結果も出してみた。

scipy.integrate: 3.43506155523
モンテカルロ積分(1): 3.4456557745
モンテカルロ積分(2): 3.43354612482

どちらの場合でもscipy.integrateとほぼ同じ結果が得られることがわかる。何回かやるたびに結果が変わるのだけれど正規分布の方がscipy.integrateの結果に近いことが多いようだ。コーシー分布は正規分布より裾野が広いため平均から外れたとんでもないサンプルが出てしまうことがあるからかな？特にサンプリング数Nが小さいときに違いが出てくる。どちらの分布からサンプリングした方がよいかの検証は次回やってみよう。