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 算法的一個缺點是：並非所有的推薦采樣都被接受，因此它浪費了寶貴的計算資源。這個缺點在對高維分布進行采樣時更明顯。這就是吉布斯采樣被引入的原因。我們在下一篇文章中將會介紹吉布斯采樣，這個采樣利用條件概率的優勢，可以保留馬爾可夫鏈中所有推薦的狀態。