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

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

下面以求∫10x2dx的定积分，定义域为[0,1], x2的值域为[0,1]。 在坐标区域0≤x≤1; 0≤y≤1 内撒点，用纵坐标j小于等于i2的所有点(i,j) 占总点数的比例近似估计定积分的值。

# -*- coding:utf8 -*-import numpy as npimport matplotlib.pyplot as pltdef 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