蒙特卡洛方法

蒙特卡洛模拟是一个广义概念，指的是使用随机数和统计抽样来解决问题的方法，而非精确建模。从这种抽样的性质来看，结果会有一些不确定性，但统计学上的 “大数定理 “会确保不确定性随着样本数量的增加而下降。

统计抽样的一个重要工具是随机数发生器。参见教程32，了解随机数生成。

激励

在边长为2的正方形中随机选取一个点，它落入正方形内切圆中的概率为$\pi/4$。因此我们可以通过大量随机选取点来估计$\pi$的值。我们甚至可以做一个简单的物理实验：假设院子中存在一个形状不规则的池塘，而院子本身是面积已知的长方形。如果现在向院子里扔卵石，使它们在任何给定的地方都有同样的可能性，那么落在池塘里的卵石和落在外面的卵石的比例就等于面积的比例。

从数学角度讲，令$\Omega\in [0，1]^2$为某个区域，函数$f(\tilde{x})$为$\Omega$的边界，即 $$ \left{\begin{array}{ll} f(\bar{x})<0 & x \notin \Omega \\ f(\bar{x})>0 & x \in \Omega \end{array}\right. $$ 现在取随机点$\tilde{x}_0,\tilde{x}_1,\tilde{x}_2在[0,1]^2$，我们可以通过计算$f(\bar{x}_i)$为正或负的频率来估计Ω的面积。

我们可以把这个想法扩展到积分。一个函数在一个区间$(a, b)$上的平均值定义为 $$ \langle f\rangle=\frac{1}{b-a} \int{a}^{b} f(x) d x $$ 另一方面，我们也可以估计其平均数为 $$ \langle f\rangle \approx \frac{1}{N} \sum{i=1}^{n} f\left(x{i}\right) $$ 如果点$x$的分布合理，并且函数$f$不是太苛刻的话。可以得出 $$ \int{a}^{b} f(x) d x \approx(b-a) \frac{1}{N} \sum{i=1}^{n} f\left(x{i}\right) $$ 统计理论告诉我们，积分中的不确定性$\sigmaI$与下列因素有关。标准差$\sigma_f$的表达为 $$ \sigma{I} \sim \frac{1}{\sqrt{N}} \sigma_{f} $$ 且服从正态分布。

吸引人的地方是什么？

到目前为止，蒙特卡洛积分看起来与经典的黎曼和积分没有什么区别，但当增加空间维度时，差别就较为明显。此时，经典积分需要在每个维度上有$N$个点，因此有$d$维上有$N^d$个点。另一方面，在蒙特卡洛方法中，点是从$d$维空间中随机抽取的，而且点的数量少得多。

在计算方面，蒙特卡洛方法很有吸引力，因为所有的函数评估都可以并行进行。

对标准差为$\sigma$的事件进行$N$次取样观察，其平均标准差为：$\sigma / \sqrt{N}$。这意味着更多采样观测将降低平均标准差以提高精度；使蒙特卡洛方法的有趣之处在于，这种准确性的提高与原始问题的维度无关。

蒙特卡洛方法是模拟具有统计性质现象的不二选择，如放射性衰变或布朗运动。而蒙特卡洛方法其他应用范畴则不属于科学计算领域，例如，用于股票期权定价的Black-Scholes模型[15]就使用了蒙特卡洛模拟。

线性方程组也可以通过蒙特卡洛方法求解，但这并不常见。下面我们将以一个用精确方法解答十分复杂的应用为例，说明蒙特卡罗方法所代表的统计抽样手段的便捷之处。

示例

伊辛模型的蒙特卡洛模拟

伊辛模型（介绍见[33]）最初是为了模拟铁磁性而提出的。磁性是原子排列其 “自旋 “方向的结果：假设自旋只能是 “向上 “或 “向下”，那么如果有更多的原子自旋向上，而不是向下，那么这种材料就具有磁性，反之亦然。这些原子被称为 “晶格 “结构中的原子。

现在想象一下，加热一种材料会使原子松动。如果在材料上施加一个外部场，原子将开始与场对齐，如果场被移除，磁性又会消失。然而在某个临界温度以下，材料将保留其磁性。我们将使用蒙特卡洛模拟来寻找保留的稳定构型。

假设晶格$\Lambda$有$N$个原子，我们把原子的构型表示为$\sigma=(\sigma1，…，\sigma_N)$，其中每个$\sigma_i=\pm1$。 晶格的能量模型为 $$ H=H(\sigma)=-J \sum{i} \sigma{i}-E \sum{i j} \sigma{i} \sigma{j} $$ 第一个项模拟单个自旋$\sigma_i$与强度为$J$的外部场的相互作用。第二项是对近邻对的求和，模拟原子对的排列：如果原子具有相同的自旋，则乘积$\sigma_i\sigma_j$为正，如果相反则为负。

在统计力学中，构型的概率是 $$ P(\sigma)=\exp (-H(\sigma)) / Z $$ 其中 “分区函数 “$Z$定义为 $$ Z=\sum_{\sigma} \exp (H(\sigma)) $$ 其中，总和在所有$2^N$配置上运行。

如果一个构型的能量在微小扰动下没有减少，那么它就是稳定的。为了探索这一点，我们在晶格上进行迭代，探索改变原子的自旋是否会降低能量。我们引入了一个偶然因素来防止人为的解决方案。(这就是Metropolis算法[155]）。

for fixed number of iterations do for each atom 𝑖 do
calculate the change Δ𝐸 from changing the sign of 𝜎𝑖
if Δ𝐸 < or exp(−Δ𝐸) greater than some random number then
accept the change

如果我们注意到与稀疏矩阵-向量乘积结构的相似性，这种算法可以被并行化。在该算法中，我们也是通过结合几个近邻的输入来计算一个局部数量。这意味着我们可以对晶格进行分区，在每个处理器收集到一个幽灵区域（ghost region）后计算局部更新。

让每个处理器在晶格中的局部点上迭代，相当于晶格的一个特定的全局排序；为了使并行计算等同于串行计算，我们还需要一个并行随机发生器（第32.3节）。