Pythonを使って円周率を求める #モンテカルロ法

13431 ワード
モンテカルロ法 Python matplotlib Python テキストリンク

初めに

この記事では
ラムダさんの動画を見てそのまま、深夜テンションで書いたモンテカルロ法を用いて円周率を求めるプログラムを解説する記事です。深夜テンションなのであまり期待しないで下さい！

コードが酷い

import matplotlib.pyplot as plt
import random

Xlist=[]
Ylist=[]
Number_list=[]
pi=0
pi_list=[]

for i in range(15000):
    Xlist.append(random.uniform(0,1))
    Ylist.append(random.uniform(0,1))
    if (i+1)%10 == 0:
        Number_list.append(i+1)
        pi_list.append((pi/(i+1))*4)
    if Xlist[i]**2+Ylist[i]**2 <= 1 and (i+1)%10 == 0:
        Number_list.append(i+1)
        pi_list.append((pi/(i+1))*4)
        pi+=1
    elif Xlist[i]**2+Ylist[i]**2 <= 1:
        pi+=1

print((pi/(i+1))*4)
plt.plot(Number_list,pi_list)
plt.show()


何これ汚ない
一応動くちゃ動くけど、同じ処理が何回も書かれいるし流石にやばいので結局朝に書き直しました。

import matplotlib.pyplot as plt
import random

Xlist=[]
Ylist=[]
Number_list=[]
pi=0
pi_list=[]

for i in range(15000):
    Xlist.append(random.uniform(0,1))
    Ylist.append(random.uniform(0,1))
    if (i+1)%10 == 0:
        Number_list.append(i+1)
        pi_list.append((pi/(i+1))*4)
    if Xlist[i]**2+Ylist[i]**2 <= 1:
        pi+=1

print((pi/(i+1))*4)
plt.plot(Number_list,pi_list)
plt.show()

少しはマシになりましたね。
それでは1行ずつコードを解説して行きましょう。

import matplotlib.pyplot as plt
import random

一行目ではmatplotlibと言うグラフを作図する時に便利なライブラリーをpltと言う名前で使える様にしました。
二行目ではrandomと言うランダムな値を生成できるライブラリーをインポートしています。

Xlist=[]
Ylist=[]
Number_list=[]
pi=0
pi_list=[]

ここでは諸々円周率の計算に必要な変数を宣言をしています。
XlistとYlistにはランダムに生成した点の座標を格納して後の計算で、点が円内部に入っているか否かを判定します。Number_listではグラフ作成に必要な等差数列を作成しています。本当はnumpyを使って先に等差数列を作った方が早いですが、今回は方法を調べるのがめんどくさいのでパスしました(サボり魔)。piには合計で円内部にプロットされた点の数の合計をpi_listにはそれぞれ計算で求めた円周率ををリスト形式で保存する時に使います。

for i in range(15000):
    Xlist.append(random.uniform(0,1))
    Ylist.append(random.uniform(0,1))
    if (i+1)%10 == 0:
        Number_list.append(i+1)
        pi_list.append((pi/(i+1))*4)
    if Xlist[i]**2+Ylist[i]**2 <= 1:
        pi+=1

ここではfor文を使って以下の処理を15,000回繰り返します(人間がやったら死んじゃう)。
XlistとYlistにramdomモジュールを使って点の座標を格納します。
説明の都合先にif Xlist[i]2+Ylist[i]2 <= 1:のブロックを解説します。
ここでは面倒臭い説明は省きますが
$x^2+y^2=r^2$
より取り敢えずXlistとYlistの最後に格納されている数をそれぞれ二乗して足して半径の二乗以下だったら、円の中に点が入っている事が分かります。これを数式として表したのが上のif文です。このif文で点が円の中にいる事が分かったので、piに1加えます。最後にif (i+1)%10 == 0:のブロックについてです。ここでは十回に一回円周率を計算してpi_listに格納し、ついでにNumber_listを等差数列にするために数を追加します。

plt.plot(Number_list,pi_list)
plt.show()
これはグラフを作成するだけのプログラムですので解説はしません(手抜き)。

終わりに

多分いらないと思いますが、PDFにして今回の内容をまとめて見ました。
プログラム歴はびっくりするほど浅いので間違っていたら気軽にコメントで教えてください！！
後よかったらコメントもくだい。