[理系大学生応援] モンテカルロ法で円周率を求める

記事更新
2019 Aug. 30: 図を追加

理系大学生(勝手に)応援企画

この記事は、自分が日々Pythonスクリプトを書いている際に気づいた備忘録でありつつも、「理系大学生だけれどもプログラミングの知識あまりなくて、課題提出がツライ」などの学生さんをサポートする記事です。

今回はタイトルにある通り「モンテカルロ法を用いて円周率3.141592...を求めなさい」という課題応援です。

モンテカルロ法とは

確率を用いてなんらかの値を見積もる方法(と私は解釈している)。

円周率を求める(アルゴリズム)

今の場合、$x \in [0, 1]$, $y \in [0, 1]$の領域に点をランダムに描画していき、それが円の4分の1の中に存在する確率から円周率を求めていきます(下図参照)。

円周率を求める(Pythonスクリプト)

pi_estimate.py

#!/usr/bin/env python


import numpy as np


class PiEstimate:
	def __init__(self):
		self.NRAND = 10000

	def pi_estimate(self):
		# set (x, y)
		x = np.random.rand(self.NRAND)
		y = np.random.rand(self.NRAND)
		# compute distance from origin
		dist = x * x + y * y
		# count sum of the number of distance < 1
		under = np.count_nonzero(dist<1)
		# estimate PI
		pi = 4 * under / self.NRAND
		# print result
		print(pi)

	def run(self):
		self.pi_estimate()

if __name__ == '__main__':
	a = PiEstimate()
	a.run()

スクリプト詳細

__init__

self.NRANDでランダムに打つ点の数を指定しています。上述のスクリプトでは10000としていますが、この値を様々に変えて、計算結果がどのように変化するか考察してみるのも面白いでしょう。

pi_estimate

np.random.rand(self.NRAND)

で、0~1のランダムな値が入った要素数self.NRAND個のリストを作成しています。

dist = x * x + y * y

で、その点の原点からの距離の2乗を求めています。

np.count_nonzero(dist<1)

で、条件dist < 1の条件を満たす要素数の合計を計算しています。細かく説明するとnp.count_nonzeroはTrueの要素数を数える関数です。Trueは1, Falseは0として扱われるので、np.sum関数を使用しても同様の結果が得られます。

最後に得られる結果under/self.NRANDは円の4分の1の面積の中に点が入る確率を表していますので、これを4倍することで円周率を求めることができます。

結言

これを応用すれば、様々な積分や物理問題を解くことが可能です。アルゴリズムを理解しいろいろなことに応用していきましょう。何よりこの記事が理系学生のプログラミングの後押しになればこれ幸いです。