確率論的に円周率πを推定してみた

久しぶりの投稿です。
Atsumaruエンジニアの渡邊です。

最近は新規プロジェクトの立ち上げ・ドライブが増えて、なかなかコードを触れていない日々が続いています><

今回は、円周率の求め方で面白いアプローチがあることを知ったので、自分でもやってみた話を書いてみようと思います。

考え方

xy座標を思い浮かべてください。

f:id:snoopy_no_sora:20200226141704p:plain

ここに、原点Oを中心とする、半径rの円（Cとします）と、一辺が2rの正方形（Sとします）を描きます。

f:id:snoopy_no_sora:20200226141708p:plain

Sの中に、ひたすらランダムで点（dot）を打ち続けます。

f:id:snoopy_no_sora:20200226141924p:plain

Sの中を埋め尽くすほど点を打ったとき、「Sの中にある点の数：Cの中にある点の数 ≒ Sの面積：Cの面積」と言えます。

f:id:snoopy_no_sora:20200226141637p:plain

Sの中にある点の数を「s_dot_num」、Cの中にある点の数を「c_dot_num」とおくと、

(c_dot_num) / (s_dot_num)
　= Cの面積 / Sの面積
　= (π * rの2乗) / (2r)の2乗
　= π / 4

となることから、

π = 4 * (c_dot_num) / (s_dot_num)

と定義ができました。
あとは、c_dot_numとs_dot_numをそれぞれ求めることができれば、πの値を推定できます。
（rの値は関係ないようですね）

コード

今回はpythonで書きました。

# 一辺の長さが1の正方形内にランダムにドットを打ち、円周率を推定する。
import random

def estimate_pi(n):
    c_dot_num = 0
    s_dot_num = 0
    for _ in range(n):
        x = random.uniform(0, 1)
        y = random.uniform(0, 1)
        distance = x**2 + y**2
        if distance <= 1:
            c_dot_num += 1
        s_dot_num += 1
    pi = 4 * c_dot_num / s_dot_num
    return pi

if __name__ == '__main__':
    print('dotが100個のとき、円周率は{}'.format(estimate_pi(100)))
    print('dotが1000個のとき、円周率は{}'.format(estimate_pi(1000)))
    print('dotが10000個のとき、円周率は{}'.format(estimate_pi(10000)))
    print('dotが100000個のとき、円周率は{}'.format(estimate_pi(100000)))