Monte Carlo Approximations

准备总结几篇关于 Markov Chain Monte Carlo 的笔记。

本系列笔记主要译自A Gentle Introduction to Markov Chain Monte Carlo (MCMC) 文章下给出的链接。

Monte Carlo Approximations

Monte Carlo Approximation for Integration

理论部分

本文主要参考 Monte Carlo Approximations

蒙特卡洛方法是用来近似计算积分的，通过数值方法也可以计算积分：最简单的近似方法是通过求小矩边梯形的面积再累加。但是当函数的维度太高，变量数太多的时候，数值方法就不适用了，蒙特卡洛方法从概率学的角度进行积分的近似计算。

我们经常需要计算下面形式的积分：

\(I = \int_{a}^b h(x)g(x)dx\)

对于某些关于随机变量 \(X\) 的函数 \(f(x)\), 其期望值可以通过下面的公式计算：

\(E[x] = \int_{p(x)}p(x)f(x)dx\)

其中 \(p(x)\) 为变量\(X\) 的密度函数，它在定义域上的积分应该为1。

事实上，如果 \(g(x)\) 满足下面的条件：

1. 在区间 \((a,b)\) 上非负

\(g(x)\geqslant0, x\in(a,b)\)

2. 在区间上的积分有限

\(\int_a^b g(x)=C<\infty\)

那么，\(g(x)\) 就可以转换为密度函数

\(p(x)=\frac{g(x)}{C}\)

只需要在积分号前乘以相应的\(C\)，就可以保证原积分不变了。那么修改成密度函数后有什么好处呢？修改成密度函数的好处就是：可以通过随机模拟符合这个密度分布的变量\(x\)，来求出原积分的近似值：

根据大数定理：

\(E[f(X)] = \int f(x)p(x)dx \approx \frac{1}{S}\sum_{n=1}^S f(x_s)\)

其中\(x_s\sim p(X)\), 这一点很关键，也就是说，生成的随机数必须要符合密度函数为 \(p(X)\) 的分布，只有这样，计算出来的平均值才是积分的近似值。

为了对这种方法进行理论分析，下面给出两个定理：

大数定理 对于一个已知分布的随机序列，当取样数趋于无穷时，其均值趋向于期望。

事实上，在日常生活中我们常常是这么做的，比如统计一年级某门成绩的期望，我们可以随机选一些学生进行抽测，用他们的平均成绩来近似估计该年级的成绩期望，选的学生越多，其均值就越逼近真实期望值。

在统计背景下，\(A(n)\) 是 \(B\) 的一个一致估计量（consistent estimator），如果满足在 \(n\) 趋于无穷时，\(A(n)\) 收敛于 \(B\)，也就是说：对任意概率 \(P[0<P<1]\)，任意值 \(\delta\)，我们都可以找到 \(k\)，使得对任意 \(n>k\)， \(A(n)\) 与 \(B\) 的差距在 \(\delta\) 内的概率大于 \(P\)。

中心极限定理 大量独立随机变量的和近似服从正态分布。

假如我们有独立同分布(i.i.d)的随机变量 \(z_i\)，均值为 \(\mu\)，方差为 \(\sigma^2\)，n 个这样的变量的和记为 Y:

\(E(Y) = E(\sum\limits_{i=1}^n z_i) = \sum\limits_{i=1}^nE(z_i) = n\mu\)

\(var(Y) = var(\sum\limits_{i=1}^n z_i) = \sum\limits_{i=1}^n var(z_i) = n\sigma^2\)

对变量 Y 进行标准化：

\(\frac{Y-E(Y)}{\sqrt{var(Y)}}=\frac{Y-n\mu}{\sqrt{n}\sigma}\sim\mathcal{N}(0,1)\ as\ n\rightarrow \infty\)

也就是：

\(\frac{\sum\limits_{i=1}^n z_i-n\mu}{\sigma\sqrt{n}}\sim\mathcal{N}(0,1)\ as\ n\rightarrow \infty\)

因此，蒙特克洛方法的误差为：

可以看出：

1. 如果 \(f(x)\) 的方差是有限的，那么 MC 估计是一致的

2. MC 估计的误差是渐进无偏的

3. MC 估计的误差近似正态分布的

4. MC 估计误差的标准差是 \(\sigma=\frac{\sqrt{Var[f(x)]}}{\sqrt{n}}\)

同时，可以发现，MC 估计对数据的维度天然免疫，无论维度多大，只需要按分布函数生成随机变量，计算随机变量的函数值，计算这些函数值的均值，即可得到对积分的估计。由于误差的标准差越小，估计越精确，因此，我们可以通过增大取样数目 \(n\) 的方式来增大 MC 估计的准确性。另一个增加准确性的思路是降低 \(f(x)\) 的方差，这类方法称为 variance-reducing techniques 这里(我)不(没)做(看)介(懂)绍。

实例

1. 通过蒙特卡洛方法计算圆周率

通常的教程中介绍这个例子时会说：正方形中均匀生成的随机点落到圆内的概率为balabala，所以圆的面积比上正方形的面积为balabala，所以圆周率为balabala。这种方法固然好让人理解，但是如何严格地用公式表示这个过程呢？

圆的面积可以用下面的积分来获得：

其中在符合内部条件\(x^2+y^2\leq r^2\)时为1，在不符合内部条件时为0。也就是说，点在圆的内部时，函数值为1，点在圆外部时函数值为0。令 \(p(x)\)，\(p(y)\) 符合 \([-r,r]\) 上的均匀分布，因此，\(p(x)=p(y)=\frac{1}{2r}\)。

所以，根据蒙特卡洛积分，\(I = \int_{a}^b f(x)g(x)dx\) ,这里的 \(f(x)\) 转换为二维函数相当于 \(1(x^2+y^2\leq r^2)\)，\(g(x)\) 相当于 1，我们把它转化为二维的密度函数，也就是\(p(x,y)=p(x)\cdot p(y)=\frac{1}{4r^2}=\frac{g(x,y)}{4r^2}\)。

所以，

也就是，只需要生成 -r 到 r 均匀分布的横坐标和纵坐标，然后带入\(1(x^2+y^2\leq r^2)\) 判断其是否在圆内，在的话记为1，不在记为0，加和后除以点的总数，计算出落到圆内的概率，最后乘以 \(4r^2\) 求出圆的面积，知道了圆的面积，又知道半径了，自然也可以求出圆周率的近似值。

% Radious of circle
r = 1;
S = 10000;
% Generate S points uniformly distributed points
x = unifrnd(-r,r,S,1);
y = unifrnd(-r,r,S,1);
% The points falled into the circle
fxy = (x.^2+y.^2-r^2<=0);
area = 4*r^2/S*sum(fxy);
pi_estimate = area/r^2;

% Visualization
figure;
inside = fxy; outside = logical(1-fxy);
scatter(x(inside),y(inside),15,'r','filled');
hold on;
scatter(x(outside),y(outside),15, 'b','filled');
axis equal;
xlim([-1,1]);
ylim([-1,1]);

2. 通过蒙特卡洛方法近似计算\(xe^x\) 的积分　　

我们想计算 \(I = \int_0^1 xe^xdx\) , 由分部积分公式，不难得到其积分值为1。下面通过蒙特卡洛方法求其近似值，看是否与1接近。

可以把\(f(x)\) 看成 \(xe^x\)，\(g(x)\) 看成 1，由于积分区间长度也是1，故概率分布的密度函数为\(p(x)=\frac{gx}{\int_0^1 g(x)dx}=1\)。

所以，根据上文的讨论，只需要在0到1中均匀取得随机数，然后带入\(f(x)=xe^x\) 中计算函数值，最后求均值即可得到积分的近似。

代码如下：

rng('default');
% 100 sample
S1 = 100;
x1 = unifrnd(0,1,S1,1);
y1 = x1.*exp(x1);
I1 = sum(y1)/S1;

% 10000 samples
S2 = 10000;
x2 = unifrnd(0,1,S2,1);
y2 = x2.*exp(x2);
I2 = sum(y2)/S2;

3. 计算beta分布的期望

关于Beta分布，详细的资料请移步百度文库，这里只介绍一点：beta分布含有两个参数 A，B。对于随机变量 \(X\sim Beta(A,B)\)，\(E(x)=\frac{A}{A+B}\)。

计算 Beta 分布的期望，使用如下公式:

\(E(x)=\int_{p(x)}p(x)xdx\), where \(x\sim Beta(\alpha_1\alpha_2)\)

1. 首先我们发现 \(f(x)=x\)，而\(p(x)\) 本身就是概率密度函数，所以不用进行转换。

2. 我们根据概率分布 \(p(x)=Beta(\alpha_1,\alpha_2)\) 来生成大量随机的点，然后带入 \(f(x)=x\)，求出其均值，即为 Beta 分布的期望。

rand('seed',12345);
alpha1 = 2; alpha2 = 10;
N = 10000;
x = betarnd(alpha1,alpha2,1,N);%生成10000个按Beta分布的点

% MONTE CARLO EXPECTATION
expectMC = mean(x);%由MC方法算出的期望近似值

% ANALYTIC EXPRESSION FOR BETA MEAN
expectAnalytic = alpha1/(alpha1 + alpha2);%Beta分布的理论期望值

% DISPLAY
figure;
bins = linspace(0,1,50);%将[0,1]区间分成50份，最后一份是大于50的部分
counts = histc(x,bins);%统计每个小区间内落入的随机点的数目
probSampled = counts/sum(counts);%Beta分布的点落入每个小区间对应的概率
probTheory = betapdf(bins,alpha1,alpha2);%每个区间点的理论概率密度值
b = bar(bins,probSampled); colormap hot; hold on;%概率分布直方图，注意这里每个竖条代表的是对应区间的概率，不是概率密度。当它除以区间长度时才是概率密度
t = plot(bins,probTheory/sum(probTheory),'r','Linewidth',2)%这里probTheory/sum(probTheory)相当于是probTheory*dx，因为sum(probTheory*dx)=1, dx=1/sum(probTheory)
m = plot([expectMC,expectMC],[0 .1],'g')%期望的近似值
e = plot([expectAnalytic,expectAnalytic],[0 .1],'b')%期望的理论值
xlim([expectAnalytic - .05,expectAnalytic + 0.05])
legend([b,t,m,e],{'Samples','Theory','$\hat{E}$','$E_{Theory}$'},'interpreter','latex');
title(['E_{MC} = ',num2str(expectMC),'; E_{Theory} = ',num2str(expectAnalytic)])
hold off

从图中我们可以看到两条竖线十分接近，也就是说由样本估计的均值十分接近真实值。

Monte Carlo approximation for optimization

蒙特卡洛方法还可以用来求解优化问题：

\(x_{opt}=\arg\max g(x)\)

如果 \(g(x)\) 满足之前提到的两条性质：1. 非负 2. 积分有限

那么，就可以将其放缩为密度函数：

\(p(x) = \frac{g(x)}{C}\)

其中 \(C = \int_x g(x) dx\)

原优化问题变形为：

\(x_{opt}=\arg\max\limits_x Cp(x)\)

以 \(p(x)\) 为概率密度函数生成随机数，那么，在随机数密度最大的地方，对应的密度函数必然也是最大的，而密度函数与目标函数之间只差一个常数的放缩，故而该位置也是目标函数最大的地方，也就是说，这一点就是优化问题的近似解。

实例

求解优化问题 \( x_{opt}=\arg\max\limits_x e^{-\frac{(x-4)^2}{2}}\)

通过求导，令导函数为零，可以求得极值点为 x=4。

同时可以观察到，被优化函数\(g(x)\)与正态分布只差了一个参数：\( \sqrt{2\pi}\)

原问题可以变形为：

其中\(C=\sqrt{2\pi}\)。

因此，可以借助matlab生成均值为4，方差为1的正态分布随机变量，然后找出其密度函数最大的点对应的横坐标，此即为优化问题的解。

代码如下：

% MONTE CARLO OPTIMIZATION OF exp(x-4)^2
randn('seed',12345)

% INITIALZIE
N = 100000;
x = 0:.1:6;
C = sqrt(2*pi);
g = inline('exp(-.5*(x-4).^2)','x');
ph = plot(x,g(x)/C,'r','Linewidth',3); hold on%画出密度函数
gh = plot(x,g(x),'b','Linewidth',2); hold on;%画出待优化函数
y = normpdf(x,4,1);

% CALCULATE MONTE CARLO APPROXIMATION
x = normrnd(4,1,1,N);
bins = linspace(min(x),max(x),100);%生成100个随机区间
counts = histc(x,bins);%统计100个区间中每一个区间落入的随机点数
[~,optIdx] = max(counts);%找到落入点最多的区间的下标
xHat = bins(optIdx);%找到优化问题的近似解

% OPTIMA AND ESTIMATED OPTIMA
oh = plot([4 4],[0,1],'k');
hh = plot([xHat,xHat],[0,1],'g');

set(gca,'fontsize',16)
legend([gh,ph,oh,hh],{'g(x)','$p(x)=\frac{g(x)}{C}$','$x_{opt}$','$\hat{x}$'},'interpreter','latex','Location','Northwest');

　　

可以看到，通过蒙特卡洛算法得到的优化问题的解与真实解基本一致。