现代Monte Carlo算法的起源：Metropolis算法

Monte Carlo算法的思想起源很早，一个很著名的例子是18世纪的布丰投针实验。人们很早就发现Monte Carlo算法能够用于求数值积分。

统计物理中的Monte Carlo算法始于20世纪40年代，由波兰数学家斯塔尼斯拉夫·乌拉姆（Stanisław Marcin Ulam）提出，用来数值求解曼哈顿计划中中子散射问题的微分方程。1945年10月，第一台通用电子计算机ENIAC诞生，乌拉姆将这一思想告知冯·诺依曼（John von Neumann），随后冯·诺依曼将Monte Carlo算法编写成程序，利用计算机求解。而Monte Carlo这个名字，则由尼古拉斯·梅特罗波利斯（Nicholas Metropolis）提出，来源与乌拉姆的叔叔经常去赌博的摩纳哥蒙特卡罗赌场。

1953年，梅特罗波利斯与罗森布卢特（Rosenbluth）夫妇、泰勒（Teller）夫妇共同发表论文Equation of State Calculations by Fast Computing Machines[1]，提出了一种利用计算机求解强关联多体系统的统计性质的随机算法。在这篇文章中，他们计算了二维硬球（rigid sphere）模型。文章核心思路是将正则系综的相空间看做概率空间。如果给定N个粒子的坐标，那么可以求这一构型的能量（势能）。因为势能只与位置有关，与速度无关，于是相空间积分当中的对动量自由度积分便可约去，所以系统的统计性质只与势能部分有关

$\begin{equation} \la A\ra=\frac{\int A e^{-E/k _ BT}d^{dN}q}{\int e^{-E/k _ BT}d^{dN}q} \end{equation}$

上式中d表示维度，N表示粒子数量，q表示（广义）坐标，A为待求物理量，$\la A \ra$表示对对物理量A的统计平均。

一个直观的思路是随机给定粒子坐标，求这一构型的物理量A及能量E，然后给定其统计权重为$e^{-E/K_BT}$，最后求物理量A的加权平均。这一方法对稀疏系统可行，但是对稠密系统而言，由于绝大多数构型的权重极小（由于重叠而使得能量极大），而权重大的构型难以被采样，算法效率会极其低下。作者提出了修正方法，即举世闻名的Metropolis算法：与其随机均匀地抽取构型再赋以权重$e^{-E/k_BT}$，不如抽取概率为$e^{-E/k_BT}$的构型再按照平均。具体地，

随机地初始化N个粒子的坐标，求其能量。
随机选择一个粒子，将其移动$q_i\to q_i+\alpha\xi_i$，这里$\alpha$是步长，可以人为设定，$xi_i$是一个$[-1,1]范围内均匀分布的随机数。为了减小界面效应，作者采用了周期性边界条件。
求新构型与原构型的能量差（$\Delta E$）。如果新构型能量更低，则新构型被“接受”（Accepted）。如果新构型能量更高，那么取$[0,1]$均匀分布的随机数$\varepsilon$，若其小于$e^{-\Delta E/k_BT}$，则新构型被“接受”；否则被“拒绝”（rejected），即恢复到本次移动之前的构型。
【重要！】无论新构型是否被接受，经过2、3两步之后的构型被视为一个抽样构型。
重复步骤2-4多次，直到抽样构型足够多以控制统计误差到接受范围内。样本容量记为$M$。
对待求物理量做无权重的样本平均，平均值就是一次“试验”的结果。有时会需要求物理量的高阶矩，同样这个高阶矩是一次“试验”的结果。

以上就是Metropolis算法，其中第4步是其核心。显然，这是一个Markov过程，每一次转移（步骤2-4）得到一个状态，从第$n$个状态得到第$n+1$个状态的转移概率$P(X_{n+1}=x|X_{n}=x_n)$与$n$无关。

我们可以进行多次“试验”，得到多次“试验”的结果，再取平均得到物理量的计算结果；对多次“试验”结果求样本标准差得到物理量的计算误差。这是一个完整的计算过程。

对Metropolis算法采得的构型权重都正比于$e^{-E/k_BT}$，即

$\begin{equation} \frac{\nu_s}{\nu_r}=\frac{e^{-E_s/k_BT}}{e^{-E_r/k_BT}} \end{equation}$

对构型权重的证明

首先，一个粒子允许在整个空间中移动，只要步数足够多，那么这个粒子一定能够到达空间的任意位置。这就是说，算法满足遍历性（ergocity）。当然，这只是申明算法能够充分采样，并不需要在实际计算过程中遍历整个空间。

其次，需要证明如果样本容量足够大，样本池中状态为r的样本数量$\nu_r$应当正比于$e^{-E_r/k_BT}$。关于这点，首先明确这个过程中涉及到

先验（a priori）概率$P_{rs}$：当构型为s时，选择下一个构型为r的（条件）概率。显然$P_{rs}=P_{sr}$，因为粒子从位置1能够移动到位置2，那么单纯的看空间位置，当然也能从位置2移动到位置1，而且在步骤2中，选择移动位置是均匀随机的。
接受（acceptance）概率$\alpha_{rs}$：即从状态r到状态s能够被“接受”的概率。若$E_r>E_s$，那么$\alpha_{rs}=1$，$\alpha_{sr}=e^{-(E_r-E_s)/k_BT}$。

这样，Markov链的转移矩阵为

$\begin{equation} p_{r\to s}=P_{rs}\alpha_{rs}=P_{rs}\times\min[1,e^{-(E_s-E_r)/k_BT}], \qquad s\neq r \end{equation}$

由随机过程理论可知，当Markov过程达到平稳时，即采样充分时，满足细致平衡条件（detailed balance）

$\begin{equation} \nu_sp_{s\to r}=\nu_rp_{r\to s} \end{equation}$

即

$\begin{equation} \frac{\nu_s}{\nu_r}=\frac{p_{r\to s}}{p_{s\to r}}=\frac{P_{rs}\min[1,e^{-(E_s-E_r)/k_BT}]}{P_{sr}\min[1,e^{-(E_r-E_s)/k_BT}]}=\frac{\min[1,e^{-(E_s-E_r)/k_BT}]}{\min[1,e^{-(E_r-E_s)/k_BT}]}=\frac{e^{-E_s/k_BT}}{e^{-E_r/k_BT}} \end{equation}$

证毕。

这里特别说明：之所以步骤4非常重要，即必须在移动被拒绝的时候将原构型加入样本，可以通过计算$p_{r\to r}$看出：由于转移矩阵的性质$\sum_s p_{r\to s}=1$，

$\begin{align} p_{r\to r}=&1-\sum_{s\neq r}p_{r\to s}=1-\sum_{s\neq r}P_{rs}\min[1,e^{-(E_s-E_r)/k_BT}]\nonumber\newline =&1-\sum_{s\neq r|E_s\leqslant E_r}P_{rs}-\sum_{s\neq r|E_s>E_r}P_{rs}+\sum_{s\neq r|E_s>E_r}P_{rs}(1-e^{-(E_s-E_r)/k_BT})\nonumber\newline =&P_{rr}+\sum_{s\neq r|E_s>E_r}P_{rs}(1-e^{-(E_s-E_r)/k_BT}) \end{align}$

即：状态r回到自身，不仅可以通过直接在第2步选择自身并直接移回（$e^{-\Delta E/k_BT}=1$），还能够通过第2步选择其他状态但是被拒绝之后移回。

Metropolis算法的局限性

Metropolis算法能够保证在足够多的Monte Carlo采样下得到足够精确的物理量，但是也有重要的不足：

一般地，这个算法不能控制达到足够精度的步数，或者说不能告诉我们样本容量需要达到多大才能满足精度要求。达到足够的精度又称为“收敛”（converge）。Markov过程需要达到平稳才能保证抽样的正确分布，也就是说通常需要弃掉一开始的多个构型使得系统达到“平衡”，再进行采样。但是系统需要多久才能达到“平衡”并非已知。有的时候（例如在临界相变点附近），会发生临界降速（critical slowing-down），达到平衡的时间变为无穷长。
对于连续空间的Metropolis算法，步长$\alpha$取值尚未确定。显然，$\alpha$与收敛速度有关，如何调节$\alpha$值以使Monte Carlo算法最有效率，是一个很重要的问题。

经过维尔弗雷德·哈斯廷（Wilfred Keith Hastings）的总结推广[2]，Metropolis算法中接受概率可以不由波尔兹曼分布得出，而可以是任意要求的概率。这样推广的Metropolis算法成为一种通用的计算数学方法，称为Metropolis-Hastings算法。哈斯廷在总结中还讨论了Metropolis算法的误差估计等。

Metropolis算法是Markov Chain Monte Carlo（MCMC）算法的一种，也是重要性采样（Importance Sampling）方法的一例。除了Metropolis-Hastings算法外，由Metropolis算法还可以导出处理高维采样、计算联合分布的Gibbs抽样法。这些方法在计算统计学、机器学习等领域中有很重要的应用。关于这些算法，可以参考：
随机采样和随机模拟：吉布斯采样Gibbs Sampling
Gibbs采样

参考资料

Equation of State Calculations by Fast Computing Machines, Metropolis N., Rosenbluth A. W., Rosenbluth M. N., Teller A. H. and Teller E., Journal of Chemical Physics, 21, 1087-1092 (1953)
Hastings, W. Keith. “Monte Carlo sampling methods using Markov chains and their applications.” (1970): 97-109.