蒙特卡罗（Monte Carlo）方法算积分

❝

蒙特卡罗（Monte Carlo）是摩纳哥最著名的一区，以豪华的赌场闻名于世，用它作为名字大概是因为随机性，就像赌博场里面的扔骰子的过程。最早的「蒙特卡罗方法」是为了解决一些难求解的积分问题。

❞

• 「问题」

• 「蒙特卡洛方法」

如果可以选择在的概率分布函数，则有：

若在之间是均匀分布时，即，那么：

这就是之前讲解的平均值法(点击跳转)，另外随机投点法(点击跳转)也是「蒙特卡洛方法」. 一般均匀分布并不是好选择，因为如果在有不少点使得，那么这些点对  的近似计算贡献很小，所以应尽可能少用这些点. 此时就需要采用「重要采样方法」选择合适的，从而提高精度，这部分内容我们后续会详细阐述，这次我们先分析「随机投点法」和「平均值法」的随机误差.

• 「误差分析」

（1）「随机投点法」

令  且 ，则 iid . 由中心极限定理知：

从而

所以

因此  的随机误差为：.

（2）「平均值法」

由中心极限定理知：

其中

因此  的随机误差为：，但其渐近方差更小.

类似的，计算高维定积分的蒙特卡罗方法的随机误差也为，所以蒙特卡罗方法计算积分和维数关系不大，但数值积分则存在「维数诅咒」问题，这也是蒙特卡罗方法的「优势」.

• 「高维积分算例」

「 以下为Python代码 」

import numpy as np
from scipy import integrate

##  (x1)^2(x2)^2(x3)^2 在 [0,1] 的积分

a1,b1 = 0,1
a2,b2 = 0,1
a3,b3 = 0,1

# 三重积分计算
def f(x1,x2,x3):
    return x1**2 * x2**2 * x3**2
    
I_exact, Error = integrate.tplquad(f,a1,b1,a2,b2,a3,b3)

# 平均值法
N = 10000
x1_sample = a1 + (b1-a1)*np.random.rand(N)
x2_sample = a2 + (b2-a2)*np.random.rand(N)
x3_sample = a3 + (b3-a3)*np.random.rand(N)
np.random.seed(1)

h_x = f(x1_sample,x2_sample,x3_sample)
I_approx_stat = (b3-a3)*(b2-a2)*(b1-a1)/N*np.sum(h_x)

# 数值积分
M = 200
h1 = (b1-a1)/(M-1)
h2 = (b2-a2)/(M-1)
h3 = (b3-a3)/(M-1)
x1 = np.linspace(a1,b1,M)
x2 = np.linspace(a2,b2,M)
x3 = np.linspace(a3,b3,M)
x1_mesh, x2_mesh, x3_mesh = np.meshgrid(x1,x2,x3)
I_approx_rec = np.sum( f(x1_mesh, x2_mesh, x3_mesh)*h1*h2*h3 )


print( '多重积分值：', I_exact )
print( '\n平均值法结果：', I_approx_stat )
print( '\n数值积分结果：', I_approx_rec )


❝

多重积分值：0.037037037037037035
平均值法结果：0.03737256369148107
数值积分结果：0.03788231093787493

（大家可尝试画出：不同数量采样点对应的结果和真实值之间的关系图）

❞