计算物理（四）Monte Carlo

本篇文章是西安交大《计算物理与程序设计》课程的笔记，程序语言为 MATLAB。主要内容包括具有特定分布随机数的生成，以及Monte Carlo积分法。

随机数

内置命令

• rand(n,m) 生成 $n$ 行 $m$ 列随机矩阵，矩阵元素为 $[0,1]$ 之间均匀分布的随机数。$a + (b-a)\cdot\mathrm{rand}$ 可生成 $[a,b]$ 之间均匀分布的随机数。
• randn(n,m) 生成 $n$ 行 $m$ 列随机矩阵，矩阵元素为 $[0,1]$ 之间标准正态分布的随机数。$\sigma \cdot \mathrm{rand} + \mu$ 可生成标准差为 $\sigma$，均值为 $\mu$ 的正态分布随机数。
• randi([a,b],n,m) 生成 $n$ 行 $m$ 列随机矩阵，矩阵元素为 $[a,b]$ 之间均匀分布的随机整数。
• exprnd(lambda,n,m) 生成 $n$ 行 $m$ 列随机矩阵，矩阵元素为服从均值为 $\lambda$ 的指数分布 $f(x)=\lambda \mathrm{e}^{-\lambda x}$ 的随机数。exprnd 需要 Statistics and Machine Learning 工具箱。

逆变换法生成随机数

定理1.1 设 $X$ 为连续型随机变量，取值于区间 $[a,b]$ ，其概率密度 $f(x)$ 在 $[a,b]$ 上取正值，$X$ 的分布函数为 $F(x) \in [0,1]$ 。若随机变量 $U$ 在 $[0,1]$ 上均匀分布，则 $Y = F^{-1}(U) \sim F(\cdot)$ ，即随机变量 $Y$ 与 $X$ 同分布。

定理1.2 设 $X$ 为离散型随机变量，取值于集合 $\left\{ a_1, a_2, a_3,\cdots \right\}$ 且 $a_i<a_{i+1}$，$F(x)$ 为 $X$ 的分布函数。$U$ 在 $[0,1]$ 上均匀分布。根据 $U$ 的值定义随机变量 $Y$ 为

$$F(a_{i-1}) < U \le F(a_i) \implies Y = a_i$$

则 $Y$ 与 $X$ 同分布。

舍选法生成随机数

舍选法 (acceptance-rejection method) 的原理是：要对一维随机变量 $X$ 进行采样，可以在二维笛卡尔坐标系中进行均匀随机抽样，并将样本保留在其概率密度函数 $f(x)$ 以下的区域中。舍选法可以扩展到 $n$ 维函数，只需在 $n+1$ 维坐标系中进行均匀随机抽样即可。

简单舍选法

假设需要生成 $[a,b]$ 区间上的随机数 $X$，其对应的概率密度函数为 $f(x)$ ，分布为 $F(x)$ 。简单舍选法的步骤如下：

1. 设 $f(x)$ 有一上界 $M$ ，即 $\forall x \in [a,b],\, f(x) \le M$ ；
2. 生成一个随机数 $y$ ，$y$ 来自对 $[a,b]$ 均匀分布 $U(a,b)$ 的抽样；
3. 生成一个随机数 $u$ ，$u$ 来自对 $[0,1]$ 均匀分布 $U(0,1)$ 的抽样；
4. 若 $u < \frac{f(y)}{M}$ ，接受 $x=y$ ，即将 $y$ 作为 $X$ 的一个样本，否则拒绝 $y$ 。
5. 回到步骤 2 并反复抽取 $(y,u)$，这样舍选得到的 $y$ 服从分布 $F(x)$ 。

证明：$(y,u)$ 的联合概率密度为

$$\varphi(y,u) = \frac{1}{b-a}\, ,(y,u)\in[a,b] \times [0,1]$$

按照简单舍选法抽出的随机数 $y$ 被成功留下且小于 $x$ 的概率为

$$P\left[y \le x\mid u \le \frac{f(y)}{M}\right] = \frac{P\left[y \le x , u \le \frac{f(y)}{M}\right]}{P\left[ u \le \frac{f(y)}{M}\right]} = \frac{\int_a^x \mathrm{d}y \int_0^{\frac{f(y)}{M}} \frac{1}{b-a} \mathrm{d}u}{\int_a^b \mathrm{d}y \int_0^{\frac{f(y)}{M}} \frac{1}{b-a} \mathrm{d}u} = M(b-a) \int_a^x \frac{f(y)}{M} \frac{1}{b-a} \mathrm{d}y = F(x)$$

$M$ 不一定需要是 $f(x)$ 的精确上界，但精确上界将使程序效率更高。 /* 注：$f(x)$ 有可能大于 $1$ ！不要跟分布 $F(x)$ 搞混了。*/

舍选法

假设需要生成随机数 $X$，其对应的概率密度函数为 $f(x)$ ，分布为 $F(x)$ 。舍选法的步骤如下：

1. 找一个可以方便生成随机数的随机变量 $Y$。假设概率密度 $g(y)$ 满足 $\frac{f(x)}{M g(x)} \le 1$ 对 $\forall x$ 都成立，$M$ 是某个常数；
2. 生成一个随机数 $y$ ，$y$ 来自对 $Y$ 的抽样；
3. 生成一个随机数 $u$ ，$u$ 来自对 $[0,1]$ 均匀分布 $U(0,1)$ 的抽样；
4. 若 $u < \frac{f(y)}{M g(y)}$ ，接受 $x=y$ ，即将 $y$ 作为 $X$ 的一个样本，否则拒绝 $y$ 。
5. 回到步骤 2 并反复抽取 $(y,u)$，这样舍选得到的 $x$ 服从分布 $F(x)$ 。

证明：$(y,u)$ 的联合概率密度为

$$\varphi(y,u) = g(y)$$

按照舍选法抽出的随机数 $y$ 被成功留下且小于 $x$ 的概率为

$$P\left[y \le x\mid u \le \frac{f(y)}{M}\right] = \frac{\int_{-\infty}^x \mathrm{d}y \int_0^{\frac{f(y)}{M g(y)}} g(y) \mathrm{d}u}{\int_{-\infty}^{+\infty} \mathrm{d}y \int_0^{\frac{f(y)}{M g(y)}} g(y) \mathrm{d}u} = M \int_{-\infty}^x \frac{f(y)}{M} \mathrm{d}y = \int_{-\infty}^xf(y) \mathrm{d}y = F(x)$$

/* 注：可以发现简单舍选法仅仅是取 $g(y) = 1$ 均匀分布的特例。 */

$M$ 不一定需要是 $f(x)/g(x)$ 的精确上界，但精确上界将使程序效率更高。

舍选法可以推广到高维，例如欲得到按照二维概率密度 $f(x,y)$ 分布的随机数，只需寻找与 $f(x,y)$ 相似且更易抽样的 $g(x,y)$ 即可。

Monte Carlo 积分

概述

Monte Carlo方法是通过抽取随机数的方式进行统计模拟的方法。Monte Carlo方法由大数定理保证

辛钦大数定理：设 $X_1,X_2,\cdots,X_N$ 是相互独立的随机变量序列，服从相同分布，且具有有限的数学期望 $\mu$ ，则对于 $\forall \varepsilon>0$ ，有

$$\lim _{N \rightarrow \infty} P\left\{\left|\frac{1}{N} \sum_{i=1}^N X_i-\mu\right|<\varepsilon\right\}=1$$

即样本均值依概率收敛于总体均值。

以及无意识统计学家定律 (Law of the unconscious statistician, LOTUS)

LOTUS：设随机变量 $X$ 的分布已知，其概率密度为 $f_X (x)$ 。则随机变量 $X$ 的函数 $g(X)$ 的期望为

$$\mathbb{E}[g(X)] = \int g(x) f_X (x) \, \mathrm{d}x ,\quad \text{for continuous } X$$

$$\mathbb{E}[g(X)] = \sum_x g(x) f_X (x) ,\quad \text{for discrete } X$$

均匀抽样法

设函数 $f(\boldsymbol{x})$ 在有限区域 $\Omega$ 上有界，区域 $\Omega$ 的体积为 $V$。取随机变量 $X \sim \mathrm{U}(\Omega)$ ，则由 LOTUS 定律，$f(X)$ 的期望为

$$\mathbb{E} f(X) = \int f(\boldsymbol{x}) \frac{1}{V} \mathrm{d}\boldsymbol{x} = \frac{1}{V} \int f(\boldsymbol{x}) \mathrm{d}\boldsymbol{x}$$

另一方面，定义随机变量 $Y = f(X)$ 。从 $X$ 中抽样得 $\boldsymbol{x}_1, \boldsymbol{x}_2, \cdots, \boldsymbol{x}_N$ ，进而得 $Y$ 的一个样本 $f(\boldsymbol{x}_1),f(\boldsymbol{x}_2),\cdots, f(\boldsymbol{x}_N)$ 。由于 $\boldsymbol{x}_i$ 独立同分布，故 $f(\boldsymbol{x}_i)$ 独立同分布。根据辛钦大数定律 $N \to \infty$ 时

$$\overline{Y} = \frac{1}{N} \sum_{i=1}^N f(\boldsymbol{x}_i) \to \mathbb{E}Y = \mathbb{E} f(X)$$

从而有

$$\frac{V}{N} \sum_{i=1}^N f(\boldsymbol{x}_i) = \int f(\boldsymbol{x}) \mathrm{d}\boldsymbol{x}$$

重要抽样法

均匀抽样法和重要抽样法法的区别在于重要抽样法中 $X$ 不再在 $\Omega$ 上均匀分布，而是有恰当的概率密度。

设函数 $f(\boldsymbol{x})$ 在有限区域 $\Omega$ 上有界，区域 $\Omega$ 的体积为 $V$。取概率密度为 $\rho(\Omega)$ 的随机变量 $X$ ，满足 $\int_\Omega \rho(\boldsymbol{x}) \mathrm{d}\boldsymbol{x} = 1$ 。定义随机变量 $X$ 的函数

$$g(X) = \frac{f(X)}{\rho (X)}$$

则由 LOTUS 定律

$$\mathbb{E} g(X) = \int \frac{f(\boldsymbol{x})}{\rho(\boldsymbol{x})} \rho(\boldsymbol{x}) \, \mathrm{d}\boldsymbol{x} = \int f(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x}$$

另一方面由大数定律

$$\overline{g(X)} = \frac{1}{N} \sum_{i=1}^N g(\boldsymbol{x}_i) \to \mathbb{E} g(X)$$

因此有

$$\int f(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x} = \lim_{N \to \infty} \frac{1}{N} \sum_{i=1}^N g(\boldsymbol{x}_i)$$

重要抽样法与均匀抽样法相比，其优势在于可以恰当地选取 $\rho(x)$ ，使得 $\frac{f(x)}{\rho(x)}$ 在 $\Omega$ 上起伏不大，保证 $\mathbb{D} (f/\rho) < \mathbb{D} f$。从而用样本均值 $\frac{1}{N} \sum_{i=1}^N g(\boldsymbol{x}_i)$ 替代总体均值的误差甚小。

Proof

例：计算 $\int_0^\infty \frac{3+4x}{1+x^2} \mathrm{e}^{-x} \, \mathrm{d}x$ 。

解：在 $(0,\infty)$ 上，$\rho(x) = \mathrm{e}^{-x}$ 是合理的概率密度，有 $\int_0^\infty \mathrm{e}^{-x} \, \mathrm{d} x=1$ ，以及分布 $\Rho(x) = 1-\mathrm{e}^{-x}$ 。假设能够生成以 $\mathrm{e}^{-x}$ 为概率密度的随机变量 $X$，则

$$\int_0^\infty \frac{3+4x}{1+x^2} \mathrm{e}^{-x} \, \mathrm{d}x = \lim_{N \to \infty} \frac{1}{N} \sum_{i=1}^N \frac{3+4x_i}{1+x_i^2}$$

其中 $x_i \sim \Rho(x)$ 。问题在于服从 $\Rho(x)$ 分布的随机变量如何生成。利用定理1.1，可以先生成 $U \sim \mathrm{U}(0,1)$ ，然后通过逆变换

$$x_i = \Rho^{-1}(u_i) = - \ln (1-u_i)$$

得到 $x_i$ 。事实上，如果事先把待求积分变形

$$\int_0^\infty \frac{3+4x}{1+x^2} \mathrm{e}^{-x} \, \mathrm{d}x = \int_0^1 \frac{3+4x}{1+x^2} \, \mathrm{d} (1-\mathrm{e}^{-x}) = \int_0^1 \frac{3-4\ln(1-u)}{1+\ln^2(1-u)} \, \mathrm{d} u$$

然后再用 $[0,1]$ 上的均匀抽样法，得到的是相同的结果。

Metropolis算法

Metropolis-Hastings 算法是一种马尔科夫链蒙特卡洛方法 (Markov Chain Monte Carlo, MCMC)，用于在难以直接采样时从某一概率分布中抽取随机样本序列，得到的序列可用于估计该概率分布。Markov 链是指 $t+1$ 步的随机变量在给定第 $t$ 步随机变量后与其余的随机变量条件独立(conditionally independent) ，用概率论的语言书写就是

$$P\left(X_{t+1} \mid X_t,X_{t-1},\cdots,X_1 \right) = P \left(X_{t+1} \mid X_t \right)$$

随机游走 (random walk) 便是一个 Markov 过程，游走者每次都从当前位置出发以固定概率随机向某个方向移动一个单位。Metropolis-Hastings 算法正是利用这种方法抽取随机样本序列的。

设 $f(x)$ 为目标概率密度，设想一个随机游走者位于 $x^t$ 的位置，$q(x|x^t)$ 为建议分布 (proposal distribution) ，定义为当前位置为 $x^t$ 时下一个位置 $x$ 的条件概率密度分布。Metropolis-Hastings 算法为

1. 从建议分布 $q(x|x^t)$ 中随机抽取候选状态 $x^*$
2. 计算是否采纳候选状态的概率 $r = \min \left\{ 1, \frac{f(x^*)}{f(x^t)}\frac{q(x^t|x^*)}{q(x^*|x^t)} \right\}$
3. 从 $\mathrm{U}(0,1)$ 中随机抽取 $u$ 。若 $u\le r$，接受该状态，并令 $x^{t+1}=x^*$ ；若 $u>r$，拒绝该状态，并令 $x^{t+1} = x^t$ 。
4. 回到步骤1反复迭代

当随机游走步数 $N$ 足够多时，$\{x^t\}$ 将趋近于概率密度为 $f(x)$ 的分布。

以下是一些说明：

• 建议分布 (proposal distribution) 为游走者提供下一步可选的状态。最常见的有均匀分布，例如规定一维随机游走者每次可等概率行进到距离当前位置 $\pm \delta$ 范围内的任意位置，则 $q(x|x^t) = \frac{1}{2\delta} , x \in [x^t -\delta,x^t+\delta]$ ；还有正态分布，例如规定一维随机游走者每次可按照以当前位置为中心，标准差 $\sigma$ 的正态分布行进到下一位置，则 $q(x|x^t) = \frac{1}{\sqrt{2 \pi} \sigma} \exp \left[-\frac{(x-x^t)^2}{2 \sigma^2}\right]$ 。这两种常见的分布都是对称的 (symmetric)，即 $q(x|x^*) = q(x^*|x)$，这使得 Metropolis-Hastings 算法的第二步概率公式化简为 $r = \min \left\{ 1, \frac{f(x^*)}{f(x^t)}\right\}$，这是 Metropolis 最初提出 Metropolis 算法时的原始形式。

• Metropolis-Hastings 算法的优点在于，目标概率密度 $f(x)$ 不必是归一化的，因为算法只需要用到两个状态之间的相对概率密度 $\frac{f(x^*)}{f(x^t)}$ 。这在统计物理中极为关键，以经典统计中的正则系综为例，能态 $E$ 出现的概率密度为

  $$\rho(E) = \frac{1}{Z} \mathrm{e}^{-\beta E}$$

  其中 $Z$ 为配分函数。但问题在于配分函数往往是不易求的，因为对态空间的积分十分困难，而且求出配分函数了还要你蒙特卡洛干什么？直接利用热力学公式便可获知体系的性质了。Metropolis-Hastings 算法提供了在不知道配分函数 $Z$ 的情况下，仅通过不同能级出现概率的比值 $\frac{\mathrm{e}^{-\beta E_1}}{\mathrm{e}^{-\beta E_2}}$，便可模拟系统的性质。