蒙特卡洛方法——玩具例子(2)定积分

蒙特卡洛方法可以计算任意一个积分的值。

下面以求$\int_{0}^{1} x^2 dx$的定积分，定义域为[0,1], $x^2$的值域为[0,1]。 在坐标区域$0 \leq x \leq 1$; $0 \leq y \leq 1$ 内撒点，用纵坐标j小于等于$i^2$的所有点$(i,j)$ 占总点数的比例近似估计定积分的值。

# -*- coding:utf8 -*-

import numpy as np
import matplotlib.pyplot as plt

def interplote_exp(up, down, N):
    x = np.random.uniform(up, down, (N,))
    y = np.random.uniform(up, down, (N,))
    incircle = np.sum(y - np.power(x, 2) <= 0.0)
    error = 1.0/3 - incircle * 1.0/N
    ######################################
    fx = np.linspace(0.0, 1.2, 3600)
    fy = np.power(fx, 2)
    ######################################
    plt.fill_betweenx(fy, 1.0, fx, where=fx <= 1.0, facecolor='red')
    rx = np.asarray([0.0, 1.0, 1.0, 0.0, 0.0], dtype = np.float)
    ry = np.asarray([0.0, 0.0, 1.0, 1.0, 0.0], dtype = np.float)
    plt.plot(rx,ry, color="g")
    plt.plot([0.0, 1.2],[0.0, 0.0], color="b")
    plt.plot([0.0, 0.0],[0.0, 1.2], color="b")
    plt.plot(fx,fy, color="b")
    cl = ["b" if j - i**2 <= 0.0 else "r" for i, j in zip(x,y)]
    plt.scatter(x, y, c=cl)
    plt.text(0.2, 1.4, s = r"$\int_{0}^{1} x^2 dx=\frac{1}{3}$") 
    plt.text(0.2, 1.2, s = "error={error}, N={N}".format(N=N, error=error))
    plt.show()

if __name__ == "__main__":
    interplote_exp(0.0, 1.0, 10000)
    pass

输出的图示为：

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