MCMC

背景

给定一个的概率分布 $p(x)$, 我们希望产生服从该分布的样本。前面介绍过一些随机采样算法（如拒绝采样、重要性采样）可以产生服从特定分布的样本，但是这些采样算法存在一些缺陷（如难以选取合适的建议分布，只适合一元随机变量等）。下面将介绍一种更有效的随机变量采样方法：MCMC 和 Gibbs采样，这两种采样方法不仅效率更高，而且适用于多元随机变量的采样。

如果一定条件下，马尔可夫链可以收敛到平稳分布。一个绝妙的想法是：如果能构造一个转移矩阵为 $P$ 的马尔可夫链，是其平稳分布刚好是 $p(x)$ 。那我们就可以从任何一个初始状态 $x_{0}$ 出发，沿着马尔可夫链转移，得到一个转移序列 $x_{0}$ , $x_{1}$ , $x_{2}$ ,…, $x_{n}$ , $x_{n+1}$ , $x_{T}$ , 如果马尔可夫链在第 $n$ 布已经收敛，我们就得到了 $p(x)$ 的样本 $x_{n},x_{n+1} \ldots x_{T}$ , 过程如下所示：

1. 设置 $t=0$
2. 生成一个随机状态 $u$， 设定初始状态 $x_{0} = u$
3. 重复
   $t=t+1$
   根据转移概率 $p(x_{t}|x_{t-1})$ , 采样得到 $u$
   设置 $x_{t} = u$
4. 直到 $t=T$

这种方法在1953年被 Metropolis 想到了，为了研究粒子系统的平稳性质，Metropolis 考虑了常见的玻尔兹曼分布的采样问题，首次提出了基于马尔可夫链蒙特卡罗的采样方法，即 Metropolis 算法，并在计算机上进行了编程实现。Metropolis 算法是首个普适的采样方法，并启发了一系列MCMC方法，所以人们把它作为随机模拟技术腾飞的起点。Metropolis 的这篇论文被收录在《统计学中的重大突破》，Metropolis 算法当选为二十世纪十个最重要的算法之一。

Metropolis 算法

假定目标概率分布为 $p(x)$ , 我们已经有一个转移矩阵为 $Q$ （$q(i,j)$ 表示从状态 $i$ 转移到状态 $j$ 的概率）的马尔可夫链。通常情况下不满足马尔可夫链的细致平稳条件：

$$p(i)q(i,j) \not = p(j)q(j,i)$$
所以 $p(x)$ 通常不是这个马尔可夫链的平稳分布。我们可否对该马氏链做一个改造，使得细致平稳条件成立呢？譬如，我们引入一个接受率 $\alpha(i,j)$, 我们希望：
$$p(i)q(i,j)\alpha(i,j)=p(j)q(j,i)a(j,i)$$
取什么样的 $\alpha(i,j)$ 才能让上式成立呢？最简单的，根据对称性，我们可以取
$$\alpha(i,j)=p(j)q(j,i)$$ $$\alpha(j,i)=p(i)q(i,j)$$
此时细致平稳条件成立：
$$p(i)\underbrace{q(i,j)\alpha(i,j)}_{q'(i,j)}=p(j)\underbrace{q(j,i)\alpha(j,i)}_{q'(j,i)}$$
于是我们将原来具有转移矩阵 $Q$ 的一个普通的马氏链，改造成了具有转移矩阵 $Q'$ （ $q'(i,j)$ 表示从状态 $i$ 转移到状态 $j$ 的概率）的马氏链。而此马氏链的平稳分布刚好就是目标概率分布 $p(x)$ ！

接受率 $\alpha(i,j)$ 的理解

由于 $\alpha(i,j) = p(j)q(j,i) \leq 1$ ，因此

$$\sum_{j \in S} q'(i,j) = \sum_{j \in S} q(i,j) \alpha(i,j) \leq \sum_{j \in S} q(i,j) =1$$ 索引矩阵 $Q'$ 不满足马尔可夫链状态转移矩阵的条件，因而不能直接作为马尔可夫链的转移矩阵。

细致平稳条件的直观含义如下图所示：

一般介绍MCMC的资料中没有说明状态自我转移永远是满足细致平稳条件的，自己转给自己，当然转入和转出相等。 因此可以增加自我转移的概率 $q'(i,i)$ ，使 $\sum_{j \in S} q'(i,j)=1$ 成立 。Metropolis算法细致平稳条件如下图所示：

从以上分析可以发现，$\alpha(i,j)$ 表示在原来转移矩阵为 $Q$ 的马氏链中，从状态 $i$ 以概率 $q(i,j)$ 转移到状态 $j$ 的时候，我们以 $\alpha(i,j)$ 接受这个转移。而以概率 $q(i,j)(1-\alpha(i,j))$ 拒绝这个转移（进行自我转移）。因此称 $\alpha$ 为接受率。

Metropolis算法流程

假定目标分布为 $p(x)$ ，我们已经有了一个转移矩阵为 $Q$ （对应元素为 $q(i,j)$ )，经典MCMC算法流程如下所示

以上分析中 $p(x)$ , $q(x|y)$ 说的都是离散的情形，事实上即便这两个分布式连续的，以上算法仍然有效（此时 $p(x)$ 表示概率密度，$q(x|y)$ 表示条件概率密度），因此可以生成服从连续的概率分布 $p(x)$ 的样本。

注：MCMC算法是以于马尔可夫链为基础的。马尔可夫链通常研究的是时间和状态都是离散的随机过程。而连续分布实际上对应的状态是连续的。那么MCMC算法可以对连续分布进行采样的原因又是什么？？

Metropolis-Hastings 算法

以上的Metropolis采样算法已经能漂亮地工作了，不过它有个小问题：马氏链在转移的过程中的接受率 $\alpha(i,j)$ 可能偏小，采样过程中容易原地踏步（自我转移），拒绝大量的跳转。这会导致马氏链收敛到稳态分布需要很长的时间。有没有办法提升接受率呢？

假设在Metropolis算法中 $\alpha(i,j)=0.1$ ，$\alpha(j,i)=0.2$，满足细致平稳条件，下式成立：

$$p(i)q(i,j) \times 0.1=p(j)q(j,i) \times 0.2$$
将上式两边同时扩大5倍
$$p(i)q(i,j) \times 0.5=p(j)q(j,i) \times 1$$
看！我们将接受率 $\alpha(i,j)$ ， $\alpha(j,i)$ 放大了5倍，而细致平稳条件仍然成立！这启发 我们可以把接受率 $\alpha(i,j)$ ，$\alpha(j,i)$ 同比例放大，最优的情况是使得两个数中较大的一个放大到1（因为接受率要求小于等于1）。所以我们可以取
$$\alpha(i,j)=\min\{\frac{p(j)q(j,i)}{p(i)q(i,j)},1\}$$
于是，通过对Metropolis算法中的接受率的微小改造，我们就得到了如下教科书中最常见的Metropolis-Hastings 算法。

注：对未归一化的概率分布 $p(x) = \frac{\tilde p(x)}{X_{p}}$ ($X_{p}$ 未知)。MH算法仍然适用。只需将 Algorithm 6 中的接受率可变换为：

α(xt,y)=min{p(y)q(xt|y)p(xt)p(y|xt),1}=min{p~(y)q(xt|y)p~(xt)p(y|xt),1}

$$\begin{equation} \begin{split} \alpha(x_{t},y) &= min\{\frac{p(y)q(x_{t}|y)}{p(x_{t})p(y|x_{t})},1 \} \\ &= min\{\frac{\tilde p(y)q(x_{t}|y)}{\tilde p(x_{t})p(y|x_{t})},1 \} \end{split} \end{equation}$$ 可以发现接受率 $\alpha$ 与归一化常数无关！

例子

采用 MH算法生成服从Cauchy分布的样本。MH算法采样的过程如下所示：

当前的状态为 $x_{t}$

利用建议分布（这里使用正态分布）进行采样得到 $y$

若接受转移，则 $x_{t+1}=y$

matlab代码如下：

%%M-H算法
T = 500;  %the maximum number of iterations
sigma = 1;%the deviation of normal proposal density
x_min=-10; x_max=10   %define a range for starting values
x = zeros(1,T);%init storage space for our samples
seed=1; rand('state',seed); randn('state',seed);    %random seed
x(1) = unifrnd(x_min,x_max);

%%Start sampling
t=1;
while t<T   %Iterate until we have T samples
    t=t+1;
    %propose a new value for theta using a normal proposal density
    y = normrnd(x(t-1),sigma);
    %Calculate the acceptance ratio
    alpha = min([1,(cauchy(y)*normpdf(x(t-1),y,sigma))/(cauchy(x(t-1))*normpdf(y,x(t-1),sigma))]);
    %draw a uniform deviate from [0,1]
    u=rand;
    %accept this proposal?
    if u<alpha
        x(t)=y;
    else
        x(t)=x(t-1);
    end
end

%%  Display histogram of our samples
figure(1); clf;
subplot(3,1,1);
nbins=200;
thetabins = linspace(x_min,x_max,nbins);
counts = hist(x,thetabins);
bar(thetabins,counts/sum(counts),'k');
xlim([x_min x_max]);
xlabel('x'); ylabel('p(x)');

%% Overlay the theoretical density
y=cauchy(thetabins);
hold on;
plot(thetabins,y/sum(y),'r--','LineWidth',3);
set(gca,'YTick',[]);

%% Display history of our samples
subplot(3,1,2:3);
stairs(x,1:T,'k-');
ylabel('t');xlabel('x');
set(gca,'YDir','reverse');
xlim([x_min x_max]);

实验结果：

参考文献

随机采样方法整理与讲解（MCMC、Gibbs Sampling等
MCMC方法的理解
Computational Statistics with Matlab