円周率をモンテカルロ法で

文科省教材：確率モデルのシミュレーションでも書いたが，$[0, 1)$ の一様乱数 $x$, $y$ が $x^2 + y^2 < 1$ を満たす割合を4倍すれば円周率 $\pi$ の概数が求められる：

import numpy as np

rng = np.random.default_rng(20210729)
totalcount = 10000

def pi1():
    x = rng.random(totalcount)
    y = rng.random(totalcount)
    incount = np.count_nonzero(x**2 + y**2 < 1)
    return incount * 4 / totalcount

1000回やってみて，平均と標準偏差を求める：

p1 = [pi1() for _ in range(1000)]
np.mean(p1), np.std(p1)

これについて，三重大学の白井先生がおもしろい話を教えてくださった：

円周率を求める単純抽出のモンテカルロ法をちょっといじると精度がよくなるという話をしました（理由は「そりゃそうか」というものです） pic.twitter.com/Sy9TagyCSd

— Nobu C. Shirai (@nobucshirai) 2021年7月28日

そりゃそうだ。で，上のコードをできるだけ生かして，次のようにしてみた。

def pi2():
    x = rng.random(totalcount)
    y = rng.random(totalcount)
    incount = np.count_nonzero((x**2 + y**2 < 1) & ((1-x)**2 + (1-y)**2 < 1))
    return incount * 2 / totalcount + 2

p2 = [pi2() for _ in range(1000)]
np.mean(p2), np.std(p2)

確かに標準偏差が小さくなっている。

[蛇足] TeX で乱数を使う では☃でプロットしている。

Last modified: