MCMC方法ppt课件

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