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

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

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

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

θ=abf(x)dx(1) \theta=\int_a^b f(x) \mathrm{d}x\tag{1}

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

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

随机投点法

Figure 1: 随机投点法

随机投点法

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

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

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

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

我们看每次投点的结果其实都是服从伯努利分布(也叫0-1分布),b(1,p)b(1,p),而nn次投点结果n0n_0则服从二项分布b(n,p)b(n,p)。显然,θ^\hat\thetaθ\theta的无偏估计: E(θ^)=E[M(ba)n0n]=M(ba)nE[n0]=M(ba)n×np=M(ba)p=θ(4) 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} 在无偏的情况下,我们可以用方差衡量精度: Var(θ^)=Var[M(ba)n0n]=M2(ba)2n2Var(n0)=M2(ba)2n2×np(1p),其中p=θM(ba)=1nθ[M(ba)θ](5) \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} 其中,θ,M,b,a\theta, M, b, a都是定值,只有nn是我们能够影响的数。也就是说,模拟次数nn提高到nn倍,方差缩小到1/n1/n,如果看和精度直接相关的标准差形如(E[x]±std(x)E[x]\plusmn \mathrm{std}(x)),那么精度每增加一位小数,实验量需要增加100倍。

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

随机投点法是一种很直觉,但是效率不高的方法,另一种效率更高的方法是利用期望值(矩)的估计。取随机变量UU(a,b)U\sim U(a,b),则将UU代入f(x)f(x),得到一个随机数Y=f(U)Y=f(U)E[Y]=E[f(U)]=abf(u)1badu=θbaθ=(ba)E[Y]=(ba)E[f(U)](6) 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} 在范围bab-a已知的情况下,我们需要知道f(U)f(U)的期望即可。我们知道一个随机变量的期望即为一阶原点矩(一阶矩),因此我们可以用矩估计的方法来近似计算。一阶矩的估计方式很简单,就是采样多个点,求函数均值,用样本均值代替一阶矩(总体期望)

我们取服从均匀分布U(a,b)U(a,b)的独立的随机变量{Ui,i=1,2,,n}\{U_i,i=1,2,\dotsb,n\},则Y=f(U)Y=f(U)的n个采样值为{Yi=f(Ui),i=1,2,,n}\{Y_i=f(U_i),i=1,2,\dotsb,n\},根据矩估计方法, Y^=1ni=1nf(Ui)E[Y]=E[f(U)] a.s.n(7) \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)可得: θ^=(ba)Y^=bani=1nf(Ui)(ba)E[f(U)]=θ a.s.n(8) \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)f(x)上取足够多的点,用采样点的平均值代替函数值的均值,再乘以定义域长度,得到面积(积分值)。用一个图描述了蒙特卡洛求定积分的思想:

蒙特卡洛求定积分.png

Figure 2: 蒙特卡洛求定积分.png

蒙特卡洛求定积分.png

显然,样本平均值法也是无偏估计。我们也可以用方差计算其精度。为了方便比较,我们令随机投点法的积分估计值为θ^1\hat\theta_1,样本平均值法的估计值为θ^2\hat\theta_2。那么可以计算Var(θ^2)\mathrm{Var}(\hat\theta_2),根据式(8),有: Var(θ^2)=Var[bani=1nf(Ui)]=(ban)2Var[i=1nf(Ui)]Uii.i.d根据独立随机变量和的方差性质有:=(ban)2×n×Var[f(U)]=(ba)2n×ab(f(U)E[f(U)])21badUE[f(U)]=θba=1n[(ba)2abf2(U)1badUθ2]=1n[(ba)abf2(U)dUθ2](9) \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}

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

Var(θ^1)Var(θ^2)=1nθ[M(ba)]=1n[(ba)abf2(U)dUθ2]=M(ba)n[θabf2(U)MdU]0f(x)M0<f(U)M1M(ba)n[θabf(U)dU]=0 \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}\{U:f(U)<M\}不是零测集,就有Var(θ^1)>Var(θ^2)\mathrm{Var}(\hat\theta_1)>\mathrm{Var}(\hat\theta_2)。由此,样本均值法比随机投点法方差更小,更有效。此外,样本均值法不要求f(x)f(x)有上界MM,可以推广至一般情形。

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

我们再回看式(1),对他进行适当变形: θ=abf(x)dx=abf(x)g(x)g(x)dx=E[f(X)g(X)]g(x)0;x[a,b](10) \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)是一个概率密度函数。那么原来的积分问题就变成了求期望,若有来自g(x)g(x)的样本(随机数){x1,x2,,xn}\{x_1,x_2,\dotsb,x_n\}(服从g(x)g(x)的分布),则θ\theta的估计可以用一阶矩估计来实现: θ^=1ni=1nf(xi)g(xi)(11) \hat\theta=\frac{1}{n}\sum_{i=1}^n \frac{f(x_i)}{g(x_i)}\tag{11} 显然,如果当g(x)=1bag(x)=\frac{1}{b-a}时,这种方法就是样本均值法,也可以说样本均值法是此类方法的一个特例

估计量θ^\hat{\theta}构造是很简单的,说白了就是一阶原点矩估计。我们对θ^\hat\theta求期望易知此估计式无偏的,即E[θ^]=θE[\hat{\theta}]=\theta。而其方差为: Var(θ^)=1n[E(f(x)g(x))2E2(f(x)g(x))]=1n[E(f(x)g(x))2θ2](12) \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)g(x)使蒙特卡洛积分的效率更高。理论上来看,当g(x)=f(x)abf(x)dxg(x)=\frac{f(x)}{\int_a^b f(x) \mathrm{d}x}时,有Var(θ^)=0\mathrm{Var}(\hat\theta)=0,此时结果最优,但是这个g(x)g(x)中必须已知式(1)abf(x)dx\int_a^b f(x) \mathrm{d}x的结果,显然是无法直接使用的。不过,他也给我们一个启示,即g(x)g(x)f(x)f(x)应该在形状相近的情况下,估计的方差更小。下面我们举一个例子: f(x)=x2+2xx[0,2]g1(x)=12;g2(x)=ex;g3(x)=12π×0.5e(x1)22×0.52 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)g(x)函数,作为重要性采用的概率密度函数。首先,我们易算出f(x)=x2+2xx[0,2]f(x)=-x^2+2x\quad x\in[0,2]的积分结果为43\frac{4}{3},我们画出各个函数的图像以及f(x)g(x)\frac{f(x)}{g(x)}的图像。

importance sampling

Figure 3: importance sampling

importance sampling

Figure 4: importance sampling

我们使用matlab分别计算三个重要性采样的方差有(取n=10n=10): Var(110i=110f(x)g1(x))=0.03556Var(110i=110f(x)g2(x))=0.1335Var(110i=110f(x)g3(x))=0.01205 \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计算结果显示g3(x)g_3(x)能够缩减方差,而与f(x)f(x)形状差异大的g2(x)g_2(x)反而使方差更大了。

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

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

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

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

我们继续以上个例子f(x)=x2+2xf(x)=-x^2+2x来说明,我们把积分区间[0,2][0,2]划分成十等分,根据式(13)有以下表格: |区间|积分值pip_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

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

分层抽样1

Figure 5: 分层抽样1

分层抽样1

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

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

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

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