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，定义建议概率密度函数为：
$$g(y|x)=\{\begin{array}{rclc}1/2& ,& y=maxx-1,2,& \\ 1/2& ,& y=maxx+1,12,& \\ 0,& ,& others& \end{array}$$

注意到：g(3|2)=g(2|2)=1/2,g(11|12)=g(12|12)=1/2，这样的定义保证了对称性条件被满足，可以看到这样的建议矩阵是对称的，这样的状态转移方式明显是能够便利的，因为只要有足够多的迭代次数，它会在各状态之间不断游走，只有这样，MCMC采能够生成正确的数据，提供一个满足目标要求的不变分布。

close all
clear
f = [0, 1, 2, 3, 4, 5, 6, 5, 4, 3, 2, 1];
d = zeros(1, 20000);
x = 5;

for i = 1:20000
    u = rand;
    if x == 2
        if u < 0.5
            y = 3;
        else 
            y = 2;
        end
    elseif x == 12
        if u < 0.5
            y = 11;
        else
            y = 12;
        end
    else
        if u < 0.5
            y = x - 1;
        else
            y = x + 1;
        end
    end
    
    h = min(1, f(y) / f(x));
    u = rand;
    if u < h
        x = y;
    end
    d(i) = x;
end

a = 1:1:12;
hist(d, a);

我们看到，从mcmc模拟结果几乎精确地再现了目标概率分布。

极大邻域法：对所有的x定义其邻域为$N_x=\Omega$，也就是说，对于任意状态x，任何其他状态都可能被建议。如果以等概率取所有可能状态，这样的建议矩阵如下：
$$G=\left[\begin{array}{lccc}1/11& 1/11& \cdots & 1/11\\ 1/11& 1/11& \cdots & 1/11\\ ⋮& ⋮& & ⋮\\ 1/11& 1/11& \cdots & 1/11\end{array}\right]$$

这种邻域方式的建议概率密度函数g(.|x)是常量函数，它当然也满足对称性条件，这一，这种邻域选择方式将产生独立的抽样序列。

close all
clear
f = [0, 1, 2, 3, 4, 5, 6, 5, 4, 3, 2, 1];
d = zeros(1, 20000);
x = 5;

for i = 1:20000
    u = rand;
    if u < 1/11
        y = 2;
    elseif u < 2 / 11
        y = 3;
    elseif u < 3 / 11
        y = 4;
    elseif u < 4 / 11
        y = 5;
    elseif u < 5 / 11
        y = 6;
    elseif u < 6 / 11
        y = 7;
    elseif u < 7 / 11
        y = 8;
    elseif u < 8 / 11
        y = 9;
    elseif u < 9 / 11
        y = 10;
    elseif u < 10 / 11
        y = 11;
    else
        y = 12;
    end
    
    h = min(1, f(y) / f(x));
    u = rand;
    if u < h
        x = y;
    end
    d(i) = x;
end

a = 1:1:12;
hist(d, a);

我们看到，模拟结果同前面是一样的，这体现了MCMC的灵活性。

2. 连续马尔科夫链的模特卡罗模拟
   考虑密度函数为$f=0.5x^2{e}^{-x}$的伽马分布，我们要涉及的马尔科夫链是连续的状态空间$\Omega =[0,+\infty )$，这里要注意，目标概率密度函数有定义域的限制，也就是说，建议概率密度函数不能产生$y<0$的状态；但同时要使得建议概率密度函数是对称的，即满足$g(y|x)=g(x|y)$，我们定义这样的建议概率密度函数：对所有的x，取半径为1，中心在x的均匀分布密度，但在任何状态不能小于0.如果$y<0$，则令$y=x$，这样类似于上面离散例子的做法避免了不对称的情况

f = @(x)(0.5 * x * x * exp(-x));
d = zeros(1, 40000);
x = 2;
for i = 1:40000
    y = unifrnd(x-1, x+1);
    
    if y < 0
        y = x;
    end
    
    h = min(1, f(y) / f(x));
    u = rand;
    
    if u < h
        x = y;
    end
    
    d(i) = x;
end

a = 0:0.2:20;
hist(d(20001:40000), a);

MCMC模拟练习

1. 我们计算在举例1中的掷骰子问题在不同初始分布$p_0$下的概率分布$p_t$,例子中已经完成了对$p_0=(1,0,0)$情况的模拟，对同样的问题，做$p_0=(0,0,1)$情况下的模拟；对两种不同选择的情况下计算$p_t,t=1,2,\cdots ,10$，再利用这些计算，计算更一般情况$p_0=(a,b,1-a-b)$其中$a,b\ge 0,a+b<1$(事实上，矩阵乘法是线性的)。试随机模拟对于任何初始起点$p_0$，分布$p_t$收敛于不变分布。

p = [0.2, 0.3, 0.5];

P = [0.3, 0.2, 0.5; 0.4, 0.2, 0.4; 0.4, 0.3, 0.3];

n = 25000;


p_final = p;

for i = 1:n
    p_final = p_final * P;
    
end

p_final

将上述p配置成[1, 0, 0], [0, 1, 0], [0, 0, 1], [0.2, 0.3, 0.5]等满足题意的初始分布，经过计算后得到的结果均为：
p_final =

0.3636    0.2397    0.3967

故该马尔科夫链的最终均收敛于平稳状态，与初始分布无关；

2. 利用Metropolis算法从掷两个均匀流面骰子之和的样本空间$\Omega =2,3,\cdots ,12$取样，首先用$\Omega =2,3,\cdots ,12$上的均匀分布作为系统邻域，其次运用转移矩阵
   $$G=\left[\begin{array}{lcccccc}1/2& 1/2& 0& \cdots & 0& 0& 0\\ 1/2& 0& 1/2& \cdots & 0& 0& 0\\ 0& 1/2& 0& \cdots & 0& 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ 0& 0& 0& \cdots & 0& 1/2& 0\\ 0& 0& 0& \cdots & 1/2& 0& 1/2\\ 0& 0& 0& \cdots & 0& 1/2& 1/2\end{array}\right]$$

给出随机游走作为系统邻域。记录下计算时间，并画出直方图和相关图，发表评论你的意见。

仿真代码参考上面的掷骰子的最大邻域方法与最小邻域方法，经过计算，
最小邻域方法计算速度更快，最大邻域方法耗时较慢

3. 继续考察上例的取样问题，探讨MCMC方案的收敛速度。首先采用由上体转移矩阵所给出的随即又走的系统邻域分布。设样本序列是$X_1,X_2,\cdots ,X_N$，计算样本均值$\overline{X}=\frac{1}{N}{\sum }_{i=1}^N{X}_i$，运行N=100的100个十里，画出样本均值的直方图，并计算样本方差。再次重复运行N=200的100个实例，以及N=1000的100个实例。这给出了样本均值估计随迭代次数增加怎样收敛的启发。现在除了最大邻域系统外，所有状态都选择等可能，比较两个不同邻域系统在其收敛性分布直方图和方差图。

从理论层面分析$\overline{X}$样本均值的理想状态是这个分布的期望值，通过计算很容易得到其均值为7，随着N的增大，样本均值越来约接近其理论均值7，在实例维度，样本均值近似于一个正态分布，仿真代码如下：

close all
%clear

N = 10000;
M = 10000;
f = [0, 1, 2, 3, 4, 5, 6, 5, 4, 3, 2, 1];
d1 = zeros(M, N);
x = 5;
tic
for j = 1:M
    for i = 1:N
        u = rand;
    if u < 1/11
        y = 2;
    elseif u < 2 / 11
        y = 3;
    elseif u < 3 / 11
        y = 4;
    elseif u < 4 / 11
        y = 5;
    elseif u < 5 / 11
        y = 6;
    elseif u < 6 / 11
        y = 7;
    elseif u < 7 / 11
        y = 8;
    elseif u < 8 / 11
        y = 9;
    elseif u < 9 / 11
        y = 10;
    elseif u < 10 / 11
        y = 11;
    else
        y = 12;
    end
    
    h = min(1, f(y) / f(x));
    u = rand;
    if u < h
        x = y;
    end
    d1(j, i) = x;
    end
end
toc

d_sum = sum(d1, 2) / N;
fprintf('sequence_length = %d, std = %f\n', N, std(d_sum));
[dNum, dCen] = hist(d_sum, 20);
bar(dCen, dNum);

>> ex3_2
Elapsed time is 0.109177 seconds.
sequence_length = 100, std = 0.280242
>> ex3_2
Elapsed time is 0.726260 seconds.
sequence_length = 1000, std = 0.088843
>> ex3_2
Elapsed time is 7.015829 seconds.
sequence_length = 10000, std = 0.028191

从输出可知，随着序列长度的增加，$\overline{X}$的方差不断减小

4. 模拟马尔科夫链转移矩阵
   $$P=\left[\begin{array}{lcccc}& c_0& \frac{1}{2}{e}^{-3}& 0& 0\\ & \frac{1}{2}& 0& \frac{1}{2}& 0\\ & 0& \frac{1}{2}{e}^{-1}& c_2& \frac{1}{2}\\ & 0& 0& \frac{1}{2}{e}^{-1}& c_3\end{array}\right]$$

其中$c_i$是互补的概率，试模拟几百次迭代以获得不变分布并画出直方图，另外手工计算其不变分布，并比较计算一下结果：
$${\omega }_i=\frac{e^{-f_i}}{Z}$$

其中$Z={\sum }_{i=0}^3{e}^{-f_i},f_0=-1,f_1=2,f_3=0$

c0 = 1 - 0.5 * exp(-3);
c2 = 1 - 0.5 * exp(-1) - 0.5;
c3 = 1 - 0.5 * exp(-1);

P = [c0, 0.5*exp(-3), 0, 0; 0.5, 0, 0.5, 0; 0, 0.5*exp(-1), c2, 0.5; 0, 0, 0.5*exp(-1), c3];

p = [1, 0, 0, 0];

N = 10000;

p_final = zeros(1, 3);
p_final = p * P;


for i = 2:N
    p_final = p_final * P;
end

p_final

从1000次开始，分布进入不变状态，

>> ex5

p_final =

    0.6439    0.0321    0.0871    0.2369

>> ex5

p_final =

    0.6439    0.0321    0.0871    0.2369

>> 

通过手工计算
$${\omega }_0=0.0321,{\omega }_1=0.6439,{\omega }_2=0.2369,{\omega }_3=0.0871$$

和通过模拟迭代得到的不变分布一致。