马尔科夫链蒙特卡洛采样(MCMC)入门 之二
目录
将概率模型应用到数据中,常需要复杂的推理过程,需要用到复杂的、高维的分布。马尔科夫链蒙特卡洛理论(MCMC)是一种通用的计算方法,通过迭代地对生成的样本进行求和代替复杂的数学推理。比较棘手的问题分析方法通常可以用某种形式的MCMC来解决,即便是高维的问题也同样如此。MCMC的发展可以说是统计学计算方法的最大进步。MCMC是一个非常活跃的研究领域,现在也有一些标准化的技术被广泛应用。我们讨论两种形式的MCMC方法:Metropolis-Hastings和Gibbs Sampling。在我们学习这些技术之前,我们首先要了解下面两种主要思想:蒙特卡洛积分和马尔科夫链(Monte Carlo integration, and Markov chains)。
2.1 蒙特卡洛积分(Monte Carlo integration)
概率推理中的很多问题需要复杂的积分计算或者对非常大的空间中的数据求和。例如,一个常见的问题是计算对于随机变量x的函数g(x)的期望值(简单起见,我们设是x单随机变量)。如果x是连续的,函数期望定义如下:
E[g(x)]=∫g(x)p(x)dx
如果x是离散变量,则积分被求和取代:
E[g(x)]=∑g(x)p(x)dx
很多情况下,我们想要计算统计量的均值和方差。例如,对g(x)=x,我们要计算分布的均值,使用积分或求和的分析技术对于某些特定分布具有很大的挑战性。例如,密度p(x)的函数可能不能进行积分。对于离散分布,可能由于结果空间太大而不能进行显式的求和。
蒙特卡洛积分方法一般的想法是使用样本近似估计复杂分布的期望。具体地,我们获得一系列样本x(t),t=1,...,N,这些样本从p(x)中独立获得。在这种情况下,我们可以使用有限样本的累加近似估计期望:
![](https://ask.qcloudimg.com/http-save/yehe-1565119/8au10c0qbg.png?imageView2/2/w/1620)
在上述过程中,我们用适当样本的求和(先求和再平均)来代替积分。一般来说,近似计算的精确度可以通过增加n来提高。另外需要注意的是,近似计算的精确度取决于样本之间的独立性。当样本相互关联的时候,有效样本的数量减少了,这是MCMC的存在的一个潜在的问题。
2.2 马尔科夫链
Markov链是一个随机过程,我们利用顺序过程从一个状态过渡到另一个状态。我们在状态x(1)开始Markov链,使用转移函数p(x(t)∣x(t−1))作为状态的转移矩阵,确定下一个状态x(2)。然后以x(2)作为开始状态使用同样的转移函数继续确定下一个状态,如此重复,得到一系列状态:
x(1)→x(2)→...→x(t)→...
这样一个状态序列就被称为马尔科夫链(Markov chain)。生成一个包含T个状态的马尔科夫链的过程如下:
步骤
设 t=1
生成一个初始值 u,并设 x(t)=u
重复下列过程
t=t+1 从转移函数 p(x(t)∣x(t−1)) 中采样一个新的值
设 x(t)=u
直到t=T
在上述迭代过程中,马尔科夫链的状态 t+1仅仅是根据上一个状态 t 产生的。因此,马尔科夫链向新的状态转移总是只依赖最后一个状态。这是MCMC中使用马尔科夫链的很重要的一个属性。
当初始化每个马尔科夫链的时候,马尔科夫链将会在状态空间中进行状态转移。如果我们开始一系列的马尔科夫链,每一个链赋予不同的开始状态,通过一系列的状态转移过程,最终每个链都将处于接近起始状态的某个状态,这个过程被称为“burn-in”。马尔科夫链的一个重要的性质是,链的起始状态经过足够次数的转换后最终的状态不会受初始状态的影响(假定满足马尔科夫链的某些条件),也就是说,马尔科夫链经过有限次的状态转移之后,最终能达到稳定的状态,被称为平稳分布,从其平稳分布中可以反映样本。当满足一定的条件时,无论从哪个状态开始,马尔科夫链最终都会到达一个确定的稳定状态。应用到MCMC中,它允许我们从一个分布中连续地采样,且序列的初始状态不会影响估计过程。
举例
举个例子:图2.1展示了一个马尔科夫链的例子,为简单起见,以单个连续变量为例。从Beta分布中进行采样,分布函数为
Beta(200(0.9x(t−1)+0.05),200(1−0.9x(t−1)−0.05))。
从四个不同的起始状态开始,每个链都有T=1000次迭代。图中的两个子图分别展示了不同时刻的状态序列,线的颜色代表不同的链。注意,前10次迭代显示了序列对初始状态的依赖,这是“burn-in”过程。接下来序列的剩余部分是稳定状态(如果不停止,那么继续保持稳定状态)。那么我们如何确切知道稳定状态到达的时间,又怎么知道链是收敛的呢?要解释上述问题并不是容易的,尤其是在高维状态空间中。
![](https://img2018.cnblogs.com/i-beta/13318/202002/13318-20200227094109510-1442044089.png)
2.3 综合讲解:马尔科夫链蒙特卡洛理论
前面两个部分讨论了MCMC理论背后的主要思想,蒙特卡洛采样和马尔科夫链。蒙特卡洛采样提供估计分布的各种特征:如均值,方差,峰值,或者其他的研究人员感兴趣的统计特征。马尔科夫链包含一个随机的顺序过程,并从平稳分布中采样状态。
MCMC的目标是设计一个马尔科夫链,该马尔科夫链的平稳分布就是我们要采样的分布,这就是所谓的目标分布。换句话说,我们希望从马尔科夫链的状态中采样等同于从目标分布中取样。这个想法是用一些方法设置转移函数,使无论马尔科夫链的初始状态是什么最终都能够收敛到目标分布。有很多方法使用相对简单的过程实现这个目标,我们在这里只讨论Metropolis,Metropolis Hastings和 Gibbs sampling。
2.4 Metropolis Sampling
Metropolis Sampling 是最简单的MCMC方法,是 Metropolis Hastings Sampling 的一个特例,Metropolis Hastings Sampling 将在下一节讨论。假设我们的目标是从目标密度函数 p(θ) 中采样,其中 −∞<θ<∞。Metropolis sampler 创建一个马尔科夫链并且产生一系列值:
θ(1)→θ(2)→...→θ(t)→...
其中 θ(t) 表示马尔科夫链在第 t 次迭代的状态。当“burn-in”过程之后,即达到平稳分布之后,从马尔科夫链中采样的样本反映了从目标分布 p(θ) 中的采样样本。
在Metropolis过程中,给第一个状态 θ(1) 初始化一些值。我们然后使用建议分布
![](https://img2018.cnblogs.com/i-beta/13318/202002/13318-20200227095108441-358482185.png)
(建议分布一般是非常简单的分布)生成一个候选节点θ*,该节点是在上一个状态的基础上生成的。下一步就是接受或者拒绝该建议。接受建议的概率是公式2.6:
![](https://img2018.cnblogs.com/i-beta/13318/202002/13318-20200227095055519-1295577584.png)
为了决定是否接受该建议,我们定义一个偏差变量u,如果u≤α,我们接受这个建议,并将下一个状态设置为该建议:
![](https://img2018.cnblogs.com/i-beta/13318/202002/13318-20200227095130382-1709262166.png)
步骤
设 t=1
生成一个初始值,并设 θ(t)=u
重复下列过程:
t=t+1
从 q(θ∣θ(t−1)) 中生成一个建议 θ∗
计算接受该建议的概率
![](https://img2018.cnblogs.com/i-beta/13318/202002/13318-20200227095243119-1998763573.png)
从均匀分布 Uniform(0,1) 中生成一个值u
如果 u≤α,接受该建议且设置 θ(t)=θ∗;否则设置 θ(t)=θ(t−1)
直到 t=T
图2.2说明了两个状态序列转移的过程。为了直观地了解该过程为什么能从目标分布中获得样本,注意,如果新的建议分布得到的新状态的建议值相比旧节点状态更有可能在目标分布之下,公式(2.6)将始终接受新的建议。因此,采样器将移动到状态空间中,目标函数密度更高区域。然而,如果新的建议没有当前的状态好,仍然有可能接受这个“更糟糕”的建议并且转移到该建议状态。该过程总是接受“好的”建议,并且偶尔接受“坏的”建议,以确保采样器能探索整个状态空间,并且确保采样的样本从分布的所有部分获得(包含尾部)。
Metropolis sampler一个关键的要求是:建议分布是必须对称的,即:
![](https://img2018.cnblogs.com/i-beta/13318/202002/13318-20200227095605459-345041787.png)
。因此,基于旧状态到新状态的建议概率,与从新状态返回到旧状态的建议概率是相同的,这样才能满足平稳细致方程(具体细节见:LDA数学八卦, 这里面有非常清晰说明平稳细致方程的含义,下图有个定义)。对称的建议分布有如下:正态分布(Normal),柯西分布(Cauchy),学生t分布(Student-t),以及均匀分布(uniform distributions)。如果建议分布不具有对称性,应该使用Metropolis-Hastings sampler,将在下一节讨论。
![](https://img2018.cnblogs.com/i-beta/13318/202002/13318-20200227095705023-1569085635.png)
图 摘自LDA数学八卦,第31页
Metropolis sampler的主要优点是,公式(2.6)只包含了一个密度函数的比例。因此在函数p(θ)中任何独立于θ的项都可以被忽略。因此,我们不需要知道密度函数的归一化常量。并且该采样规则允许从非标准分布中采样,这是非常重要的,因为非标准分布在贝叶斯模型中经常用到。
![](https://img2018.cnblogs.com/i-beta/13318/202002/13318-20200227095846995-985181057.png)
举例
举个例子:假设我们希望从柯西分布中随机采样,给出柯西分布的概率密度如下面公式(2.7):
![](https://img2018.cnblogs.com/i-beta/13318/202002/13318-20200227095936943-619811470.png)
如何使用Metropolis sampler来模拟这个分布了,采样得到符合这个分布的样本?
答:
因为在Metropolis sampler过程中我们不需要标准化常量,所以将上式重写为公式2.8,:
![](https://img2018.cnblogs.com/i-beta/13318/202002/13318-20200227095955843-1949718266.png)
因此,Metropolis 接受概率变为:
![](https://img2018.cnblogs.com/i-beta/13318/202002/13318-20200227100010056-1576540093.png)
我们使用正态分布作为建议分布,即从Normal(θ(t),σ)分布中产生建议。因此该分布的均值集中在当前状态,参数σ表示标准差需要人为设定,以控制建议分布的可变性。Listing 2.1展示了MATLAB函数,该函数返回非标准的柯西分布的密度。Listing 2.2展示了MATLAB代码实现的Metropolis sampler工具。图2.3展示了单个链经过500次迭代的模拟结果。其中,图2.3的第一个图展示了虚红线的理论密度,直方图展示了500个样本的分布。第二个图展示了样本的一个链。代码很简单,大家可以看看,不懂的可以在评论下方一起讨论。
1 function y = cauchy( theta )
2 %% Returns the unnormalized density of the Cauchy distribution
3 y = 1 ./ (1 + theta.ˆ2);
Listing 2.1: Matlab function to evaluate the unnormalized Cauchy.
1 %% Chapter 2. Use Metropolis procedure to sample from Cauchy density
2
3 %% Initialize the Metropolis sampler
4 T = 500; % Set the maximum number of iterations
5 sigma = 1; % Set standard deviation of normal proposal density
6 thetamin = −30; thetamax = 30; % define a range for starting values
7 theta = zeros( 1 , T ); % Init storage space for our samples
8 seed=1; rand( 'state' , seed ); randn('state',seed ); % set the random seed
9 theta(1) = unifrnd( thetamin , thetamax ); % Generate start value
10
11 %% Start sampling
12 t = 1;
13 while t < T % Iterate until we have T samples
14 t = t + 1;
15 % Propose a new value for theta using a normal proposal density
16 theta star = normrnd( theta(t−1) , sigma );
17 % Calculate the acceptance ratio
18 alpha = min( [ 1 cauchy( theta star ) / cauchy( theta(t−1) ) ] );
19 % Draw a uniform deviate from [ 0 1 ]
20 u = rand;
21 % Do we accept this proposal?
22 if u < alpha
23 theta(t) = theta star; % If so, proposal becomes new state
24 else
25 theta(t) = theta(t−1); % If not, copy old state
26 end
27 end
28
29 %% Display histogram of our samples
30 figure( 1 ); clf;
31 subplot( 3,1,1 );
32 nbins = 200;
33 thetabins = linspace( thetamin , thetamax , nbins );
34 counts = hist( theta , thetabins );
35 bar( thetabins , counts/sum(counts) , 'k' );
36 xlim( [ thetamin thetamax ] );
37 xlabel( '\theta' ); ylabel( 'p(\theta)' );
38
39 %% Overlay the theoretical density
40 y