概率统计随机过程之蒙特卡洛方法

• 蒙特卡洛求定积分的基本方法
• 随机投点法
• 样本平均值法
• 重要性采样
• 分层抽样法
• 最大期望算法（Expectation Maximum, EM）算法
• 马尔科夫链蒙特卡洛MCMC

概率统计随机过程之蒙特卡洛方法

蒙特卡方法洛（Monte Carlo Method）简称MC，又叫随机模拟方法，本质是随机抽样。MC方法的使用框架一般为

1. 针对实际问题建立一个便于实现的概率统计模型，使所求解恰好是所建模型的概率分布或其某个数字特征，比如某个事件的概率或期望值。
2. 对模型中随机变量建立抽样，在计算机上进行随机试验，抽取足够的随机数，并对有关事件进行统计；
3. 对模拟实验结果加以分析，给出所求解的估计及精度（方差）的估计；
4. 必要时，还应改进模型以降低估计方差和减少实验复杂度，提供模拟计算效率。

蒙特卡洛求定积分的基本方法

蒙特卡洛方法近年来受到了广泛运用，不少统计问题到最后到可以归结为定积分的计算，如计算概率、各阶矩、贝叶斯机器学习等，因此我们首先介绍定积分的蒙特卡洛计算。考虑定积分的一个通用形式：

\[ \theta=\int_a^b f(x) \mathrm{d}x\tag{1} \]

我们先考虑两种基础的蒙特卡洛积分方法，然后再给出一些给复杂而有效的办法。

随机投点法

随机投点法底层原理是伯努利大数定理（用频率接近概率），可以简单理解成几何概型的应用。简单起见，假设\(a,b\)有限，式（1）中的函数满足\(0\leq f(x)\leq M\)，令\(\Omega=\{(x,y):a\leq x\leq b,0\leq y\leq M\}\)。那么根据几何概型，在\(\Omega\)内取均匀分布的点，该点落在\(f(x)\)下方的概率为\(p=\frac{S(\Omega)}{\theta}\)，其中\(S(\Omega)\)是整体的面积，值为\(M(b-a)\)，\(\theta\)是\(f(x)\)函数下方的面积，根据积分的定义可知\(\theta=\int_a^b f(x) \mathrm{d}x\)。

随机投点法

显然，这就是高中学过的几何概型，概率等于面积的比值，由此，我们也可以反推\(f(x)\)下方的面积，即积分的值为： \[ \theta=\int_a^b f(x) \mathrm{d}x=M(b-a)p\tag{2} \] 现在不知道的就是概率\(p\)。根据大数定律，如果我们往\(\Omega\)中投很多点\(n\)，若有\(n_0\)个点落在\(f(x)\)下方，那么可以通过频率逼近概率，即\(\hat p=\frac{n_0}{n}\rightarrow p (a.s.\; n\rightarrow\infty)\)，代入式（2）中有 \[ \hat{\theta}=M(b-a)\frac{n_0}{n}\rightarrow \theta(a.s.\; n\rightarrow\infty) \tag{3} \] 上述思想是容易实现的，步骤如下

1. 独立地产生\(n\)对服从\(U(0,1)\)的独立随机数\((u_i,v_i),i=1,2,\dotsb,n\)；
2. 根据随机变量的函数关系有\(x_i=a+u_i(b-a),y_i=Mv_i\)；
3. 统计\(y_i\leq f(x_i)\)的个数\(n_0\)(该点落在\(f(x)\)下方)；
4. 使用式（3）估算积分值。

随机投点法用了类似于舍选法的做法， 在非随机问题中引入随机性时用了二维均匀分布和二项分布，靠求二项分布概率来估计积分，随机投点法容易理解，但是效率较低。

那么，随机投点法精度提高一位数需要多大的代价呢？

我们看每次投点的结果其实都是服从伯努利分布（也叫0-1分布），\(b(1,p)\)，而\(n\)次投点结果\(n_0\)则服从二项分布\(b(n,p)\)。显然，\(\hat\theta\)是\(\theta\)的无偏估计: \[ E(\hat\theta)=E[M(b-a)\frac{n_0}{n}]=\frac{M(b-a)}{n}E[n_0]\\ =\frac{M(b-a)}{n}\times np=M(b-a)p=\theta\tag{4} \] 在无偏的情况下，我们可以用方差衡量精度： \[ \begin{aligned} Var(\hat\theta)&=Var[M(b-a)\frac{n_0}{n}]=\frac{M^2(b-a)^2}{n^2}Var(n_0)\\ &=\frac{M^2(b-a)^2}{n^2}\times np(1-p),其中p=\frac{\theta}{M(b-a)}\\ &=\frac{1}{n}\theta[M(b-a)-\theta] \end{aligned}\tag{5} \] 其中，\(\theta, M, b, a\)都是定值，只有\(n\)是我们能够影响的数。也就是说，模拟次数\(n\)提高到\(n\)倍，方差缩小到\(1/n\)，如果看和精度直接相关的标准差形如（\(E[x]\plusmn \mathrm{std}(x)\)），那么精度每增加一位小数，实验量需要增加100倍。

随机模拟积分的精度一般都服从这样的规律。

样本平均值法

随机投点法是一种很直觉，但是效率不高的方法，另一种效率更高的方法是利用期望值(矩)的估计。取随机变量\(U\sim U(a,b)\)，则将\(U\)代入\(f(x)\)，得到一个随机数\(Y=f(U)\)： \[ E[Y]=E[f(U)]=\int_a^b f(u) \frac{1}{b-a} \mathrm{d}u=\frac{\theta}{b-a}\\ \Rightarrow\theta=(b-a)E[Y]=(b-a)E[f(U)]\tag{6} \] 在范围\(b-a\)已知的情况下，我们需要知道\(f(U)\)的期望即可。我们知道一个随机变量的期望即为一阶原点矩（一阶矩），因此我们可以用矩估计的方法来近似计算。一阶矩的估计方式很简单，就是采样多个点，求函数均值，用样本均值代替一阶矩（总体期望）。

我们取服从均匀分布\(U(a,b)\)的独立的随机变量\(\{U_i,i=1,2,\dotsb,n\}\)，则\(Y=f(U)\)的n个采样值为\(\{Y_i=f(U_i),i=1,2,\dotsb,n\}\)，根据矩估计方法， \[ \hat Y={1\over n}\sum_{i=1}^nf(U_i) \rightarrow E[Y]=E[f(U)]\ a.s. n\rightarrow \infty\tag{7} \] 代入式（6）可得： \[ \hat\theta=(b-a)\hat Y=\frac{b-a}{n}\sum_{i=1}^nf(U_i) \rightarrow (b-a)E[f(U)]=\theta\ a.s. n\rightarrow\infty\tag{8} \]

简单来说，就是从\(f(x)\)上取足够多的点，用采样点的平均值代替函数值的均值，再乘以定义域长度，得到面积（积分值）。用一个图描述了蒙特卡洛求定积分的思想：

蒙特卡洛求定积分.png

显然，样本平均值法也是无偏估计。我们也可以用方差计算其精度。为了方便比较，我们令随机投点法的积分估计值为\(\hat\theta_1\)，样本平均值法的估计值为\(\hat\theta_2\)。那么可以计算\(\mathrm{Var}(\hat\theta_2)\)，根据式（8），有： \[ \begin{aligned} \mathrm{Var}(\hat\theta_2)&=\mathrm{Var}[\frac{b-a}{n}\sum_{i=1}^nf(U_i)]=(\frac{b-a}{n})^2\mathrm{Var}[\sum_{i=1}^nf(U_i)]\\ &\because U_i\quad i.i.d\quad \therefore 根据独立随机变量和的方差性质有：\\ &=(\frac{b-a}{n})^2\times n \times \mathrm{Var}[f(U)]\\ &=\frac{(b-a)^2}{n}\times \int_a^b (f(U)-E[f(U)])^2 \frac{1}{b-a} \mathrm{d}U\\ &\because E[f(U)]=\frac{\theta}{b-a}\\ &=\frac{1}{n}[(b-a)^2\int_a^b f^2(U)\frac{1}{b-a}\mathrm{d}U-\theta^2]\\ &=\frac{1}{n}[(b-a)\int_a^b f^2(U)\mathrm{d}U-\theta^2] \end{aligned}\tag{9} \]

直观来看无法比较随机投点法和样本均值法的差异，我们将\(\hat\theta_1,\hat\theta_2\)做差，进行比较，在\(0\leq f(x) \leq M\)时（随机投点法要求），可以证明\(\mathrm{Var}(\hat\theta_1)\geq \mathrm{Var}(\hat\theta_2)\)，二者相减结果如下：

\[ \begin{aligned} \mathrm{Var}(\hat\theta_1)-\mathrm{Var}(\hat\theta_2)&=\frac{1}{n}\theta[M(b-a)]-\\ &=\frac{1}{n}[(b-a)\int_a^bf^2(U)\mathrm{d}U-\theta^2]\\ &=\frac{M(b-a)}{n}[\theta-\int_a^b \frac{f^2(U)}{M}\mathrm{d}U]\\ &\because 0\leq f(x) \leq M \therefore 0<\frac{f(U)}{M}\leq 1\\ &\geq \frac{M(b-a)}{n}[\theta-\int_a^b f(U)\mathrm{d}U]=0\\ \end{aligned} \] 由此可见，只要\(\{U:f(U)<M\}\)不是零测集，就有\(\mathrm{Var}(\hat\theta_1)>\mathrm{Var}(\hat\theta_2)\)。由此，样本均值法比随机投点法方差更小，更有效。此外，样本均值法不要求\(f(x)\)有上界\(M\)，可以推广至一般情形。

重要性采样

我们根据随机投点法和样本均值法已经能够进行蒙特卡洛积分的无偏估计，并且随着采样次数的增加最终能以\(O(n^{-1/2})\)收敛。我们现在需要的是优化蒙特卡洛积分的计算方法，使估计的方差尽量小。

我们再回看式（1），对他进行适当变形： \[ \theta=\int_a^b f(x) \mathrm{d}x=\int_a^b \frac{f(x)}{g(x)} g(x) \mathrm{d}x=E[\frac{f(X)}{g(X)}]\\ g(x)\neq 0;x\in[a,b]\tag{10} \] 其中，\(g(x)\)是一个概率密度函数。那么原来的积分问题就变成了求期望，若有来自\(g(x)\)的样本（随机数）\(\{x_1,x_2,\dotsb,x_n\}\)（服从\(g(x)\)的分布），则\(\theta\)的估计可以用一阶矩估计来实现： \[ \hat\theta=\frac{1}{n}\sum_{i=1}^n \frac{f(x_i)}{g(x_i)}\tag{11} \] 显然，如果当\(g(x)=\frac{1}{b-a}\)时，这种方法就是样本均值法，也可以说样本均值法是此类方法的一个特例。

估计量\(\hat{\theta}\)的构造是很简单的，说白了就是一阶原点矩估计。我们对\(\hat\theta\)求期望易知此估计式无偏的，即\(E[\hat{\theta}]=\theta\)。而其方差为： \[ \mathrm{Var}(\hat\theta)=\frac{1}{n}[E(\frac{f(x)}{g(x)})^2-E^2(\frac{f(x)}{g(x)})]=\frac{1}{n}[E(\frac{f(x)}{g(x)})^2-\theta^2]\tag{12} \] 那么根据式（12），除了增加采样的次数，我们又多了一个可以用来降低方差的工具\(g(x)\)，我们可以通过合理设计\(g(x)\)使蒙特卡洛积分的效率更高。理论上来看，当\(g(x)=\frac{f(x)}{\int_a^b f(x) \mathrm{d}x}\)时，有\(\mathrm{Var}(\hat\theta)=0\)，此时结果最优，但是这个\(g(x)\)中必须已知式（1）\(\int_a^b f(x) \mathrm{d}x\)的结果，显然是无法直接使用的。不过，他也给我们一个启示，即\(g(x)\)与\(f(x)\)应该在形状相近的情况下，估计的方差更小。下面我们举一个例子： \[ f(x)=-x^2+2x\quad x\in[0,2]\\ g_1(x)=\frac{1}{2};g_2(x)=e^{-x};g_3(x)=\frac{1}{\sqrt{2\pi}\times 0.5}e^{-\frac{(x-1)^2}{2\times 0.5^2}} \] 分别用3个\(g(x)\)函数，作为重要性采用的概率密度函数。首先，我们易算出\(f(x)=-x^2+2x\quad x\in[0,2]\)的积分结果为\(\frac{4}{3}\)，我们画出各个函数的图像以及\(\frac{f(x)}{g(x)}\)的图像。

我们使用matlab分别计算三个重要性采样的方差有（取\(n=10\)）： \[ \mathrm{Var}(\frac{1}{10}\sum_{i=1}^{10} \frac{f(x)}{g_1(x)})=0.03556\\ \mathrm{Var}(\frac{1}{10}\sum_{i=1}^{10} \frac{f(x)}{g_2(x)})=0.1335\\ \mathrm{Var}(\frac{1}{10}\sum_{i=1}^{10} \frac{f(x)}{g_3(x)})=0.01205 \] matlab计算结果显示\(g_3(x)\)能够缩减方差，而与\(f(x)\)形状差异大的\(g_2(x)\)反而使方差更大了。

可以从第二个图中看出，如果\(g(x)\)与\(f(x)\)的形状相似，相除之后的函数图像更加贴近积分结果\(\theta=\frac{4}{3}(\frac{f(x)}{g_3(x)})\)，这样我们在抽样时结果接近\(\theta\)的概率也更大；如果形状相差太大，如\(g_2(x)\)积分的结果反而会变得更加不好统计。

综合来说，重要性采样让新的函数值\(h(x)=\frac{f(x)}{g(x)}\)取到接近于积分结果的概率变大了。我们现在的问题就转变成了如何找到一个和\(f(x)\)形状相似的函数。

分层抽样法

在重要性采样的结论中，我们知道当\(f(x)\)与\(g(x)\)比值为常数时，可以得到理解最小的方差0，那么我们可不可以近似的构造一个这样的函数呢？基于这个思想，我们有了分层抽样法。

分层抽样法首先把样本空间\(D\)分成一些小区间\(D_1,\dotsb,D_m\)，且诸\(D_i\)不交，\(\cup D_i =D\)，然后在各小区间内的抽样数由其“贡献”大小决定。这里的“贡献”我们定义为在小区间\(D_i\)内的积分值，即 \[p_i=\int_{D_i}f(x)\mathrm{d}x\tag{13}\] 在区间\(D_i\)抽样数与\(p_i\)成正比。显然，分层抽样法是利用小区间构造一个离散的形状类似于原函数\(f(x)\)的分段函数。

我们继续以上个例子\(f(x)=-x^2+2x\)来说明，我们把积分区间\([0,2]\)划分成十等分，根据式（13）有以下表格： |区间|积分值\(p_i\)|归一化比例| |:-:|:-:|:-:| |[0,0.2)|0.0373|0.0280| |[0.2,0.4)|0.1013|0.0760| |[0.4,0.6)|0.1493|0.1120| |[0.6,0.8)|0.1813|0.1360| |[0.8,1.0)|0.1973|0.1480| |[1.0,1.2)|0.1973|0.1480| |[1.2,1.4)|0.1813|0.1360| |[1.4,1.6)|0.1493|0.1120| |[1.6,1.8)|0.1013|0.0760| |[1.8,2.0]|0.0373|0.0280

我们将归一化\(p_i\)值与\(f(x)\)画到一张图上(图中黄线与红线)。此外，为了表现构造的分段函数与原函数形状是相似的，我们将分段函数线性放大6倍（蓝线）。

分层抽样1

可以看出，分层抽样构造出的函数与目标积分函数是类似的，我们区间划分的越细，分段函数越相似。根据分段函数比例分配抽样次数，可以大大降低方差。理论上，给定分层抽样函数后，每个区间的分配的抽样数最优方式为: \[ n_i=nl_i\sigma_i/(\sum_{i=1}^m l_i\sigma_i)\tag{14} \] 其中，\(n_i\)是区间\(D_i\)的抽样数，\(n\)为总抽样数，\(m\)为总区间数目，\(l_i\)是区间\(D_i\)的长度，\(\sigma_i\)是区间\(D_i\)的标准差，此时方差最小，为\(\frac{1}{n}(\sum_{i=1}^m l_i\sigma_i)^2\)。然而实际计算中，每个区间的标准差\(\sigma_i\)总是未知的，式(14)无法直接使用。即使如此，最简单的分配方案\(n_i=nl_i/\sum l_i\)方差也比样本均值法低。

此外，其他缩减方差技术还有关联抽样法、控制变量法、对立变量法、条件期望法等等。

最大期望算法（Expectation Maximum, EM）算法

我们在参数估计中常用最大似然估计，不过实际问题中，除了需要估计的未知参数，还有别的

马尔科夫链蒙特卡洛MCMC

马尔可夫链蒙特卡洛（英语：Markov chain Monte Carlo，MCMC）方法（含随机游走蒙特卡洛方法）是一组用马氏链从随机分布取样的算法，之前步骤的作为底本。步数越多，结果越好。