MCMC方法

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,MCMC,方法,Markov,链在统计计算中的应用,随机过程,.,.,Markov,链,下面我们进行一个独立重复掷色子实验,假设掷得,1,点的概率为,p,1,(0,p,1,1) ,设,X,n,表示投掷,n,次后掷得,1,点的累计数目。显然，,X,n,之间并不相互独立，但是，若给定,X,n,的值，如,X,n,=,i,，则,X,n+,1,的值只能取,i,或,i,+1,，,对应的概率分别是,1-,p,1,和,p,1,。,Markov,性,:,若已知现在的状态，将来与过去无关。,具有,Markov,性的离散时间随机过程称为,Markov,链。,离散时间,Markov,链,定义,：考虑只取有限个或可数个状态的随机过程，有,上式右边的条件概率称为,Markov,链在时刻,m,处于状态,i,条件下，在时刻,m+n,转移到状态,j,的（,n,步）转移概率，记为,p,ij,(,n,),。,由于链在时刻,m,从任何一个状态,i,出发，到另一时刻,m+n,，必然转移,j,1,，,j,2,，,诸状态中的某一个，所以,由转移概率组成的矩阵为马氏链的转移概率矩阵。由上式知，此矩阵的每一行元素之和等于,1,。,Chapman-,Kolmogorov,方程,C-K,方程,:,利用,C-K,方程我们容易确定,n,步转移概率,事实上，在上式中令,m,1,，,n,n-1,，得递推关系：,P(n,)=P(1)P(n-1)=PP(n-1)=,P,n,就是说，对齐次马氏链而言，,n,步转移概率矩阵是一步转移概率矩阵的,n,次方。,进而可知，一个,Markov,链的概率分布完全由它的一步概率矩阵与初始分布决定。,遗传学中的,Hardy-Weinberg,平衡定律,考虑一个生物群体，其中每一个个体带有一特定的基因型，假定有,A,和,a,两个等位基因，从而基因型有三种可能组合：,AA,，,Aa,，,aa,。,设群体在开始观测时（第,0,代）雄性与雌性具有相同的基因型频率分布：,d,:,2h,:,r,，则,A,和,a,的基因频率分别为,p,=,d,+,h,q,=,r,+,h,。,下面用,Markov,链来描述遗传过程，用,1,，,2,，,3,分别表示,AA,，,Aa,，,aa,，用,p,ij,表示给定一个亲代（父或母）的基因型,i,时子代出现基因型,j,的概率，在随机交配的前提下,p,12,=P(,子代基因型为,Aa,|,母亲基因型为,AA),=,P(,父亲基因型含,a,|,母亲基因型为,AA),=,P(,父亲基因型含,a,)=,q,类似地,可算得其他一步转移概率，一步转移概率矩阵为,第,0,代具有基因型分布,(,d,2h,r),从而第,1,代基因型分布为,(,d,2h,r)P=(p,2, 2pq, q,2,),第,2,代基因型分布为,(p,2, 2pq, q,2,) P=(p,2, 2pq, q,2,),由此可见，第,2,代的基因型分布与第,1,代相同，都为。按类似计算，可知第,3,代、第,4,代,的基因型分布仍为,(p,2, 2pq, q,2,).,Hardy-Weinberg,平衡定律，即不论父母基因型频率是什么数值，在随机交配的假定下，第,1,代继承者将有基因型频率,(p,2, 2pq, q,2,),，,而且这样的频率将保持永远恒定。,平稳分布,若存在一概率分布,j,使得转移概率矩阵满足,则称此,Markov,链具有平稳分布,j,在遗传学中的,Hardy-Weinberg,平稳定律一例中，,Markov,链具有平稳分布,(p,2, 2pq, q,2,),。,对于马尔可夫链的最终状态平稳分布常常是存在的，但是不一定是唯一的。我们利用马尔可夫链希望产生唯一的平稳分布。这需要产生具有遍历性的马尔可夫链。遍历性由非周期性和不可约性两个条件决定。我们下面对这三个概念做一个简要的介绍。,假定我们用如下方法产生一个随机变量序列,X,(0), X,(1),在任一时刻,t,，,序列中下一时刻,t,十,1,处的,X,(t+1),由条件分布,P(,x,|X,(t),),产生，例如从,N(0.5X,(t),),中抽取,X,(t+1),，,这样一个随机变量序列为,Markov,链,从,N(0.5X,(t),),中抽取,X,(t+1),，,X,(0),对,X,(t),的,影响,不管链的初始分布如何，经过长时间的转移后，链的分布收敛到同一分布（极限分布）。,Markov,链的,转移核,p,(,x,),对于离散状态,Markov,链，,有一步转移概率,p,(,x,x,),p,(,x, x,)=P(,X,(t+1),=,x,|,X,(t),=,x,),p,(,x, x,),作为,x,的函数可看成一离散型概率分布列。,连续状态,Markov,链，有一步转移概率密度,p,(,x,x,),p,(,x,B)=P(,X,(t+1),B,|,X,(t),=,x,)=,p,(,x,x,),作为,x,的函数可看成一连续型概率密度函数。,MCMC,方法,-,Markov,链在统计计算中的应用,在统计计算中，我们经常需要计算某函数关于一概率分布的期望,(,如均数、方差等）,其中 为,k,维向量。,较简单时：直接计算、数值积分、静态,Monte Carlo,。,较复杂时：,Markov chain Monte Carlo (MCMC),方法。,MCMC,方法实施步骤：,MCMC,方法常按如下步骤实施：,1.,构造一个马氏链，使它具有 极限分布 ；,2.,从某个 出发，用上述马氏链产生点列,；,3.,选取适当的 ，计算 估计值,MCMC,方法在实施时，有时同时产生整个,是困难的，,(,降维的思想,),于是便要用到满条件分布，即形如 的条件分布,其中,注意到，,例,:,设,则满条件分布,即,类似可得 降维的思想。,在构造马氏链时，转移核的构造很重要，不同的,MCMC,方法，往往也是转移核的构造方法不同。下面介绍两种常用的,MCMC,方法。,Gibbs,抽样,给定,x,-T,的条件下，如下产生,x,=(,x,T,x,-,T,):,(1),x,-,T,=,x,-T,;,(2),构造转移核,即由分布 产生,当只含有一个元素时称为单元素,Gibbs,抽样,此时从满条件分布中抽样成为单变量抽样,是最简单的,MCMC,。,单元素,Gibbs,抽样的步骤,在给出初始值 后，假定第,t,步迭代开始时的估计值为 ，则第,t,步迭代分为如下,k,步：,(1),由满条件分布 抽取, ,(i),由满条件分布 抽取, ,(k),由满条件分布 抽取,记,Gibbs,抽样收敛性的判断（,m,大小的确定）通常采取如下两种办法来判断：,方法,1,：从不同的初始点出发，同时产生多个,Markov,链，把每条链的演变情况作成散点图，从而可直观地观察链何时稳定下来，与初始点情况基本无关。,方法,2,：对抽样形成的,Markov,链每隔一段取一个样本，得到一个样本子列，然后对此子列计算参数的样本均数，当这样得到的均值稳定后，便可认为收敛了。,如用,方法,1,同时产生,9,条链，其中一个参数的实现值的散点图：,例：,Gelfand,和,Smith,考虑如下的例，设某实验可能有五个结果，其出现的概率分别为,现有,22,次结果的观测值,y=(y,1,y,2,y,3,y,4,y,5,)=(14,1,1,1,5).,请估计参数。,不防将,y,1,与,y,4,都分割成两部分,y,1,z,1,+(y,1,-z,1,), y,4,=z,2,+(y,4,-z,2,),则似然函数,若参数的先验分布为：,则后验分布,从而满条件分布,由此进行,Gibbs,抽样。,Metropolis-Hastings,方法,Metropolis-Hastings,方法转移核的构造如下：,潜在的转移核 作为 的函数是一个概率密度或概率分布，称为建议分布,.,建议分布可以取各种形式，常把它取为易于抽取的分布。,称为接受概率。要使平稳分布为，接受概率最常用的一个选择是,如果链在时刻,t,处于状态,x,，,则首先由,q,(,x,x,),产生一个潜在的转移,x,x,，,然后以概率,(,x,x,),接受转移到,x,，,即作为链在下一时刻,t,+1,时的状态值；而以概率,1-,(,x,x,),拒绝转移到,x,，,即仍为链在下一时刻,t,+1,时的状态值。因此，具体实施：,产生一个区间,0,，,1,上均匀分布的随机数,u,令,相应地有基于满条件分布情形。,固定 ，由转移核产生一个潜在的转移 ，接受概率,Gibbs,抽样是一种特殊的,Metropolis-Hastings,方法，其中建议分布为 ，接受概率 为常数,1,。,在,Gibbs,抽样中 可能很难抽取，,而,Metropolis,方法具有很大的灵活性，它,可取 为易于抽取的分布。,