MCMC-蒙特卡洛算法

1.马尔可夫链

     有这么一种链，就是当前点的状态的概率只与前一个点的状态有关，这就是一阶马尔可夫链。多阶的就是将与前一个点改成与前k个点。

这里有一个很经典的列子：

就是将一天的天气分成三种状态：出太阳，下雨，阴天。

首先明确点是什么？点就是某一天，可以设为x,

然后就是一条链了，链就是..., 昨天，今天，明天，后天,...

今天的天气状态和今天以前的我们都知道的，那么明天呢？我们现在要做一个天气预报，预报明天或后天的天气情况该怎么办呢？当然，这可以观测天象，但是，观察天象后也还是不能确切的知道什么情况明天会发生。所以，这里就有一个概率计算在里面。现在，假设连续天的天气情况是一个马尔可夫链，且是一阶的，即我们知道今天的情况，那么明天的情况就会由今天的情况所决定。记链为：x（t+1），t+1表示是第t+1天。

假设，今天的天气情况为出太阳。则可以将三种状态今天的概率表示为：(1,0,0)，即出太阳概率为1，其它情况为0.这是因为我们知道了今天是出太阳的。那么，明天天气是什么情况呢？这就需要知道转移概率了，所谓转移概率就是一个条件概率：P(明天出太阳/今天出太阳),P(明天下雨/今天出太阳)，P(明天阴天/今天出太阳)。是不是，这就直接的知道了明天各种天气情况的概率？当我们去推算后天的时候又该怎么样呢？明天具体是什么天气是不知道的，因此，只能用概率的形式去推算了。这样可以列举出下面的递推公式了：

π(t+1)，π（t）是我们每天天气的状况，这里t为时间，P为转移矩阵。对于三种天气状况，我们可得到3x3的转移矩阵，一个例子如下：

这里是借用了文献中的源数据：

这样的话，我们就可以一直的计算下去了，想要哪天的天气，一计算就可以了。但是在实际中是存在误差的，就是这个转移概率矩阵我们不能够精确的得到，所以，天气预报不敢报的很远了。

好了，这里来理一理这个算法的要素：1.假设一个链式马尔可夫链，即要满足马尔可夫的性质，当前状态只和前一状态有关；2.需要一个已知的转移概率矩阵。

以上，只是MCMC方法的基础。但有了这个以后就好办了。

2.MCMC生成随机样本

很多时候，我们都只是在做这样一件事情：就是我们的样本群都已经得到，然后去估计这个样本群服从什么分布。

那么，现在如果我们编写了一个算法，现在要去测试我们所编写的算法是否正确，那么我们是需要测试数据了。当我们的算法是去估算一个样本群的分布时，那么我们应该要去产生某分布的样本数据。这个想到该怎么产生吗？例如，现在要产生一个服从均匀分布的样本群，数量是在10000以上，如果想靠手动去产生，那是不可能的。所以得给用计算机了，这个要是在MATLAB中我们有现成的函数rand()来产生，其它方面也是能够有相应的函数调用。rand函数怎么实现就不说了。这里说一下，如果没有相应的函数来生成随机数该怎么办？例如，我们现在考虑二维高斯分布，那么现在要产生这个分布的数据该怎么产生呢？如果是要产生的数据有一百维呢，一千维呢，一万维呢（这个在图像中很常见的）？虽然知道了这些数据的具体分布，但是要怎么产生样本还是特别的困难。所以就有了MCMC方法了。

2.1Metropolis-Hastings algorithm 

首先这个算法的目的是：产生服从p(x)分布的样本数据，其中p(x) = f(x)/K.这里的K为一常数，但是不知道且很难计算，对于f(x)其表达式是已知的。

这里为Metropolis 算法的主要流程：

1.随机的选取初始样本x0.这里f(x0)>0,当然了，作为概率是应该要大于0的。

2.利用当前产生的x0去产生下一个样本x1,要得到下一个样本我们需要一个转移概率或是称作转移分布q(x0,x1)，这个分布可以自己选定。

3.对于给定的x1,计算α=p(x1)/p(x0) = f(x1)/f(x0)；对于Metropolis-Hasting算法，转移概率是需要加入到α的计算中去的，这样α表达式为：α=p(x1)*q(x0,x1)/(p(x0)*q(x1,x0) = f(x1)/f(x0),对于q如果是对称的，就转换成了Metropolis算法了。

4.决定x1是否接受：若α>1,则我们直接接受x1，若α<1,则我们就以α的概率接受这个x1.然后将x0= x1,返回到2继续。

下面是用MATLAB编写的算法：

%% 利用Metropolis-Hastings 算法产生样本
% 这里的样本分布为：p(x) = C* x.^(-n/2)*exp(-a/(2x));
% 对于参数n=5，参数a=4.
xlen = 200000;
x = zeros(1,xlen);
x0 = 1 ;
len = length(x);
k = 1 ;
while k <= len
    nextx = 100*rand();
    %compute the p from x0
    pn = (nextx^(-2.5))*(exp(-2/(nextx)));
    p0 = (x0^(-2.5))*(exp(-2/(x0)));
    p = min(pn/p0,1);
    if p >= 1
        x(k) = nextx ;
        x0 = nextx ;
        k = k + 1 ;    
    else
        pp = rand();
        if pp < p
            x(k) = nextx ;
            x0 = nextx   ;
            k = k + 1 ;      
        end
       
    end
   
end
hist = histc(x(10000:k-1),0:0.1:100);
plot(hist);  
xx = 0:0.1:100;
n = 5 ;
a = 4 ;
C = 1 ;
y = C .* xx.^(-n/2).*exp(-a./(2*xx));
figure;
plot(y);

当然，这里有个地方是不那么明确，就是关于q(x,y)该怎么选择？关于q(x,y)的选择主要有两种：1，随机漫步，即y的得到为y=x+z，其中z为一随机变量，可以为任何分布产生的随机变量，为了简单我们可以直接选取z为均匀分布，当然，均匀分布的区间要看x分布的而定。2，独立选取一个分布，即q(x,y)=g(y)，g(y)为我们选取的一个分布，当g(y)不是对称的时（g(y) != g(-y)）,在式中计算α时用Metropolis-Hastings算法，当然，可以选一个最简单的分布---均匀分布，这样g(x) = g(y) = C。

 

这里完成了一半，对于收敛度的证明不再讲述。文献参考为《Markov Chain Monte Carlo and Gibbs Sampling》Lecture Notes for EEB 581, version 26 April 2004°c B. Walsh 2004。