モンテカルロ法で円周率を求める - Pythonでいろいろやってみる

モンテカルロ法で円周率を求めます。
モンテカルロ法はシミュレーションや数値計算を乱数を用いて行う手法の総称です。
図のように一辺がrの正方形とそれに内接する半径rの4分割円があり、その正方形の中にランダムに点を打ちます。すべての点の数(青い点と赤い点の合計)と4分割円の中の点の数(青い点)から、4分割円の面積が
4分割円の面積 = 4分割円の中の点の数 ÷ すべての点の数 × 正方形の面積( r × r )
と近似されます。4分割円の面積は円周率をπとすると
4分割円の面積 = π × r × r ÷ 4
となり、以上より
4分割円の中の点の数 ÷ すべての点の数 × 正方形の面積( r × r ) = π × r × r ÷ 4
π = 4分割円の中の点の数 ÷ すべての点の数 × 4
で円周率の近似値が求められます。
f:id:T_A_T:20201215181706p:plain

環境

windows10 home
Anaconda 3/ jupyter notebook 5.6.0
Python 3.7.0

コード

正方形の一辺および円の半径を1として、点のx座標、y座標として、0以上1未満の乱数を生成するrandom.random()関数を用いてランダムな値を与えます。求めた座標が4分割円の中にあるかどうかは原点(0,0)と座標(x,y)の距離が半径r(ここでは1)以下かどうかで判定します。試行回数を10回、100回・・・10000000回と変え、試行回数により求まる円周率がどのように変わるかを確認します。

import random
import math

# モンテカルロ法での円周率計算関数
def paicalc(testn):
    incircle = 0  # 点が円内の場合の集計用変数
    for i in range(testn):
        x = random.random()  # 0以上1未満の乱数
        y = random.random()  # 0以上1未満の乱数
        distance = math.sqrt(x**2+y**2)  # (0,0)と(x,y)の距離
        if distance<=1:  # 距離が1以下なら(円の中なら)
            incircle += 1  # 点が円内ならカウント        
    print('試行回数:',testn,' 円周率:',incircle/testn*4)  # 試行回数と円周率の表示 

tests = [10**(i+1) for i in range(7)]  # 試行回数のリスト
for j in tests:  
    paicalc(j)    # paicalc呼び出し

実行結果

試行回数が増えるにつれ値が3.141592・・に近づくことが分かります。

試行回数: 10 円周率: 2.8
試行回数: 100 円周率: 3.16
試行回数: 1000 円周率: 3.204
試行回数: 10000 円周率: 3.1296
試行回数: 100000 円周率: 3.15456
試行回数: 1000000 円周率: 3.142624
試行回数: 10000000 円周率: 3.142068