MCMC方法

MCMC方法Markov链在统计计算中的应用随机过程见仅与一概率及其规律。由此可间内接到的病人个数的的一段时，率，而是要研究诊接到病人的个数的概内门研究的不仅是单位时间但在实际问题中所需要。表示，且有离散随机变量，可用一条件下接到的病人个数一门诊在正常条件工作间内书中曾看到：在单位时我们在概率论中或统计个数的问题。人工作条件下，来到的病试研究某一门诊在正常例：t00,,2,1,0,!kekkXPXk称为随机过程。关的随机变量有量有关，而且还与另一参但与这种不随机变量就应表达为，机现象。在这种情况下变量，才能表达这类随有关的随机与时间参量机现象，就必须用一个现象了。要反映这类随点就不能反映这类随机无关的随机间无关的随机变量或与时时间tXTtTttX,.,,,ttt集合。线的交点的纵坐标值的与所有的这类阶梯曲它的取值范围就是直线表示了一随机变量，面，对每一固定的或函数）。另一方表示了一族阶梯曲线（之，即另一条阶梯曲线。总示另一次试验的结果，就表，中取与上图不同的若在。它的取值可用图表出（即一次试验），对固定的则。例如若令有关的一个数值的集合值，而是还与有关的数就不再是仅与可见，这种试验的结果此到的病人数的试验。由对门诊做一次观察其接的出现就表示确定后对每一在本例中，在000,,,,,1,,0,tttttXttXtXtXtXtX时所处的状态。为随机过程于时，称＝的状态空间。当，称为随机过程的取值范围记为集。数集的子称为参数空间，它是实表示时间，常中参数称为参数，在实际问题。其中或，；简记为为一随机过程。记为的随机变量族参数与之对应，则称依赖于上的一个随机变量，，均有定义在，若对每一数集及一参，定义：设已给概率空间0001,,E,TttX,,,,,t,,,TtRT,txEtxtXtXTttXTttXtXtXPFPF道或现实。的轨为随机过程对应于函数，有时也称的一个样本对应于为随机过程称，＝记为：或复值函数），此函数有随机性）函数（实值上的普通确定性（不具是定义在，（即对每完成的试验）一确定的上的随机变量，对于每，是定义在，对于每一确定的的二元函数中的数和参数集中的元空间，随机过程是概率有随机过程的定义可知0000000000,,,,,T,,,,tTtxtXtxtxtXtXPFtXTtMarkov链•下面我们进行一个独立重复掷色子实验,假设掷得1点的概率为p1(0p11),设Xn表示投掷n次后掷得1点的累计数目。显然，Xn之间并不相互独立，但是，若给定Xn的值，如Xn=i，则Xn+1的值只能取i或i+1，对应的概率分别是1-p1和p1。•Markov性:若已知现在的状态，将来与过去无关。•具有Markov性的离散时间随机过程称为Markov链。离散时间Markov链•定义：考虑只取有限个或可数个状态的随机过程，有上式右边的条件概率称为Markov链在时刻m处于状态i条件下，在时刻m+n转移到状态j的（n步）转移概率，记为pij(n)。由于链在时刻m从任何一个状态i出发，到另一时刻m+n，必然转移j1，j2，…诸状态中的某一个，所以由转移概率组成的矩阵为马氏链的转移概率矩阵。由上式知，此矩阵的每一行元素之和等于1。}|{},,,|{1111iXjXPiXiXiXjXPmnmmmmnm.,2,1,1,jnmmPjij见黑板）的一步转移概率矩阵。或由它们组成＝＝率概论特别重要的一步转移次马氏链。我们以下讨齐时齐的。我们限于讨论同时也称此是齐次的或稳性。时，称转移概率具有平＝即有关时，及时间只与当转移概率(1,,,1imjmijijijijijaXaXPPPnPnmmPnjinmmPChapman-Kolmogorov方程•C-K方程:•利用C-K方程我们容易确定n步转移概率•事实上，在上式中令m＝1，n＝n-1，得递推关系：•P(n)=P(1)P(n-1)=PP(n-1)=Pn•就是说，对齐次马氏链而言，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，用pij表示给定一个亲代（父或母）的基因型i时子代出现基因型j的概率，在随机交配的前提下•p12=P(子代基因型为Aa|母亲基因型为AA)=P(父亲基因型含a|母亲基因型为AA)=P(父亲基因型含a)=q•类似地,可算得其他一步转移概率，一步转移概率矩阵为•第0代具有基因型分布(d,2h,r)•从而第1代基因型分布为(d,2h,r)P=(p2,2pq,q2)•第2代基因型分布为(p2,2pq,q2)P=(p2,2pq,q2)•由此可见，第2代的基因型分布与第1代相同，都为。按类似计算，可知第3代、第4代……的基因型分布仍为(p2,2pq,q2).•Hardy-Weinberg平衡定律，即不论父母基因型频率是什么数值，在随机交配的假定下，第1代继承者将有基因型频率(p2,2pq,q2)，而且这样的频率将保持永远恒定。平稳分布•若存在一概率分布{j}使得转移概率矩阵满足•则称此Markov链具有平稳分布{j}•在遗传学中的Hardy-Weinberg平稳定律一例中，Markov链具有平稳分布(p2,2pq,q2)。•对于马尔可夫链的最终状态平稳分布常常是存在的，但是不一定是唯一的。我们利用马尔可夫链希望产生唯一的平稳分布。这需要产生具有遍历性的马尔可夫链。遍历性由非周期性和不可约性两个条件决定。我们下面对这三个概念做一个简要的介绍。PiPlimP,S.30,.2i,1i0:1,1.ijjjnijijijjinijnijnnSaapnjipnSiijiijj＝或不依赖于极限存在，转移概率如果对于所有的，马氏链的状态空间为遍历性：一般，设齐次的。。即所有的状态是互通得使都存在正整数态不可约性：对于任意状期的。是非周则称状态的最大公因子是对于每个状态若正整数集非周期性：对于链的平稳分布。马氏链的极限分布又是为链的极限分布。，，＝，则同时称＝又若则陈此链具有遍历性。或121j212121jjjjnnPnP•假定我们用如下方法产生一个随机变量序列{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的函数可看成一连续型概率密度函数。xdxxpB),(MCMC方法-----Markov链在统计计算中的应用在统计计算中，我们经常需要计算某函数关于一概率分布的期望(如均数、方差等）其中为k维向量。较简单时：直接计算、数值积分、静态MonteCarlo。较复杂时：MarkovchainMonteCarlo(MCMC)方法。dxxxffE)()(),,,(21kxxxx)(xMCMC方法实施步骤：•MCMC方法常按如下步骤实施：1.构造一个马氏链，使它具有极限分布；2.从某个出发，用上述马氏链产生点列；3.选取适当的，计算估计值)(x)0(X)()2()1(,,,nXXX)(,nmnmfEnmttXfmnfE1)()(1ˆ•MCMC方法在实施时，有时同时产生整个是困难的，(降维的思想)于是便要用到满条件分布，即形如的条件分布,其中注意到，),,,(21kXXXX)|(TTxx},,,2,1{kT},{},,{TixxTixxiTiTTTTdxxxxx)()()|()(x常数。差一个比例以及满条件分布可以相时，，在应用方法的一个显著优点是法计算，而往无验分布的正则化常数往积项。同时，复杂的后些乘分布密度函数通常是一上式很重要，因为后验有关的项需保留。的乘积项中，只有与即，在xMCMCMCMCxxT•例:设则满条件分布即类似可得降维的思想。•在构造马氏链时，转移核的构造很重要，不同的MCMC方法，往往也是转移核的构造方法不同。下面介绍两种常用的MCMC方法。),(21xx})1()1(21exp{2221xx)|(21xx})1()1(21exp{2221xx))1(,1(22xN)|(12xxGibbs抽样给定x-T的条件下，如下产生x=(xT,x-T):(1)x-T=x-T;(2)构造转移核即由分布产生•当只含有一个元素时称为单元素Gibbs抽样,此时从满条件分布中抽样成为单变量抽样,是最简单的MCMC。)|(ˆ)|,(TTTTTxxxxxp)|(TTxxTx单元素Gibbs抽样的步骤•在给出初始值后，假定第t步迭代开始时的估计值为，则第t步迭代分为如下k步：(1)由满条件分布抽取…………(i)由满条件分布抽取…………(k)由满条件分布抽取记),,,()0()0(2)0(1)0(kxxxx)1(tx),,|()1()1(21tktxxx)(1tx),,,,,|()1()1(1)(1)(1tktititixxxxx)(tix),,|()(1)(1tktkxxx)(tkx),,,()()(2)(1)(tktttxxxx•Gibbs抽样收敛性的判断（m大小的确定）通常采取如下两种办法来判断：•方法1：从不同的初始点出发，同时产生多个Markov链，把每条链的演变情况作成散点图，从而可直观地观察链何时稳定下来，与初始点情况基本无关。•方法2：对抽样形成的Markov链每隔一段取一个样本，得到一个样本子列，然后对此子列计算参数的样本均数，当这样得到的均值稳定后，便可认为收敛了。•如用方法1同时产生9条链，其中一个参数的实现值的散点图：•例：Gelfand和Smith考虑如下的例，设某实验可能有五个结果，其出现的概率分别为现有22次结果的观测值y=(y1,y2,y3,y4,y5)=(14,1,1,1,5).请估计参数。不防将y1与y4都分割成两部分y1＝z1+(y1-z1),y4=z2+(y4-z2),则似然函数)1(21,834,4,4,81452423112121834814),|,(yzyzyzyyzzyL若参数的先验分布为：则后验分布从而满条件分布52321)1(),|,(yzyyzzy