Statistics Simulation (Monte Carlo Sampling)

实例：生成具有如下概率分布的随机变量 $X$ :
$$p_1=0.20,p_2=0.15,p_3=0.25,p_4=0.40\text{where}{p}_i=P\{X=i\}$$
用逆变换抽样的算法如下：
首先产生一个 $(0,1)$ 上的随机数 $U$ ，然后：
如果 $U<0.10$ 则令 $X=1$ 并停止；
如果 $U<0.35$ 则令 $X=2$ 并停止；
如果 $U<0.60$ 则令 $X=3$ 并停止；
其他情况，令 $X=4$ 。

实例：生成具有如下概率密度函数的随机变量 $X$ :
$$f(x)=\{\begin{array}{ll}8x& \text{if}0\le x<0.25\\ \frac{8}{3}-\frac{8}{3}x& \text{if}0.25\le x<1\\ 0& \text{otherwise}\end{array}$$

用逆变换抽样(Inverse transform sampling)的方法如下：
首先，由概率密度函数求得对应的分布函数：
$$\begin{array}{rl}F(x)& ={\int }_{-\mathrm{\infty }}^{x}f(t)\mathrm{d}x\\ & =\{\begin{array}{ll}0& \text{if}x<0\\ 4x^2& \text{if}0\le x<0.25\\ \frac{8}{3}x-\frac{4}{3}{x}^2-\frac{1}{3}& \text{if}0.25\le x\le 1\\ 1& \text{if}x>1\end{array}\end{array}$$

然后，生成 $(0,1)$ 上的随机数 $U$ ，利用上面分布函数的反函数可得 $X$ 的一个样本值：
$$X=F^{-1}(U)=\{\begin{array}{ll}\frac{\sqrt{U}}{2}& \text{if}0\le U<0.25\\ 1-\frac{\sqrt{3(1-U)}}{2}& \text{if}0.25\le U\le 1\end{array}$$

摘自：
Sampling Methods, by Patrick Lam: http://www.people.fas.harvard.edu/~plam/teaching/methods/sampling/sampling.pdf

随机变量 $X$ 的概率密度如下：
$$f(x)=20x(1-x{)}^3,0<x<1$$
请用拒绝抽样方法来生成它的样本。

我们考虑一个更简单的随机变量： $(0,1)$ 上的随机数，即满足：
$$g(x)=1,0<x<1$$

然后，寻找一个最小的常数 $c$ 满足 $f(x)/g(x)\le c$ ，这可以通过微积分对函数求极值的知识来求得，即求下式的最大值：
$$\frac{f(x)}{g(x)}=20x(1-x{)}^3$$
用求导法，得上式在 $x=\frac{1}{4}$ 时取得最大值，这时有：
$$\frac{f(x)}{g(x)}\le 20\left(\frac{1}{4}\right){\left(\frac{3}{4}\right)}^3=\frac{135}{64}\equiv c$$
从而：
$$\frac{f(x)}{cg(x)}=\frac{256}{27}x(1-x{)}^3$$

对于这个问题，拒绝抽样的步骤如下：
第 1 步，生成两个 $(0,1)$ 上的随机数 $U_1,U_2$
第 2 步，如果 $U_2\le \frac{256}{27}{U}_1(1-U_1{)}^3$ 成立，则 $X=U_1$ ；如果不成立，则回到第 1 步接着执行。

说明：这个问题中，产生一个随机变量样本值的上面算法的平均迭代次数为 $c=\frac{135}{64}\approx 2.11$

实例：请生成下面 Gamma 分布的样本：
$$f(x)=K{x}^{1/2}{e}^{-x},x>0$$
上式中 $K=1/\mathrm{\Gamma }(\frac{3}{2})=2/\sqrt{\pi }$

Gamma 分布的分布函数不便于直接求反函数，故不方便直接用逆变换抽样。可以用拒绝抽样法，思路是先产生 Exponential 分布（可以直接用逆变换抽样可对指数分布进行抽样）的样本。

详细过程略。可参考 Simulation, Sheldon M Ross, 5th, 5.2 The Rejection Method

如果使用的 $q(x)$ 和 $c$ 不合适（概率密度和 $p(x)$ 不相近），则被拒绝的点会很多，从而产生样本值的效率会很低。

Adaptive rejection sampling 可以克服上述的拒绝抽样缺点，具体介绍略。

假设问题是：模拟具有概率密度函数为 $f(x)$ 的随机变量 $X$ 。

SIR 算法步骤如下：
(1) 寻找另外一个容易抽样的概率函数为 $g(x)$ 的随机变量 $Y$ ，模拟随机变量 $Y$ 产生 $N$ 个样本值，并记为：
$$\{y_1,y_2,\cdots ,y_N\}$$
(2) 定义权重 $w_i,i=1,2,\cdots ,N$ 为：
$$w_i=\frac{f(y_i)}{g(y_i)}$$
(3) 模拟随机变量 $X$ ，使其满足：
$$P\{X=y_j\}=\frac{w_j}{\sum _{i=1}^N{w}_i},j=1,2,\cdots ,N$$
那么，当 $N$ 充分大时，按上面方法产生的随机变量 $X$ 具有概率函数 $f$ （证明略）。

说明 1：上面过程中，步骤(3)称为 “重抽样(Resample)” ，因为它是在步骤(1)产生的样本集合中再次进行抽样。
说明 2：在步骤(2)中计算 $w_i=\frac{f(y_i)}{g(y_i)}$ 时， $f(y_i)$ 和 $g(y_i)$ 中的常数因子可以省略。假设 $f(y_i)=C_1{f}_0(y_i),g(y_i)=C_2{g}_0(y_i)$ ，其中 $C_1,C_2$ 为常数，则 $w_i$ 可以这样计算 $w_i=\frac{f_0(y_i)}{g_0(y_i)}$ ，这是因为常数因子会同时出现在步骤(3)的分子和分母中，不会影响步骤(3)的计算结果。所以，我们也可以把 $w_i$ 记为：
$$\begin{array}{rl}{w}_i& =\eta \frac{f(y_i)}{g(y_i)}\\ & \propto \frac{f(y_i)}{g(y_i)}\end{array}$$
说明 3：上面描述中，一般 $f(x)$ 称为 target 分布，$g(x)$ 一般称为 proposal 分布。

下面是用 R 语言实现 SIR 算法的一个例子。

假设需要模拟的随机变量概率密度为： $f(x)=\frac{1}{c}h(x)$ ，其中 $c$ 为未知的常数， $h(x)=e^{0.4\ast (x-0.4{)}^2-0.08\ast x^4}$ 。尽管 $f(x)$ 中有未知常数 $c$ ，我们也可以用 SIR 算法对 $f(x)$ 进行抽样。

# SIR example
k <- function(x, a=.4, b=.08){exp(a*(x-a)^2 - b*x^4)}  # k(x)本身不是概率密度函数，它是概率密度函数f(x)乘以某个未知常数
x <- seq(-4, 4, 0.1)

g <- function(x){rep.int(0.125,length(x))}   # g(x)为proposal分布，这里是(-4,4)上的均匀分布
N <- 1000
x.g <- runif( N, -4, 4 )       # 第1步，对g(x)进行抽样
ww <- k(x.g) / g(x.g)          # 第2步，计算每个样本的权重
qq <- ww / sum(ww)
x.acc <- sample(x.g, prob=qq, replace=T)  # # 第3步，resample

pdf(file="myplot_SIR.pdf")
par(mfrow=c(2,2), mar=c(2,2,1,1))
plot(x,k(x),type="l")
plot(x,g(x),type="l")
barplot(table(round(x.acc,1))/length(x.acc))

上面程序生成的 myplot_SIR.pdf 如图 6 所示。

注：这个例子比较简单，也可以用拒绝抽样的方法。

参考：http://people.math.aau.dk/~kkb/Undervisning/Bayes14/sorenh/docs/sampling-notes.pdf

用“样本均值”来估计“总体均值”时，样本值的数量当然是越多越好。但产生更多的样本值意味更程序运行更长的时间，何时可以停止产生新数据呢？

结合前面的分析，我们可以采用下面的方法来确定何时可以停止产生新数据。
第 1 步，选择一个合适的值 $d$ 作为估计量的标准误差。
第 2 步，产生至少 100 个数据。
第 3 步，继续产生新数据，直到 $k$ 个数据时 $\frac{S_k}{\sqrt{k}}<d$ 成立就可以停止，其中 $S_k$ 是这 $k$ 个数据确定的样本标准差（即 $S_k=\sqrt{S_k^2}=\sqrt{\frac{1}{k-1}\sum _{i=1}^k(X_i-\bar{X}{)}^2}$ ）。
第 4 步，总体均值 $\theta$ 的估计量就是样本均值 $\bar{X}=\sum _{i=1}^k{X}_i/k$ 。

说明 1：在第 3 步中，每次产生一个新数据，都要计算 $S_k$ 并验算 $\frac{S_k}{\sqrt{k}}<d$ 是否成立。在编程实现时为了更快速度，可以采用迭代的方式来计算 $S_k$ ，这样可以利用 $S_j$ 快速计算 $S_{j+1}$ ，具体迭代公式（证明略）如下：
$$\begin{array}{c}{\bar{X}}_0=0,S_1^2=0\\ {\bar{X}}_{j+1}={\bar{X}}_j+\frac{X_{j+1}-{\bar{X}}_j}{j+1}\\ S_{j+1}^2=(1-\frac{1}{j})S_j^2+(j+1){({\bar{X}}_{j+1}-{\bar{X}}_j)}^2\end{array}$$

说明 2：如果随机变量的方差是期望的函数（如Bernoulli distribution），这时可以调整上面算法。因为在推导上面算法时用样本方差 $S^2$ 来估计总体方差 ${\sigma }^2$ ，而对于 Bernoulli distribution 来说没有必要用样本方差估计总体方差。因为：
$$X_i=\{\begin{array}{l}1\text{with probability p}\\ 0\text{with probability 1-p}\end{array}$$
Bernoulli distribution 的均值为 $p$ ，方差为 $p(1-p)$ 。
我们直接用 ${\bar{X}}_n(1-{\bar{X}}_n)$ 来估计总体方差 ${\sigma }^2$ 即可。

所以，对于估计 Bernoulli distribution 的期望，其算法为：
第 1 步，选择一个合适的值 $d$ 作为估计量的标准误差。
第 2 步，产生至少 100 个数据。
第 3 步，继续产生新数据，直到 $k$ 个数据时 $\sqrt{\frac{{\bar{X}}_k(1-{\bar{X}}_k)}{k}}<d$ 成立就可以停止。
第 4 步，Bernoulli distribution 的期望为样本均值 ${\bar{X}}_k=\sum _{i=1}^k{X}_i/k$ 。

估计 $\theta ={\int }_0^1{e}^x\mathrm{d}x$ 。

容易知道 $E[e^U]={\int }_0^1{e}^x\mathrm{d}x$ （其中 $U\sim U(0,1)$ ），即求随机变量 $e^U$ 的均值即可估计 $\theta$ 。
说明：容易精确计算出结果为 $\theta =e-1$ ，在这里我们仅仅是通过该例子来说明对偶变量法可以提高模拟精度。

(1)如果产生相互独立的两个随机数 $U_1,U_2$ ，则样本均值的方差（也就是估计的均方误差）为：
$$Var\left(\frac{e^{U_1}+e^{U_2}}{2}\right)=\frac{Var(e^U)}{2}$$
其中， $Var(e^U)=E[e^{2U}]-(E[e^U]{)}^2={\int }_0^1{e}^{2x}\mathrm{d}x-(e-1{)}^2=0.2420$ ，所以：
$$Var\left(\frac{e^{U_1}+e^{U_2}}{2}\right)=\frac{0.2420}{2}=0.1210$$

(2)如果使用对偶变量 $U,1-U$ ，容易验证 $Cov(e^U,e^{1-U})=E[e^U{e}^{1-U}]-E[e^U]E[e^{1-U}]=e-(e-1{)}^2=-0.2342$ ，则样本均值的方差（也就是估计的均方误差）为：
$$Var\left(\frac{e^U+e^{1-U}}{2}\right)=\frac{Var(e^U)}{2}+\frac{Cov(e^U,e^{1-U})}{2}=0.0039$$

由此可知，使用对偶变量时，样本均值的方差减少了： $(0.1210-0.0039)/0.1210=96.7\mathrm{\%}$

先回顾一下前面介绍的用模拟方法计算积分的知识。比如，要计算下面积分：
$$\theta ={\int }_0^{1}g(x)\mathrm{d}x$$
解决办法：假设 $U$ 是 $(0,1)$ 上的均匀分布，即它的概率分布为：
$$f(x)=\{\begin{array}{l}1,0<x<1\\ 0,\text{otherwise}\end{array}$$
则：
$$E[g(U)]={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}f(x)g(x)\mathrm{d}x={\int }_0^{1}g(x)\mathrm{d}x$$

从而我们把积分的计算转换成为求随机变量函数（仍是随机变量） $g(U)$ 的均值问题。随机变量的均值（总体均值）问题可以通过模拟生成一些样本，利用样本均值来对总体均值进行估计，即用 $\sum _{i=1}^k\frac{g(U_i)}{k}$ 估计 $E[g(U)]$ 。

上节介绍的方法有个缺点： 如果不同随机数 $U$ 对应的 $g(U)$ 相差很大，则“样本均值”的方差会比较大（估计值不稳定），这时用“样本均值”做估计值是不合适的。当然我们可以通过增加样本数来使方差减小（估计值变稳定），重要性抽样法提供了另外的方法来使估计值变稳定（方差减小）。

重要性抽样法的思想是：模拟具有概率密度函数 $f(x)$ 的随机变量 $X$ ，设 $Y=\frac{g(X)}{f(X)}$ ，这时由随机变量函数期望的性质可得：
$$E(Y)=\int \frac{g(x)}{f(x)}f(x)\mathrm{d}x=\int g(x)\mathrm{d}x$$
所以，我们通过模拟生成 $Y$ 的样本，利用 $Y$ 的样本均值即可估计原始问题 $\int g(x)\mathrm{d}x$ 。其中， 函数 $f(x)$ 称为 重要性概率密度函数(importance density function)，或简称重要性函数。选择重要性概率密度函数 $f(x)$ 的一个要求是当 $g(x)>0$ 时 $f(x)>0$ （此时 $f(x)$ 不能等于 0，否则上式不成立）；另一个要求是 $Y=\frac{g(X)}{f(X)}$ 的方差要小（这是重要性抽样法具有价值的关键所在），从而在模拟时可以更快地得到一个稳定的估计值。

说明：上节介绍的简单模拟方法其实就是重要性抽样法中重要性概率密度函数为均匀分布时特例。

下面例子摘自：http://www.aw-bc.com/greenwell/markov.pdf

在社会学中，人按其收入可分为下层(lower-class)，中层(middle-class)和上层(upper-class)这三个阶层。社会学发现“当前这一代的所在阶层”严重影响着“下一代的阶层”。下表是父代和子代阶层转移的概率矩阵如图 7 所示。

也可以记为矩阵的形式，如：
$$P=\left[\begin{array}{ccc}0.65& 0.28& 0.07\\ 0.15& 0.67& 0.18\\ 0.12& 0.36& 0.52\end{array}\right]$$
一般，我们用 $p_{ij}$ 表示从状态 $i$ 转移到状态 $j$ 的概率。

这个转移概率矩阵很有用。比如，我们可以利用转移概率矩阵计算“已知父亲为上层，孙子为中层的概率”。计算过程如图 8 所示。

已知父亲为上层，孙子为中层的概率为： $p_{31}\cdot p_{12}+p_{32}\cdot p_{22}+p_{33}\cdot p_{32}=0.0336+0.2412+0.1872=0.4620$ 。

一般地，假设上一代“下中上”三个阶层的比例为概率分布向量 ${\mathit{\pi }}_{\mathit{n}\mathbf{-}\mathbf{1}}=[{\pi }_{n-1}(1),{\pi }_{n-1}(2),{\pi }_{n-1}(3)]$ ，则这一代“下中上”三个阶层的比例将为： ${\mathit{\pi }}_{\mathit{n}}={\mathit{\pi }}_{\mathit{n}\mathbf{-}\mathbf{1}}P$ 。

下面我们将会观察到一个有趣的现象——不管初始时“下中上”三个阶层的所占比例是多少，经过多次迭代后，处于“下中上”三个阶层的比例会稳定不变。
例如，“下中上”三个阶层初始比例为概率分布向量 ${\mathit{\pi }}_{\mathbf{0}}=[0.210,0.680,0,110]$ ，多次迭代后的比例为（由于有舍入误差，有时比例和会不为 1）图 9 所示。

又如，“下中上”三个阶层初始比例为概率分布向量 ${\mathit{\pi }}_{\mathbf{0}}=[0.75,0.15,0,1]$ ，多次迭代后的比例为图 10 所示。

最终的概率分布不仅收敛，而且和初始概率分布无关。初始概率分布不同时，都收敛到同一个概率分布 $\mathit{\pi }=[0.286,0.489,0.225]$ 。 最后收敛到的这个分布称为马尔可夫链的极限分布。

但并不是所有的马尔可夫链都存在极限分布。比如转移概率矩阵为：
$$P=\left[\begin{array}{cccc}0& 1/2& 0& 1/2\\ 1/2& 0& 1/2& 0\\ 0& 1/2& 0& 1/2\\ 1/2& 0& 1/2& 0\end{array}\right]$$
的马尔可夫链不会收敛，不存在极限分布。

满足什么条件时，马尔可夫链存在极限分布呢？
定理：
设马尔可夫链的状态空间为 $I=\{a_1,a_2,\cdots ,a_N\}$ ， $P$ 为马尔可夫链的转移概率矩阵（矩阵中元素用 $p_{ij}$ 表示）。记 $p_{ij}(n)$ 为从状态 $i$ 经过 $n$ 次转移到达状态 $j$ 的概率。如果存在正整数 $m$ ，使对任意 $a_i,a_j\in I$ ，都有：
$$p_{ij}(m)>0,i,j=1,2,\cdots ,N$$
则称此马尔可夫链具有 遍历性 ，且极限分布（又称平稳分布）存在。且极限分布 $\mathit{\pi }=[{\pi }_1,{\pi }_2,\cdots ,{\pi }_N]$ 可以从下面方程组中求得：
$$\{\begin{array}{ll}& {\pi }_j=\sum _{i=1}^N{\pi }_i{p}_{ij},j=1,2,\cdots ,N\\ & {\pi }_j>0,j=1,2,\cdots ,N\\ & \sum _{j=1}^N{\pi }_j=1\end{array}$$
这个定理的证明略。

由这个定理知，要证明马尔可夫链是遍历的，只需要找一正整数 $m$ ，使 $m$ 步转移概率矩阵 $P^m$ 中无零元素（全部大于零）。

前面抽样算法可以工作，但如果接受率 $\alpha (i,j)$ 比较小，则采样时马尔可夫链容易原地踏步，拒绝大量的跳转，这使得马尔可夫链遍历整个状态空间要花费太长的时间，收敛到平稳分布 $\pi (x)$ 的速度太慢。下面我们将提高想办法提高接受率。

如果 $\alpha (i,j)=0.1,\alpha (j,i)=0.2$ ，此时的细致平稳条件为：
$$\pi (i)q(j\mid i)\times 0.1=\pi (j)q(i\mid j)\times 0.2$$
显然， $\alpha (i,j)$ 和 $\alpha (j,i)$ 同比例放大，使得两数中最大的一个放大到 1 时，细致平稳条件仍然会满足。如把上式两边扩大 5 倍，有：
$$\pi (i)q(j\mid i)\times 0.5=\pi (j)q(i\mid j)\times 1$$
“两数中最大的一个放大到 1”用公式表达就是：
$$\alpha (i,j)=min\{1,\frac{\pi (j)q(i\mid j)}{\pi (i)q(j\mid i)}\}$$

这就是 Metropolis-Hastings算法 。
(1) 寻找另一个辅助的概率分布（称为建议分布） $q(y\mid x)$ ，初始化马尔可夫链 $X_0=x_0$ 。
(2) 对 $t=0,1,2,\cdots$ ，用下面方法迭代进行抽样。
(2.1) 记 $t$ 时刻状态为 $X_t=x_t$ ，从建议分布中抽样得到样本 $y$ ，即 $y\sim q(x\mid x_t)$ 。
(2.2) 产生一个(0,1)之间的随机数 $u$ 。
(2.3) 如果 $u<\alpha (x_t,y)=min\{1,\frac{\pi (y)q(x_t\mid y)}{\pi (x_t)q(y\mid x_t)}\}$ ，则接受转移 $x_t\to y$ ，即 $X_{t+1}=y$ ；否则不接受转移，即 $X_{t+1}=x_t$ 。

这样，得到的样本 $X$ 经过一定迭代次数后达到稳定状态，这时它就符合概率分布 $\pi (x)$ 。

说明 1：Metropolis-Hastings 算法得到的样本并不独立，但是我们所要求的是得到的样本符合给定的概率分布，并不要求独立。
说明 2：Metropolis-Hastings 算法中的“接受”和拒绝抽样或 SIR 中的“接受”不是同一回事。Metropolis-Hastings 算法中的“接受”是指接受这个状态转移（如果不接受则意味着采用当前状态作为下一个状态，不管“接受”还是“不接受”都得到了一个样本值），而拒绝抽样或 SIR 中的“接受”是指接受样本值（如果不接受则意味着相关程序还要运行一次才可能得到一个样本值）。
说明 3：对于高维空间的 $\pi (\mathit{x})$ ，如果满足细致平稳条件
$$\pi (\mathit{x})Q\prime (\mathit{x}\to \mathit{y})=\pi (\mathit{y})Q\prime (\mathit{y}\to \mathit{x})$$
以上的 Metropolis-Hastings 算法一样有效。

参考：http://cos.name/2013/01/lda-math-mcmc-and-gibbs-sampling/

这里采用前面介绍 SIR 算法时所使用的实例。

假设需要模拟的随机变量概率密度为： $f(x)=\frac{1}{c}h(x)$ ，其中 $c$ 为未知的常数， $h(x)=e^{0.4\ast (x-0.4{)}^2-0.08\ast x^4}$ 。

k <- function(x, a=.4, b=.08){exp(a*(x-a)^2 - b*x^4)}  # k(x)本身不是概率密度函数，它是概率密度函数f(x)乘以某个未知常数
x <- seq(-4, 4, 0.1)

N <- 100000
x.acc5 <- rep.int(NA, N)
u <- runif(N)
acc.count <- 0
std <- 0.05        ## Spread of proposal distribution
xc <- 0;           ## Starting value
for (ii in 1:N){
    xp <- rnorm(1, mean=xc, sd=std)   ## proposal distribution
    alpha <- min(1, (k(xp)/k(xc)) * (dnorm(xc, mean=xp,sd=std)/dnorm(xp, mean=xc,sd=std)))
    x.acc5[ii] <- xc <- ifelse(u[ii] < alpha, xp, xc)
}

pdf(file="myplot_MH.pdf")
par(mfrow=c(2,2), mar=c(2,2,1,1))
plot(x,k(x),type="l")
barplot(table(round(x.acc5,1))/length(x.acc5))

上面程序生成的 myplot_MH.pdf 如图 11 所示。

参考：http://people.math.aau.dk/~kkb/Undervisning/Bayes14/sorenh/docs/sampling-notes.pdf