この記事は何
3 行の関数で円周率を求める方法とその解説です。
ルールとして無理矢理コードを圧縮するのではなく、あくまで自然なやり方で書くこととします。
バージョン情報
python: 3.7.3
numpy: 1.17.4
コード
calcPi.py
import numpy as np
def calcPi(n):
points = np.random.uniform(0.0, 1.0, (n, 2))
inner = np.hypot(points[:,0], points[:,1]) < 1
return inner.mean() * 4
※最後にmean()を使うやり方は@konandoiruasaさんのコメントを参考にしました。
解説
モンテカルロ法で円周率を求める方法自体は既知とします。
1行目
np.random.uniformにより 0 から 1 の範囲の一様乱数を n 行 2 列の array に格納しています。
各行が平面上の点に対応しています。
1行目
points = np.random.uniform(0.0, 1.0, (n, 2))
"""
例えば n=10 のとき
[[ 0.6906296 0.20549271]
[ 0.13386813 0.77204275]
[ 0.5970941 0.49659941]
[ 0.92884413 0.37740529]
[ 0.49212498 0.13915062]
[ 0.69357975 0.23229706]
[ 0.14287715 0.14076891]
[ 0.20199753 0.49663344]
[ 0.90105166 0.87612407]
[ 0.19636323 0.39813228]]
"""
2行目
np.hypot(x, y)は点(x, y)の原点からの距離を返します。
points[:,0]はpointsの 0 列目を抽出した array で、つまりは各点のx座標の大きさです。
points[:,1]も同様に各点のy座標になっています。
最後に< 1とやっているので、arrayの中身が1未満ならTrue, 1以上ならFalseに変わります。
つまり原点からの距離が1未満ならTrue, 1以上ならFalseになっている array が得られます。
2行目
inner = np.hypot(points[:,0], points[:,1]) < 1
"""
例えば、n=10のとき、
inner == [True True True True False True False True True False]
"""
3行目
原点からの距離が1未満の点の個数を、点の総数で割るとπ/4の近似値が求まるため、
innerの平均を取り4倍したものを返します。
（コメントにもある通りTrue, Falseはそれぞれ計算時1, 0に変換されます）
4行目
return inner.mean() * 4
実行結果
seedを 0 にした場合の実行結果です。
print(calcPy(10)) # => 2.8
print(calcPi(100)) # => 3.32
print(calcPi(1000)) # => 3.302
print(calcPi(10000)) # => 3.1544
print(calcPi(100000)) # => 3.13228
print(calcPi(1000000)) # => 3.142204
print(calcPi(10000000)) # => 3.1421468
print(calcPi(100000000)) # => 3.14170808
以上です。間違い等がありましたら教えてください。