LoginSignup
4

More than 1 year has passed since last update.

posted at

updated at

Organization

3 行の関数で円周率を求める【Python・モンテカルロ法】

この記事は何

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

以上です。間違い等がありましたら教えてください。

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
What you can do with signing up
4