モンテカルロ積分[python]

pythonでモンテカルロ積分をしてみたいと思います。
まずは基本的なところで、円の面積を計算してみたいと思います。

モンテカルロ積分は、Pythonのライブラリで関数として定義されていないため、
今回は自分で書いていきます。

モンテカルロ積分は、今回の範囲を円と仮定すると、
ランダム関数の適用範囲と、足しこんでいく範囲（加算範囲）に分けられます。
以下、概念図です。

■ソースコード

import math as mt
import matplotlib.pyplot as plt
import random
import scipy.integrate as si
import numpy as np

list_X = [0]
list_Y = [0]
list_XO = [1]
list_YO = [1]
S = 0
SS = 0
Nmax = 1000000
for i in range(Nmax):
    px = random.uniform(0,1)
    py = random.uniform(0,1)
    if np.abs(mt.sqrt(px**2+py**2))<=1:
        S = S + 1
        list_X.append(px)
        list_Y.append(py)
    else:
        list_XO.append(px) 
        list_YO.append(py) 

#円の1/4の範囲で計算したため、円の面積を出すときは4倍した値を返します。
print(4*S/Nmax)
plt.scatter(list_X,list_Y,color = "b")
plt.scatter(list_XO,list_YO,color = "r")
plt.show()

100万回まわしてみると、3.1415・・・という円周率に近づいていくのが確認できます。
当たり判定も見ていきます。

簡単なモンテカルロ積分を試してみました。
random関数も、範囲を決めて使用できるため、使い勝手が良いと感じました。
また、今回は、matplotlibも使用し、pyplotを使ってみました。

従来だと、CSVやdatファイルに書き出してからgnuplotに食わせていたため、
可視化までの時間がかかっていましたが、ソースコード上から直接グラフを描くことができるので、
手軽さがあっていいと思っています。

by　さくらえび

関連記事：

・積分プログラム(台形積分) — 思い立ったら作ってみる (y-fam.com)