蒙特卡洛方法——玩具例子(1)圆周率

近似计算圆周率

正方形内部有一个相切的圆，它们的面积之比是π/4。

由圆面积和正方形面积公式可得：

$\frac{S_{圆面积}}{S_{正方形}}=\frac{\pi r^2}{(2r)^2}=\frac{\pi}{4}$

现在向上面的正方形区域均匀的随机撒点，通过计算落在圆内的点数与总的点数来估计$\pi$

# -*- coding:utf8 -*-

import numpy as np
import matplotlib.pyplot as plt

def pi_exp(R, N):
    x = np.random.uniform(-R * 1.0, R * 1.0, (N,))
    y = np.random.uniform(-R * 1.0, R * 1.0, (N,))
    incircle = np.sum(np.power(x, 2) + np.power(y, 2) <= (R * 1.0)**2)
    error = np.pi - (4.0 * incircle / N)
    print error
    ######################################
    theta = np.linspace(0, 2*np.pi, 3600)
    cx = R * np.cos(theta)
    cy = R * np.sin(theta)
    plt.plot(cx,cy, color="r")
    ######################################
    rx = np.asarray([-1, 1, 1, -1, -1] *R, dtype = np.float)
    ry = np.asarray([1, 1, -1, -1, 1] *R, dtype = np.float)
    plt.plot(rx,ry, color="b")
    ######################################
    cl = ["r" if i**2 + j**2 <= R else "b" for i, j in zip(x,y)]
    plt.scatter(x, y, c=cl)
    plt.text(-1.2, 1.2, s = "R={R}, N={N}, error={error}, pi={pi}".format(R=R, N=N, error=error, pi=np.pi))
    plt.show()

if __name__ == "__main__":
    pi_exp(1, 10000)

运行结果如下图：

随着采样数的增大，预测的$\pi$与真实$\pi$之间的误差会越来越小。

参考：http://www.ruanyifeng.com/blog/2015/07/monte-carlo-method.html