2020年4月1日水曜日

モンテカルロ積分

Monte Carlo 積分は乱数を利用して積分を近似する手法です。

「Rによるモンテカルロ法入門」によれば、古典的な「大数の法則」の応用と、重点サンプリングが紹介されています。
前者のもっとも単純な考え方は、「経験平均」を取ること。例えば、確率密度関数から乱数を用いて確率密度を求め、その平均値(高さ)×積分区間(1次元なら長さ、2次元なら面積)で確率(面積や体積)を求める手法です。
数式を変形すれば区間a~bをn区分した際に個々の値を積分するのと同じ形です。

わかりやすいサイト↓
https://www.scratchapixel.com/lessons/mathematics-physics-for-computer-graphics/monte-carlo-methods-in-practice/monte-carlo-integration


で、実施。

ある規則に従った乱数を用意し、

2D KDE をかけて確率密度関数に変換。
さらに81*81のグリッド化してコンター表示。

これを scipy.integrate.nquad で積分すると、
答え:0.9999259355414779
中身は Fortran library の QUADPACK。 21点 Gauss-Kronrod 求積法らしいのですが(by Wikipedia)、遅い。15分かかりました。


これを 81*81=6561データで MC 積分。
答え:0.9754075714917496
甘い。


乱数で x,y を1万データ用意。同様に MC 積分。
答え:0.996469877448263
何度か試しても同程度でした。OKです。
時間は30秒。


6561データを分割。
MC 積分結果は 0.47001458203258695 + 0.505311481041339 = 0.9753 でした。6561データ全体の積分結果とあっています。データを増やすことで精度よく計算可能でしょう。



MC 積分は単純かつコーディングが用意。使えます。

0 件のコメント:

コメントを投稿