MCMC抽样---Metropolis算法

MCMC抽样—Metropolis算法

马尔科夫链蒙特卡洛抽样方法可追溯到1953年N.Metropolis等人在研究原子和分子的随机性运动问题时所引入的随机模拟方法。该方法被命名为Metropolis模拟算法，这个算法已被列为影响科学和工程技术发展的最伟大的十大算法之首。

Metropolis算法是MCMC的核心。MCMC的基本思想是构造一个遍历的马尔科夫链，使得其不变分布成为人们所需要的抽样分布。做到这点似乎相当复杂，但实际上由于人们可以非常灵活地选择简单的转移概率，所以构造该算法并不困难。

Metropolis算法的动机是接受-拒绝抽样方法。在接受-拒绝抽样方法中，我们有布标概率密度函数f，建议概率密度函数g和一个接受准则h(x,y)。同样地我们假设：f是全空间$\Omega$上的目标概率密度函数，人们需要在$\Omega$上产生样本马尔科夫链$x_1,x_2,\cdots ,x_t,\cdots$，使得它的稳定分布恰好是这个目标概率密度函数f的分布。类似于接受-拒绝方法，在当前状态$x_t$下产生下一个状态分两个步骤进行：
（1）从建议概率密度函数$g(.|x_t)$产生一个随机数y作为建议的下一个状态（注意，这是依赖当前状态的条件概率）；
（2）然后按均匀分布生成一个随机数r，如果$r<h(x,y)$则接受该建议的随机数，即$x_{t+1}=y$，否则放弃y而采用原状态$x_{t+1}=x_t$。重复这个过程产生随机序列。很显然，这个过程所产生的随机序列做到了下一随机数仅依赖于当前的随机数，而和以前产生的随机数无关。

显然，在算法中建议的概率密度函数g主要目的是用来产生状态转移，即为每个状态构造出一个领域，并使邻域中的某个邻居倍选中。要使得它产生的序列是马尔科夫链，这要求建议概率密度函数g必须在整个$\Omega$上都有定义，这通常是容易做到的。接受-拒绝抽样方法是最简单的情况：对所有的$x,g(.|x)=g(.)$，单MCMC推广了它，使它变成了一个条件概率密度函数。当然必须要求在x的某个邻域$N_x$里有$g(y|x)>0,\forall y\in N_x$，这样才能够使状态产生遍历性的转移。一种特殊的Metropolis算法需要g满足对称性，即满足条件$g(y|x)=g(x|y)$，这种对称性在大多数情况下可以很自然地做到，当然它也可以稍微减弱成下面的条件$g(y|x)>0⟺g(x|y)>0$，即要求状态转移是可逆的。

算法接下来要做的是，将选中的状态和当前状态进行比较，以一定的概率接受其中之一为随机序列链的下一状态。这里与接受-拒绝方法的另一不同之处在于，接受概率不但取决于下一状态y，而且还与当前状态x有关。具体的方法是所谓的Metropolis-Hastings算法，该算法的接受概率定义为
$$h(x,y)=min1,\frac{f(y)g(x|y)}{f(x)g(y|x)}$$

党对称性条件满足时，上式简化成
$$h(x,y)=min(1,\frac{f(y)}{f(x)})$$

这就是说，如果建议的下一状态y具有比当前状态x的概率较大，即$f(y)\ge f(x)$则肯定接受它；否则，接受概率是$\frac{f(y)}{f(x)}$。这里需要指出：MCMC算法只需知道f的相对值，即只要给出一个正比于f的函数即可，这种方便也是MCMC的优势之一，因为有些应用问题难以将f归一化。虽然这个Metropolis-Hastings算法看起来简单，但它却非常有用，连同所有改进的算法一起，它们已在许多学科领域里有着重要的应用。

下面给出Metropolis-Hastings算法的伪代码：

如果建议概率密度函数g满足对称性条件，则概率函数h简化为$h(x,y)=min(1,\frac{f(y)}{f(x)})$，此时算法被称为对称Metropolis-Hastings算法

MCMC算法举例

1. 模拟掷双骰子的游戏，它的状态空间为$\Omega =2,3,\cdots ,12$下面的表给出了未归一的目标概率密度函数f，为了比较采用不同的建议概率密度函数g对于MCMC的影响，我们用两种不同的方法来构造建议概率密度函数：极大邻域法和极小邻域法

x	2	3	4	5	6	7	8	9	10	11	12
f(x)	1	2	3	4	5	6	5	4	3	2	1

极小邻域法：这个方法为每个状态构造的邻域是由与它最相邻的状态组成，即邻域只有一个或者两个元素。具体说，对于每个状态x，定义建议概率密度函数为：
g ( y ∣ x ) = {