蒙特卡罗积分

在数学中，蒙特卡罗积分（Monte Carlo integration）是一种使用随机数进行数值积分的技术。它是一种特殊的蒙特卡罗方法，可对定积分进行数值计算。其他算法通常在规则网格上评估被积函数，而^[1]蒙特卡洛随机选择被积函数评估的点。 ^[2]该方法对于高维积分特别有用。 ^[3]

蒙特卡罗积分示意图。在该示例中，D是和正方形内切的圆，域E则是正方形。因为正方形的面积 (值是4) 很容易计算，所以圆的面积 (π*1.0²) 可以通过圆内的点 (40个红点) 与总点数 (50个点) 的比率(0.8) 来估计，得出圆面积的近似值 4*0.8 = 3.2 ≈ π。

进行蒙特卡罗积分的方法有多种，例如均匀采样、分层采样、重要性采样、顺序蒙特卡罗（也称为粒子滤波器）和平均场粒子方法。

概述

在数值积分中，梯形法则等方法使用确定性方法。另一方面，蒙特卡罗积分采用非确定性方法：每种实现都提供不同的结果。在蒙特卡罗中，最终结果是带有各自误差线的正确值的近似值，并且正确值很可能在这些误差线内。

考虑函数$${\displaystyle I=\int _{\Omega }f({\overline {\mathbf {x} }})\,d{\overline {\mathbf {x} }}}$$其中 Ω 是${\displaystyle \mathbb {R} ^{m}}$的一个子集，它的体积是：$${\displaystyle V=\int _{\Omega }d{\overline {\mathbf {x} }}}$$

朴素的蒙特卡罗方法是在Ω上均匀采样点： ^[4] 给定N个均匀样本，$${\displaystyle {\overline {\mathbf {x} }}_{1},\cdots ,{\overline {\mathbf {x} }}_{N}\in \Omega ,}$$I可以被近似为：$${\displaystyle I\approx Q_{N}\equiv V{\frac {1}{N}}\sum _{i=1}^{N}f({\overline {\mathbf {x} }}_{i})=V\langle f\rangle .}$$ 这是因为大数定律确保以下结论成立：$${\displaystyle \lim _{N\to \infty }Q_{N}=I.}$$ 我们使用I给出Q_N的I估计， Q_N的误差线（error bars）可以使用方差的无偏估计通过样本方差来估计。$${\displaystyle \mathrm {Var} (f)\equiv \sigma _{N}^{2}={\frac {1}{N-1}}\sum _{i=1}^{N}\left(f({\overline {\mathbf {x} }}_{i})-\langle f\rangle \right)^{2}.}$$由此可得：$${\displaystyle \mathrm {Var} (Q_{N})={\frac {V^{2}}{N^{2}}}\sum _{i=1}^{N}\mathrm {Var} (f)=V^{2}{\frac {\mathrm {Var} (f)}{N}}=V^{2}{\frac {\sigma _{N}^{2}}{N}}.}$$只要以下序列$${\displaystyle \left\{\sigma _{1}^{2},\sigma _{2}^{2},\sigma _{3}^{2},\ldots \right\}}$$是有界的，该方差就会逐渐减小到零，即 1/N 。 那么，Q_N的误差估计就是：$${\displaystyle \delta Q_{N}\approx {\sqrt {\mathrm {Var} (Q_{N})}}=V{\frac {\sigma _{N}}{\sqrt {N}}},}$$

我们看到，Q_N的误差是自身的${\displaystyle {\tfrac {1}{\sqrt {N}}}}$，这是平均值的标准误差乘以${\displaystyle V}$。该结果不依赖于积分的维数，这是蒙特卡洛积分相对于大多数指数依赖维数的确定性方法所具有的优势。 值得注意的是，与确定性方法不同，误差的估计不是严格的误差界限；随机抽样可能无法揭示被积函数的所有重要特征，从而导致误差的低估。

虽然朴素蒙特卡罗适用于简单的示例，但对确定性算法的改进只能通过使用特定于问题的采样分布的算法来实现。通过适当的样本分布，可以利用以下事实：几乎所有高维被积函数都是非常局部化的，并且只有小子空间对积分有显着贡献。 ^[5]蒙特卡罗文献的很大一部分致力于开发改进误差估计的策略。特别是，分层采样（将区域划分为子域）和重要性采样（从非均匀分布中采样）是此类技术的两个示例。

例子

使用蒙特卡洛方法最典型的例子是估计π。

我们考虑以下函数：$${\displaystyle H\left(x,y\right)={\begin{cases}1&{\text{if }}x^{2}+y^{2}\leq 1\\0&{\text{else}}\end{cases}}}$$即正方形集合Ω=[−1,1]×[−1,1]的体积是V = 4。 因此，使用蒙特卡罗积分计算π值的粗略方法是在Ω上选取N个随机数并计算：$${\displaystyle Q_{N}=4{\frac {1}{N}}\sum _{i=1}^{N}H(x_{i},y_{i})}$$右图中，相对误差${\displaystyle {\tfrac {Q_{N}-\pi }{\pi }}}$被认为是一个关于N的函数，由此确认${\displaystyle {\tfrac {1}{\sqrt {N}}}}$ 。

C 示例

请记住，以下的程序需要真正的随机数生成器。

int i, throws = 99999, insideCircle = 0;
double randX, randY, pi;

srand(time(NULL));

for (i = 0; i < throws; ++i) {
  randX = rand() / (double) RAND_MAX;
  randY = rand() / (double) RAND_MAX;
  if (randX * randX + randY * randY < 1) ++insideCircle;
}

pi = 4.0 * insideCircle / throws;

Python示例

from numpy import random
import numpy as np

throws = 2000
inside_circle = 0
i = 0
radius = 1
while i < throws:
    # Choose random X and Y centered around 0,0
    x = random.uniform(-radius, radius)
    y = random.uniform(-radius, radius)
    # If the point is inside circle, increase variable
    if x**2 + y**2 <= radius**2:
        inside_circle += 1
    i += 1

# Calculate area and print; should be closer to Pi with increasing number of throws
area = (((2 * radius) ** 2) * inside_circle) / throws
print(area)

参见

• 拟蒙特卡罗方法
• 蒙特卡洛辅助场
• 统计物理中的蒙特卡罗方法
• 蒙特卡罗方法
• 减小方差（英语：Variance reduction）