モンテカルロ法

モンテカルロ法で円周率を求めよう というソレを見つつ Python 実装をヒリ出してみることに。本当はお仕事対応しなきゃ、なんですがもうへろへろなのでへろへろついでに、ってことで (ぇ

1/4 円な扇形

一辺の長さが 1 の正方形の中に 1/4 円な扇形があるとして、

正方形の面積 : 扇形の面積 = 1 : π/4

とのこと。n 個落として r 個が扇形の中だとすると π は以下で求められるのか。

π = 4 * r/n

つうか

扇形の中かどうか、は半径 1 なので x 座標の二乗と y 座標の二乗を加えた値が 1 より小さいかどうかで判定すれば良いのか。
で、空行入れててハマりましたが以下がでっち上がってます。

import random

def square(n):
    return n * n

def monte_carlo_sim_pi(num):
    r = 0
    for i in range(num):
        x = random.random()
        y = random.random()
        if 1 >= (square(x) + square(y)):
            r += 1
    return 4.0 * r / num

実行してみたら以下な出力。

>>> import pi_simulate
>>> monte_carlo_sim_pi(100)
3.04
>>> monte_carlo_sim_pi(1000)
3.144
>>> monte_carlo_sim_pi(10000)
3.146
>>> monte_carlo_sim_pi(100000)
3.13792
>>> monte_carlo_sim_pi(1000000)
3.139704
>>> monte_carlo_sim_pi(10000000)
3.141104

リキがあれば演習問題 5-11 に着手の方向にて。