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)} ightarrow x^{(t)}) 的转移概率必须等于 (x^{(t)} ightarrow x^{(t-1)}) 的转移概率，这个条件被称为可逆性或者细致平稳。然而，一个对称的分布或许对很多问题并不适合，比如我们想对定义在 ([0,infty)) 的分布进行采样。

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

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

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

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

1. set t = 0

2. generate an initial state (x^{(0)} hicksimpi^{(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)})} imes c ight))

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

　　　　if (uleqalpha) 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 ightarrow 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)}cdotfrac{pi(b)}{pi(a)} ight))

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

(pi(a)cdot T(a,b) = pi(a)Q(a,b)cdotalpha)

(=pi(a)Q(a,b)cdotmathrm{min}left(1,frac{Q(b,a)}{Q(a,b)}cdotfrac{pi(b)}{pi(a)} ight))

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

(=pi(b)Q(b,a)cdotmathrm{min}left(frac{pi(a)Q(a,b)}{pi(b)Q(b,a)},1 ight))

(=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) eq pi(b)cdot Q(b,a))

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

(pi(a)cdot Q(a,b)cdotalpha(a,b) = pi(b)cdot Q(b,a)cdotalpha(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，那么这就大大提高了接受率而不改变细致平衡等式。因此，我们的接受率 ( ilde{alpha}(a,b)) 可以表示成：

( ilde{alpha}(a,b) = mathrm{min}left(frac{pi(b)Q(b,a)}{pi(a)Q(a,b)}, 1 ight))

上面式子的含义是：

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 ight)=1)，即，接受率为1，我们接受(a ightarrow 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 ight)=frac{pi(b)Q(b,a)}{pi(a)Q(a,b)})，接受率为 (frac{pi(b)Q(b,a)}{pi(a)Q(a,b)})，以这个接受率接受 (a ightarrow b) 的转移。

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

Example: Sampling from a Bayesian posterior with improper prior

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

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

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

(p( heta|y)=frac{p(y| heta)p( heta)}{p(y)})

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

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

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

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

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

( heta = [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 算法的一个缺点是：并非所有的推荐采样都被接受，因此它浪费了宝贵的计算资源。这个缺点在对高维分布进行采样时更明显。这就是吉布斯采样被引入的原因。我们在下一篇文章中将会介绍吉布斯采样，这个采样利用条件概率的优势，可以保留马尔可夫链中所有推荐的状态。