MCMC: The Metropolis-Hastings Sampler

本文主要译自：MCMC:The Metropolis-Hastings Sampler

上一篇文章中，我们讨论了Metropolis 采样算法是如何利用马尔可夫链从一个复杂的，或未归一化的目标概率分布进行采样的。Metropolis 算法首先在马尔可夫链中基于上一个个状态 \(x^{(t-1)}\) 推荐一个新的状态 \(x^*\)，这个新状态是根部建议分布 \(q(x^*|x^{(t-1)})\) 进行采样得到的。算法基于目标分布函数在 \(x^*\) 上的取值接受或者拒绝 \(x^*\)。

Metropolis 采样方法的一个限制条件是推荐分布 \(q(x^*|x^{(t-1)})\) 必须是对称的。这个限制来源于使用马尔可夫链采样：从马尔可夫链的稳态分布进行采样的一个必要条件是在时刻 t，\(x^{(t-1)}\rightarrow x^{(t)}\) 的转移概率必须等于 \(x^{(t)}\rightarrow x^{(t-1)}\) 的转移概率，这个条件被称为可逆性或者细致平稳。然而，一个对称的分布或许对很多问题并不适合，比如我们想对定义在 \([0,\infty)\) 的分布进行采样。

为了能够使用非对称的推荐分布，Metropolis-Hastings 算法引入了一个附加的修正因子\(c\)，由推荐分布定义：

\(c = \frac{q(x^{(t-1)}|x^*)}{q(x^*|x^{(t-1)})}\)

修正因子调整了转移算子，从而保证了\(x^{(t-1)}\rightarrow x^{(t)}\) 的转移概率等于 \(x^{(t)}\rightarrow x^{(t-1)}\) 的转移概率，不管推荐分布是什么。

Metropolis-Hasting 算法的实现步骤与 Metropolis 采样完全相同，除了在计算接受概率 \(\alpha\) 时需要用到修正因子。为了得到 M 个采样点，Metropolis-Hasting 算法如下：

1. set t = 0

2. generate an initial state \(x^{(0)}\thicksim\pi^{(0)}\)

3. repeat until \(t=M\)

　　set \(t=t+1\)

　　generate a proposal state \(x^*\) from \(q(x|x^{(t-1)})\)

　　calculate the proposal correction factor \(c=\frac{q(x^{(t-1)}|x^*)}{q(x^*|x^{(t-1)})}\)

　　calculate the acceptance probability \(\alpha = \mathrm{min}\left(1, \frac{p(x^*)}{p(x^{(t-1)})}\times c\right)\)

　　draw a random number \(u\) from \(\mathrm{Unif}(0,1)\)

　　　　if \(u\leq\alpha\) accept the proposal state \(x^*\) and set \(x^{(t)}=x^*\)

　　　　else set \(x^{(t)}=x^{(t-1)}\)

很多人认为 Metropolis-Hastings 算法是 Metropolis 算法的一个推广。这是因为如果推荐分布是对称的，修正因子是1，就得到了 Metropolis 采样。

译者按：

这里文章只给出了算法，但是没有讲原理，本人大致说一下自己的理解： 

假如我们要用马尔可夫链对目标分布 \(\pi(x)\) 进行采样，我们需要保证马尔可夫链的稳态分布正是目标分布 \(\pi(x)\)。定义马尔可夫链的转移算子 \(T(a,b)\)，表示 \(a\rightarrow b\) 的转移概率。那么根据稳态分布的细致平稳条件，我们想达到这样的效果：

\(\pi(a)\cdot T(a,b) = \pi(b)\cdot T(b,a)\)

和 Metropolis 算法一样，我们引入推荐转移算子 \(Q(a,b)\) 并定义新的接受概率:

\(\alpha = \mathrm{min}\left(1, \frac{Q(b,a)}{Q(a,b)}\cdot\frac{\pi(b)}{\pi(a)}\right)\)

因此，\(a\rightarrow b\)  的转移概率可以这样计算：

\(\pi(a)\cdot T(a,b) = \pi(a)Q(a,b)\cdot\alpha\)

\(=\pi(a)Q(a,b)\cdot\mathrm{min}\left(1,\frac{Q(b,a)}{Q(a,b)}\cdot\frac{\pi(b)}{\pi(a)}\right)\)

\(=\mathrm{min}(\pi(a)Q(a,b), \pi(b)Q(b,a))\)

\(=\pi(b)Q(b,a)\cdot\mathrm{min}\left(\frac{\pi(a)Q(a,b)}{\pi(b)Q(b,a)},1\right)\)

\(=\pi(b)T(b,a)\)

果然，达到了细致平衡，这里面并不要求 \(Q(a,b) = Q(b,a)\) 也就是之前 Metropolis 算法中所要求的对称性条件。

上面的解释是说明了为啥这个算法可以得到细致平衡条件，下面我们正着想一下：

我们希望马尔可夫链的稳态分布是目标分布 \(\pi(x)\)，也就是说在 \(\pi(x)\) 时达到细致平衡。但是通常情况下，推荐转移算子 \(Q(a,b)\)不满足细致平衡（否则也就不用再计算接受率了）即：

\(\pi(a)\cdot Q(a,b) \neq \pi(b)\cdot Q(b,a)\)

于是，需要一个修正因子，也就是接受概率 \(\alpha\)，使得：

\(\pi(a)\cdot Q(a,b)\cdot\alpha(a,b) = \pi(b)\cdot Q(b,a)\cdot\alpha(b,a)\)

怎样取 \(\alpha\) 才能保证等式成立呢？最简单的取法是令：

\(\alpha(a,b) = \pi(b)\cdot Q(b,a)\), and

\(\alpha(b,a) = \pi(a)\cdot Q(a,b)\)

这里的 \(\alpha\) 之所以叫接受概率，是因为我们根据推荐分布 \(Q(a,b)\) 采样得到状态 b 后，有 \(\alpha(a,b)\) 的概率接受这个采样。这样，推荐分布和修正因子构成了收敛于目标分布的马尔可夫链。

但是，问题在于，\(\alpha(a,b)\) 有可能很小，这将导致马尔可夫链在采样过程中原地踏步，推荐分布推荐了一个，被拒绝了，又推荐了一个，还被拒绝了...马氏链每次采样都保留了之前的采样点。针对这种情况，我们可以对细致平衡等式两边的 \(\alpha(a,b)\) 和 \(\alpha(b,a)\) 同时放缩相同的倍数，使其中一个达到1，那么这就大大提高了接受率而不改变细致平衡等式。因此，我们的接受率 \(\tilde{\alpha}(a,b)\) 可以表示成：

\(\tilde{\alpha}(a,b) = \mathrm{min}\left(\frac{\pi(b)Q(b,a)}{\pi(a)Q(a,b)}, 1\right)\)

上面式子的含义是：

1. 如果 \(\pi(b)Q(b,a)>\pi(a)Q(a,b)\)， 那么以相同比例放缩 \(\alpha(a,b)\) 比 \(\alpha(b,a)\) 先达到1，\( \mathrm{min}\left(\frac{\pi(b)Q(b,a)}{\pi(a)Q(a,b)},1\right)=1\)，即，接受率为1，我们接受\(a\rightarrow b\) 的转移。

2. 如果 \(\pi(b)Q(b,a)<\pi(a)Q(a,b)\)，那么以相同比例放缩 \(\alpha(b,a)\) 比 \(\alpha(a,b)\) 先达到1，\( \mathrm{min}\left(\frac{\pi(b)Q(b,a)}{\pi(a)Q(a,b)},1\right)=\frac{\pi(b)Q(b,a)}{\pi(a)Q(a,b)}\)，接受率为 \(\frac{\pi(b)Q(b,a)}{\pi(a)Q(a,b)}\)，以这个接受率接受 \(a\rightarrow b\) 的转移。

经过这样的改造，我们就把 Metropolis 算法改造为了 Metropolis-Hastings 算法，从而不必寻找严格对称的推荐转移算子。

Example: Sampling from a Bayesian posterior with improper prior

(通过不合适的先验概率来对贝叶斯后验概率取样) 

 在很多应用中，包括回归和密度估计，通常我们需要确定一个假设模型 \(p(y|\theta)\) 的参数 \(\theta\)，使得这个模型最大程度地符合观测数据 \(y\)。模型函数 \(p(y|\theta)\) 通常被称为似然函数。在贝叶斯方法中，模型参数中通常有一个显式的先验分布 \(p(\theta)\) 来控制参数可能的取值。

模型的参数是由后验分布 \(p(\theta|y)\) 决定的，这是一个基于观测值的所有可能的参数概率分布。后验分布可以由贝叶斯定理确定：

\(p(\theta|y)=\frac{p(y|\theta)p(\theta)}{p(y)}\)

其中，\(p(y)\) 是一个归一化常数，通常很难直接求解。因为它需要计算参数和 \(y\) 所有可能的取值然后再对概率进行累加。

假如我们采用下面的模型（似然函数）：

\(p(y|\theta) = \mathrm{Gamma}(y;A,B)\)，其中

\(\mathrm{Gamma}(y;A,B)=\frac{B^A}{\Gamma(A)}y^{A-1}e^{-By}\)

\(\Gamma()\) 是伽马函数。因此，模型的参数为：

\(\theta = [A, B]\)

参数 A 控制分布的形状，参数 B 控制放缩。B=1, A从 0 遍历到 5 的似然面如下图所示。

Likelyhood surface and conditional probability p(y|A=2,B=1) in green

条件概率分布 \(p(y|A=2,B=1)\) 在似然面上被用绿色表示出来，在 matlab 中可以用下面的语句把它画出来：

plot(0:.1:10,gampdf(0:.1:10,2,1)); %GAMMA(2,1)

现在，我们假设模型的参数具有下面的先验概率：

\(p(B=1)=1\)

并且

\(p(A)=\mathrm{sin}(\pi A)^2\)　

第一个先验概率声明 B 只取一个值 (i.e. 1)，因此我们可以把它看作常数。第二个（非常规）先验概率声明，A 的概率变化符合正弦函数。（注意这两个先验概率分布都被称为不当先验(improper priors)，因为它们的积分值都不是1）。由于 B 是一个常数，我们只需要估计 A 的值。

事实证明，尽管归一化常数 \(p(y)\) 可能很难计算，我们仍然可以用 Metropolis-Hastings 算法从 \(p(A|y)\) 中取样，即使不知道 \(p(y)\)。尤其，我们可以忽略归一化常数 \(p(y)\) 直接从未归一化的后验概率抽样:

\(p(A|y)\propto p(y|A)p(A)\)

\(y\) 从 0 变化到 10 的后验分布平面如下图所示。图像的右侧蓝线表示参数 A 的先验概率 \(p(A)\)。假如我们有一个数据点 \(y=1.5\)，想通过 Metropolis-Hasting 算法估计后验分布 \(p(A|y=1.5)\)。下图中的品红色曲线表示这个特定的目标分布。

Posterior surface, prior distribution (blue), and target distribution (pink)

使用对称的建议分布，例如正态分布，从 \(p(A|y=1.5)\) 采样效率十分低，因为后验分布定义域为正实数。一个非对称的，相同定义域的推荐分布将会使得算法更好地收敛于后验分布。定义在正实数上的分布函数之一是指数分布。

\(q(A) = \mathrm{Exp}(\mu) = \mu e^{-\mu A}\),

这个分布含有一个参数 \(\mu\)，这个参数控制概率分布函数的位置和缩放。下图展示了目标后验分布和推荐分布（\(\mu=5\)）。

Target posterior p(A|y) and proposal distribution q(A)

我们发现推荐概率分布对后验分布有一个很好的覆盖。运行本文底部 Matlab 代码块的 Metropolis-Hastings 采样算法，得到下图中的马尔可夫链轨迹和采样结果：

Metropolis-Hastings Markov chain and samples

顺便说一句，注意到这个采样方法中，推荐分布函数不取决于上一个采样，而仅仅取决于参数 \(\mu\)（看下方的 Matlab 代码第88行）。每一个推荐状态 \(x^*\) 都是完全独立与上一个状态采集得到的。因此这是一个独立采样器 (independence sampler) 的例子，一种特殊的 Metropolis-Hastings 采样算法。独立采样器以其表现的极端性而闻名，效果要么很好，要么很差。效果的好坏取决于建议分布的选择以及建议分布的覆盖面。在实践中找到这样一个建议分布通常是困难的。

下面是 Metropolis-Hastings 采样器的代码

% METROPOLIS-HASTINGS BAYESIAN POSTERIOR
rand('seed',12345)
 
% PRIOR OVER SCALE PARAMETERS
B = 1;
 
% DEFINE LIKELIHOOD
likelihood = inline('(B.^A./gamma(A)).*y.^(A-1).*exp(-(B.*y))','y','A','B');
 
% CALCULATE AND VISUALIZE THE LIKELIHOOD SURFACE
yy = linspace(0,10,100);
AA = linspace(0.1,5,100);
likeSurf = zeros(numel(yy),numel(AA));
for iA = 1:numel(AA); likeSurf(:,iA)=likelihood(yy(:),AA(iA),B); end;
 
figure;
surf(likeSurf); ylabel('p(y|A)'); xlabel('A'); colormap hot
 
% DISPLAY CONDITIONAL AT A = 2
hold on; ly = plot3(ones(1,numel(AA))*40,1:100,likeSurf(:,40),'g','linewidth',3)
xlim([0 100]); ylim([0 100]);  axis normal
set(gca,'XTick',[0,100]); set(gca,'XTickLabel',[0 5]);
set(gca,'YTick',[0,100]); set(gca,'YTickLabel',[0 10]);
view(65,25)
legend(ly,'p(y|A=2)','Location','Northeast');
hold off;
title('p(y|A)');
 
% DEFINE PRIOR OVER SHAPE PARAMETERS
prior = inline('sin(pi*A).^2','A');
 
% DEFINE THE POSTERIOR
p = inline('(B.^A/gamma(A)).*y.^(A-1).*exp(-(B.*y)).*sin(pi*A).^2','y','A','B');
 
% CALCULATE AND DISPLAY THE POSTERIOR SURFACE
postSurf = zeros(size(likeSurf));
for iA = 1:numel(AA); postSurf(:,iA)=p(yy(:),AA(iA),B); end;
 
figure
surf(postSurf); ylabel('y'); xlabel('A'); colormap hot
 
% DISPLAY THE PRIOR
hold on; pA = plot3(1:100,ones(1,numel(AA))*100,prior(AA),'b','linewidth',3)
 
% SAMPLE FROM p(A | y = 1.5)
y = 1.5;
target = postSurf(16,:);
 
% DISPLAY POSTERIOR
psA = plot3(1:100, ones(1,numel(AA))*16,postSurf(16,:),'m','linewidth',3)
xlim([0 100]); ylim([0 100]);  axis normal
set(gca,'XTick',[0,100]); set(gca,'XTickLabel',[0 5]);
set(gca,'YTick',[0,100]); set(gca,'YTickLabel',[0 10]);
view(65,25)
legend([pA,psA],{'p(A)','p(A|y = 1.5)'},'Location','Northeast');
hold off
title('p(A|y)');
 
% INITIALIZE THE METROPOLIS-HASTINGS SAMPLER
% DEFINE PROPOSAL DENSITY
q = inline('exppdf(x,mu)','x','mu');
 
% MEAN FOR PROPOSAL DENSITY
mu = 5;
 
% DISPLAY TARGET AND PROPOSAL
figure; hold on;
th = plot(AA,target,'m','Linewidth',2);
qh = plot(AA,q(AA,mu),'k','Linewidth',2)
legend([th,qh],{'Target, p(A)','Proposal, q(A)'});
xlabel('A');
 
% SOME CONSTANTS
nSamples = 5000;
burnIn = 500;
minn = 0.1; maxx = 5;
 
% INTIIALZE SAMPLER
x = zeros(1 ,nSamples);
x(1) = mu;
t = 1;
 
% RUN METROPOLIS-HASTINGS SAMPLER
while t < nSamples
    t = t+1;
 
    % SAMPLE FROM PROPOSAL
    xStar = exprnd(mu);
 
    % CORRECTION FACTOR
    c = q(x(t-1),mu)/q(xStar,mu);
 
    % CALCULATE THE (CORRECTED) ACCEPTANCE RATIO
    alpha = min([1, p(y,xStar,B)/p(y,x(t-1),B)*c]);
 
    % ACCEPT OR REJECT?
    u = rand;
    if u < alpha
        x(t) = xStar;
    else
        x(t) = x(t-1);
    end
end
 
% DISPLAY MARKOV CHAIN
figure;
subplot(211);
stairs(x(1:t),1:t, 'k');
hold on;
hb = plot([0 maxx/2],[burnIn burnIn],'g--','Linewidth',2)
ylabel('t'); xlabel('samples, A');
set(gca , 'YDir', 'reverse');
ylim([0 t])
axis tight;
xlim([0 maxx]);
title('Markov Chain Path');
legend(hb,'Burnin');
 
% DISPLAY SAMPLES
subplot(212);
nBins = 100;
sampleBins = linspace(minn,maxx,nBins);
counts = hist(x(burnIn:end), sampleBins);
bar(sampleBins, counts/sum(counts), 'k');
xlabel('samples, A' ); ylabel( 'p(A | y)' );
title('Samples');
xlim([0 10])
 
% OVERLAY TARGET DISTRIBUTION
hold on;
plot(AA, target/sum(target) , 'm-', 'LineWidth', 2);
legend('Sampled Distribution',sprintf('Target Posterior'))
axis tight

结语

这里我们探索了如何从 Metropolis 算法一般化得到 Metropolis-Hastings 算法，从而对复杂的（未归一化）的概率分布利用非对称推荐分布进行采样。Metropolis-Hastings 算法的一个缺点是：并非所有的推荐采样都被接受，因此它浪费了宝贵的计算资源。这个缺点在对高维分布进行采样时更明显。这就是吉布斯采样被引入的原因。我们在下一篇文章中将会介绍吉布斯采样，这个采样利用条件概率的优势，可以保留马尔可夫链中所有推荐的状态。