MCMC抽样-Metropolis算法

肖柳青上海交通大学数学系随机模拟方法与应用StochasticSimulationMethodsandItsApplications肖柳青博士[e-mail]:[e-mail]肖柳青上海交通大学数学系第6章神奇的马尔科夫链蒙特卡罗方法6.1马尔科夫链6.2MCMC抽样―Metropolis算法6.3几个MCMC的例子6.4为什么Metropolis算法能有效工作？肖柳青上海交通大学数学系36.1马尔可夫链与转移概率随机试验结果{Xt，tT}，参数T={1,2,},状态空间={x1,x2,}•定义:若随机过程{Xt，tT}，对任意tT和x1,,xn+1，条件概率P{Xn+1=ixn+1|X1=x1,X2=x2,,Xn=in}=P{Xn+1=xn+1|Xn=xn}，则称{Xn，nT}为马尔可夫链，简称马氏链。若t1,t2,,tn-2表示过去，tn-1表示现在，tn表示将来，马尔可夫过程表明：在已知现在状态的条件下，将来所处的状态与过去状态无关，所谓“遗忘性”。肖柳青上海交通大学数学系()Xtt01234k1k…………肖柳青上海交通大学数学系一些简单例子假设甲乙两人以抛硬币的方式进行赌博，每次抛同一枚硬币；若出现正面，则甲付给乙一元钱，若出现反面，则乙付给甲一元钱.记为第局之后甲赢的总数.则是马尔可夫链.()Xnn{()0}Xnn，（简单传染病模型）设个人中有部分人感染了某种传染病，假设：当一个病人接触了一个健康者时，健康者被感染的概率为；所有的接触都是两人之间的接触；任意两个人的接触都是等可能的；在一个单位时间内只有一次接触发生.记为时刻时患病的人数，则是马尔科夫链.n{()0}Xnn，p()Xnn肖柳青上海交通大学数学系马尔科夫链的模型广泛应用性•例如，一个随时间变化的热力学系统；•改变一个物种的DNA序列的突变，蛋白质分子一步一步折叠的序列；•一天又一天的股价波动；•一个赌徒的赌博资金等等。肖柳青上海交通大学数学系7马尔可夫过程通常分为三类：(1)时间、状态都是离散的，称为马尔可夫链(2)时间连续、状态离散的，称为连续时间马尔可夫链(3)时间、状态都是连续的，称为马尔可夫过程肖柳青上海交通大学数学系86.1.1马尔可夫链的表示方式•定义:称条件概率pij(n)=P{Xn+1=j|Xn=i}为马尔可夫链{Xn，nT}在时刻n的一步转移概率，简称转移概率，其中i,j。•定义:若对任意的i,j，马尔可夫链{Xn,nT}的转移概率pij(n)与n无关，则称马尔可夫链是齐次的，并记pij(n)为pij。•齐次马尔可夫链具有平稳转移概率，状态空间={1,2,3,}，一步转移概率为肖柳青上海交通大学数学系9一步转移概率矩阵：•转移概率性质(1)(2)111212122212nnnnnnppppppPppp0,,ijpij11nijjpnnP肖柳青上海交通大学数学系例如：肖柳青上海交通大学数学系•这个拥有4个状态的马尔科夫链，其转移情况由图5.1所示，则相应的转移矩阵为•我们设初始的状态是，那么下一步的状态将是哪个呢？1112142123313443440000000pppppPpppp1x肖柳青上海交通大学数学系•生成服从均匀分布一个随机数𝑟∼𝑈(0,1)，看𝑟落在[0,𝑝11]、(𝑝11,𝑝11+𝑝12]或(𝑝11+𝑝12,1]的哪个区间里？•则下一步状态就转移到它们分别对应的于状态𝑥1、𝑥2或𝑥4。•依次类推不断地进行转移，则产生了该马尔科夫链的一条状态演化路径。肖柳青上海交通大学数学系•很显然，状态的演化是随机的!•我们以向量p𝑡=(𝑝1,𝑝2,𝑝3,𝑝4)表示经过𝑡步转移后状态的概率分布（注：这里用行向量），这向量被称为概率向量。•概率向量就是初始时状态处于状态𝑥1的概率分布。•根据概率公式，我们可得概率分布的转移关系式：因此，马尔科夫链的状态演化完全由其转移概率所决定。22110tttppPpPpP0(1,0,0,0)p肖柳青上海交通大学数学系6.1.2一个股市的马氏链模型•考虑一个股票市场的模型。•因为根据金融市场有效性假设：股价的下一次变化只取决于当前状态，具有马尔科夫性，所以股价的变化可以用马尔科夫链来建模描述。肖柳青上海交通大学数学系•例6.2.众所周知，明天的股价相对于今天来说有“涨”、“持平”或者“跌”三种可能的状态，分别以数字1，2，3表示之；假设其状态转移的变化规律是：•如果股价今天是涨，那么明天再涨的概率是0.3，持平的概率是0.2，而跌的概率是0.5；•如果股价今天是持平，那么明天上涨的概率是0.4，再持平的概率是0.2，而跌的概率是0.4；•如果股价今天是跌，那么明天会涨的概率是0.4，持平的概率是0.3，而再跌的概率是0.3。肖柳青上海交通大学数学系•我们用图（见图6.2）来描述这个马尔科夫链，其中节点表示那三个状态：涨（up）、持平（same）和跌（down）；有向边和标出的数值表示状态转移方向及其概率，即转移概率矩阵为0.30.20.50.40.20.40.40.30.3P肖柳青上海交通大学数学系•模拟股市的变化状况:•首先需要给定一个起始状态:可随机选择或确定一个状态。•如果我们假定今天状态是“涨”，接下来我们需要生成服从于转移概率分布的随机变量𝑥，•为此，用均匀分布的随机数𝑟∼𝑈(0,1)来产生下面的随机数：111213{0.3,0.2,0.5}ppp1111111211121,02,3,1rpxprppppr肖柳青上海交通大学数学系•对于其它初始状态情形，譬如状态2，则相应的转移概率分布为：•类似地我们可以生成随机数。事实上，我们可以定义关于转移概率的累积分布矩阵：其中元素是累积概率。212223{0.4,0.2,0.4}ppp0.30.510.40.610.40.71C1jijikkcp肖柳青上海交通大学数学系下面是Matlab的模拟程序肖柳青上海交通大学数学系•我们选择初始状态𝑥0=1为例，程序运行结果以各状态的出现频率的直方图形式呈示。•存在不变分布，那也被称为“稳定分布”。程序给出的不变概率分布是：•结果：今后股价出现上涨的比例约为36%，而出现下跌的比例却要高于上涨的，约为40%。(0.36415,0.23950,0.39635)p肖柳青上海交通大学数学系这例子有两个问题值得注意：•1.人们如何获得转移概率？•2.为什么要假设该马尔科夫链是齐次的，即转移概率不随时间变化？肖柳青上海交通大学数学系6.1.3不变分布•上例让我们看到不变分布的出现，它是经过长时间迭代后呈现出的一种不变性，也就说，这是一个稳定的极限分布。•对于一般的马尔可夫链，有迭代关系式：•上式对𝑡取极限，假设概率向量的极限存在并它设为𝜋，则就有这个𝜋就被称为不变分布。1ttppPtpP肖柳青上海交通大学数学系定理6.1(Perron–Frobenius)•如果概率转移矩阵𝑃具有𝑝𝑖𝑗𝛿0,∀𝑖,𝑗，则1.𝑃存在有1的特征值，且对应的左特征向量𝑤严格为正的，它还是唯一的。•2.如果此特征向量被归一化，则进一步有：•3.特别地，对于任意概率向量𝛼有：•其中𝑤是归一的。这里1是分量均为1的列向量，左特征向量𝑤是行向量。•定理的证明此略，请参考有关的教材。lim1nnPwlimnnPw肖柳青上海交通大学数学系•显然，我们现在有两种方法来求不变分布：•一种是从任意概率向量出发迭代计算；•另一种是直接求转移概率矩阵的特征值为1的左特征向量（还要归一化）；•我们推荐前一种方法!因为这是一个迭代的极限过程，可借助Matlab软件轻松获得其近似解，人们不需要判别转移矩阵是否满足那些存在不变分布的条件。•不难求得其结果是：(0.3636,0.2397,0.3967)肖柳青上海交通大学数学系6.2MCMC抽样―Metropolis算法•马尔科夫链蒙特卡罗抽样方法（简记为MCMC）的想法可追溯到1953年N.Metropolis等人在研究关于原子和分子的随机性运动问题时所引入的随机模拟方法。该方法被命名为“Metropolis模拟算法”。•事实上，这个算法已被列为影响科学和工程技术的发展最伟大的十大算法之首。肖柳青上海交通大学数学系•Metropolis算法是MCMC的核心。•MCMC背后的基本思想是构造一个遍历的马尔可夫链，使得其不变分布成为人们所需要的抽样分布。做到这点似乎相当复杂，但实际上由于人们可以非常灵活地选择简单的转移概率，所以构造该算法并不困难。•Metropolis算法的动机是推广接受-拒绝抽样方法。回忆在接受-拒绝抽样方法中，我们有目标密度函数、建议密度函数和一个接受准则（概率函数）。同样地我们假设：是全空间上的离散概率密度函数，人们需要在上按照密度函数产生样本马尔可夫链.肖柳青上海交通大学数学系•类似接受-拒绝抽样方法,在当前状态下产生下一个状态将分两个步骤进行：•（1）从建议密度函数产生一个随机数作为建议的下一个状态；•（2）然后按均匀分布产生一个随机数，如果,则接受该建议的随机数，即，否则放弃𝑦而采用原状态•重复这个过程产生随机序列。很显然，这样•过程所产生的随机序列做到了下一随机数仅依赖于当前的随机数而和以前产生的随机数无关。tx(|)tgxyr(,)rhxy1txy1ttxx肖柳青上海交通大学数学系•一种特殊的Metropolis算法要求𝑔满足对称性，即满足条件：𝑔(𝑦|𝑥)=𝑔(𝑥|𝑦)•这种对称性条件在大多数情况下可以很自然地做到。当然它也可以稍微减弱成下面的条件：•𝑔(𝑦|𝑥)0⇔𝑔(𝑥|𝑦)0•即要求状态转移是可逆的。肖柳青上海交通大学数学系•算法接下来要做的是，将选中的状态和当前状态进行比较，以一定的概率接受其中之一为随机序列链的下一状态。这里与接受-拒绝抽样方法的另一不同之处在于接受概率不但取决于下一步状态𝑦而且还与当前状态𝑥有关。所谓的Metropolis-Hastings算法，该算法的接受概率定义为()(|)(,)min1,()(|)fygxyhxyfxgyx肖柳青上海交通大学数学系•当对称性条件满足时，上式简化成•这就是说，如果建议的下一状态𝑦具有比当前状态𝑥的概率较大，𝑓(𝑦)≥𝑓(𝑥)，则肯定接受它否则，接受概率是𝑓(𝑦)/𝑓(𝑥)。这里需要指出：MCMC算法只需知道𝑓的相对值，即只要给出一个正比于𝑓的函数即可，这种方便也是MCMC的优势之一，因为有些应用问题难以将𝑓归一。•虽然这个Metropolis-Hastings算法看起来很简单，但它却非常有用，连同所有改进的算法在一起它们已在许多学科领域里有着重要的应用。()(,)min1,()fyhxyfx肖柳青上海交通大学数学系肖柳青上海交通大学数学系•如果建议密度函数𝑔满足对称性条件，则概率函数ℎ简化成（5.8）式，此时算法被称为对称Metropolis-Hastings算法。•最后需指出，Metropolis-Hastings算法能够直接推广到可数多状态和连续状态空间上去，下面的例子将展示这种算法的具体运用方法。肖柳青上海交通大学数学系5.3几个MCMC的例子•例1:模拟掷双骰子的游戏。•它的状态空间为{2,3,···,12}，•表5.1尚未归一的目标概率密度函数𝑓。肖柳青上海交通大学数学系6.3.1极小邻域法•极小邻域法:每个状态构造的邻域是由与它最相邻的状态组成，即邻居只有一个或者两个。具体说，对于每个状态𝑥，定义建议密度函数为1/2,max{1,2}(|)1/2,min{1,12}0,yxgyxyx其它肖柳青上海交通大学数学系