【ベイズ統計学再入門③】モンテカルロ法

復習がてら、ベイズ統計の解説とPythonでの実装を行います。
第三回は「モンテカルロ法」についてです。

モンテカルロ法とは

まずはモンテカルロ法とは何かについて説明します。

例えば以下の一辺が1の正方形の領域にランダムに点を打って円周率を求めることを考えます。

ランダムに打った点の数を$n$，円の中に入った点の数を$p$とおいて

\frac{p}{n}

を計算します。ランダムに打つ点を大きくしたとき、この計算は$\frac{\pi}{4}$（半径1の円の面積の4分の1）に近づくはずです。

\frac{p}{n}\approx \frac{\pi}{4}

このように、乱数を使った数値計算法のことをモンテカルロ法といいます。

実装

モンテカルロ法も用いて円周率の近似値を求めてみます。

モンテカルロ法のステップを変えるスライドするバーを作成し、インタラクティブな描画ができるようなプログラムを書いてみます。

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
# インタラクティブな描画
from IPython.html.widgets import interact
# グラフを埋め込むために以下を入力
%matplotlib inline

# ggplotというスタイルでグラフを表示する（お好みで）
plt.style.use("ggplot")
# 乱数のseedを固定
np.random.seed(123)

NMC = 100
xmc = np.random.rand(NMC)
ymc = np.random.rand(NMC)
# スライドするバーの設定（0からNMCまで1刻み）
@interact(mcs=(0,NMC,1))
def animation(mcs=0):
    plt.figure(figsize=(6,6))
# グラフの描画範囲を明示的に指定
    plt.xlim([0,1])
    plt.ylim([0,1])
# 円を描画
    x = np.arange(0,1,0.001)
    y = (1-x**2)**0.5
    y2 = np.ones(x.shape[0])#すべての要素が1のarray
    plt.plot(x,y)
# 色を塗りつぶす（α値で透明に）
    plt.fill_between(x,y,alpha=0.3)
    plt.fill_between(x,y,y2,alpha=0.3)
    r = (xmc[:mcs]**2 + ymc[:mcs]**2)**0.5
# 円の中に入ったかどうかの判定（r<=1ならば1，それ以外なら0）
    accept = np.where(r<=1, 1, 0)
    accept_ratio = np.sum(accept) / mcs
    plt.scatter(xmc[:mcs], ymc[:mcs], color="black", marker=".")
    plt.show()
    print("Monte Carlo: ",accept_ratio)
    print("pi/4: ",np.pi / 4.0)

100程度では精度が微妙ですね。点を打つ数を増やせばもっと近似してくれますが、ランダムに点を打つので、目的の値に近づいたり遠ざかったりとふらついて近づくことが知られています。

次回は棄却サンプリングやMCMCについて書きたいと思います。