本文主要譯自 MCMC: The Metropolis Sampler
正如之前的文章討論的，我們可以用一個馬爾可夫鏈來對目標分布 \(p(x)\) 進行采樣，通常情況下對於很多分布 \(p(x)\)，我們無法直接進行采樣。為了實現這樣的目的，我們需要為馬爾可夫鏈設計一個狀態轉移算子（transition operator），是的這個馬爾可夫鏈的穩態分布與目標分布吻合。Metropolis 采樣算法（更通常的是 Metropolis-Hastings 采樣算法）采用簡單的啟發式方法實現了這樣的狀態轉移算子。
Metropolis 采樣
從一個隨機的初始狀態\(x^{(0)}\thicksim \pi^{(0)}\)開始，算法首先從建議分布（proposal distribution）\(q(x|x^{(t-1)})\)生成一個可能的抽樣\(x^*\)。（譯者按：這里的建議分布（proposal distribution）區別於目標分布（target distribution），通常比目標分布要簡單，可以通過計算機直接生成按建議分布的隨機變量）和普通的馬爾可夫鏈轉移算子相同，建議分布僅僅取決於馬爾可夫鏈當前的狀態。然而，Metropolis 算法的狀態轉移算子有一個附加步驟，那就是評估目標分布在建議狀態（proposed state 譯者按：通過建議分布采樣得到的狀態）下是否有足夠大的密度函數來支持這個建議狀態是合理的，從而將其設置為馬爾可夫鏈的下一個狀態。如果密度函數\(p(x)\) 在建議狀態下很小，那么，它就很有可能（但不一定）被拒絕。接受或者拒絕建議狀態的評判標准由下面的啟發式算法定義：
1. 如果\(p(x^*)\geq p(x^{(t-1)})\)，那么建議狀態 \(x^*\) 將被保留為一個采樣點，從而設置為馬爾可夫鏈上的下一個狀態（i.e. 始終向 \(p(x)\) 不變或者變大的方向移動馬爾可夫鏈的狀態）。
2. 如果\(p(x^*)<p(x^{(t-1)})\)，表示 \(p(x)\) 在 \(x^*\) 附近的密度值小於前一個狀態，但是這個建議的狀態 \(x^*\) 仍然有可能被接受，但是是隨機的，接受的概率是\(\frac{p(x^*)}{p(x^{(t-1)})}\)
這些啟發步驟可以被整合到一起，為建議狀態計算接受概率（acceptance probability），
\(\alpha = \mathrm{min}\left(1, \frac{p(x^*)}{p(x^{(t-1)})}\right)\)
掌握了接受概率，metropolis 算法的轉移算子是這樣工作的：如果一個隨機均與分布 \(u\) 小於等於 \(\alpha\)，那么我們就接受 \(x^*\) 作為下一個狀態，反之，這個建議狀態被拒絕， 我們將繼續生成下一個建議狀態。為了用 metropolis 算法得到 \(M\) 個抽樣，我們運行下面的算法：
1. set t = 0 設置 t = 0，表示初始時刻一次迭代都沒發生
2. generate an initial state \(x^{(0)}\) from a prior distribution \(\pi^{(0)}\) over initial states
3. repeat until \(t = M\)
set t=t+1
generate a proposal state \(x^*\) from \(q(x|x^{(t-1)})\)
calculate the acceptance probability \(\alpha = \mathrm{min}\left(1,\frac{p(x^*)}{p(x^{(t-1)})}\right)\)
draw a random number \(u\) from Unif(0,1)
if \(u\leq\alpha\), accept the proposal and set \(x^{(t)} = x^*\)
else set \(x^{(t)}=x^{(t-1)}\)
Example: Using the Metropolis algorithm to sample from an unknown distribution
假如我們有這樣一個神秘的概率分布函數：
\(p(x) = (1+x^2)^{-1}\)
我們想生成符合這個密度函數的隨機采樣。為了使用 Metropolis 采樣算法，我們需要定義兩件事：（1）馬爾可夫鏈初始狀態的先驗分布\(\pi^{(0)}\)，以及（2）建議分布函數 \(q(x|x^{(t-1)})\)（譯者按：這里再強調一下，我們的建議分布函數為 proposal distribution，是比較好采樣的，目標函數 target distribution 是不好采樣的，我們通過好采樣的函數采樣，然后判斷采到的樣本是否被接受）。對於這個例子，我們定義：
\(\pi^{(0)}\thicksim \mathcal{N}(0, 1)\)
\(q(x^{(t)}|x^{(t-1)})\thicksim \mathcal{N}(x^{(t-1)}, 1)\)，
這兩者都是簡單的正態分布，一個以0為中心，另一個以上一個狀態的取值為中心。下面的一塊 MATLAB 代碼采用了上面的推薦分布和先驗分布，運行了 Metropolis 采樣。
% METROPOLIS SAMPLING EXAMPLE
randn('seed',12345);
% DEFINE THE TARGET DISTRIBUTION
p = inline('(1 + x.^2).^-1','x'); % 定義目標分布
% SOME CONSTANTS
nSamples = 5000; % 抽取5000個樣本
burnIn = 500;
nDisplay = 30; % 將展示前30個抽樣過程
sigma = 1; % 建議分布的標准差定義為1
minn = -20; maxx = 20;
xx = 3*minn:.1:3*maxx; % 抽樣的區間
target = p(xx); % 目標函數在區間上的取值，用於下面畫函數圖
pauseDur = .8; % 動態展示過程中每個抽樣需要暫停0.8秒
% INITIALZE SAMPLER
x = zeros(1 ,nSamples); % 初始化抽樣序列
x(1) = randn; % 初始化第一個抽樣
t = 1; % 抽樣編號
% RUN SAMPLER
while t < nSamples % 當沒有抽滿指定數目的樣本時
t = t+1;
% SAMPLE FROM PROPOSAL
xStar = normrnd(x(t-1) ,sigma); % 按照建議抽樣函數進行采樣
proposal = normpdf(xx,x(t-1),sigma); % 建議抽樣函數在整個區間上的取值
% CALCULATE THE ACCEPTANCE PROBABILITY
alpha = min([1, p(xStar)/p(x(t-1))]); % 接受概率
% ACCEPT OR REJECT?
u = rand; % Metropolis 采樣的核心部分，按照接受概率決定是否接受當前采樣
if u < alpha
x(t) = xStar;
str = 'Accepted';
else
x(t) = x(t-1);
str = 'Rejected';
end
% DISPLAY SAMPLING DYNAMICS
if t < nDisplay + 1 % 動態展示前 nDisplay 個采樣過程
figure(1);
subplot(211);
cla
plot(xx,target,'k'); % 目標函數圖像
hold on;
plot(xx,proposal,'r'); % 建議函數圖像
line([x(t-1),x(t-1)],[0 p(x(t-1))],'color','b','linewidth',2) % 目標函數在上一個采樣點的函數值
scatter(xStar,0,'ro','Linewidth',2) % 將當前采樣點表示在x軸上
line([xStar,xStar],[0 p(xStar)],'color','r','Linewidth',2) % 目標分布函數在當前采樣點的函數值
plot(x(1:t),zeros(1,t),'ko') % 將已經采樣過的點表示在x軸上
legend({'Target','Proposal','p(x^{(t-1)})','x^*','p(x^*)','Kept Samples'})
switch str
case 'Rejected'
scatter(xStar,p(xStar),'rx','Linewidth',3) %如果拒絕，就用x表示這個點
case 'Accepted'
scatter(xStar,p(xStar),'rs','Linewidth',3) %如果接受，用方塊表示這個點
end
scatter(x(t-1),p(x(t-1)),'bo','Linewidth',3) % 目標函數在上一個采樣點的函數值
title(sprintf('Sample % d %s',t,str))
xlim([minn,maxx])
subplot(212);
hist(x(1:t),50); colormap hot;
xlim([minn,maxx])
title(['Sample ',str]);
drawnow
pause(pauseDur);
end
end
% DISPLAY MARKOV CHAIN
figure(1); clf
subplot(211);
stairs(x(1:t),1:t, 'k');
hold on;
hb = plot([-10 10],[burnIn burnIn],'b--');
ylabel('t'); xlabel('samples, x');
set(gca , 'YDir', 'reverse');
ylim([0 t]);
axis tight;
xlim([-10 10]);
title('Markov Chain Path');
legend(hb,'Burnin');
% DISPLAY SAMPLES
subplot(212);
nBins = 200;
sampleBins = linspace(minn,maxx,nBins);
counts = hist(x(burnIn:end), sampleBins);
bar(sampleBins, counts/sum(counts), 'k');
xlabel('samples, x' ); ylabel( 'p(x)' );
title('Samples');
% OVERLAY ANALYTIC DENSITY OF STUDENT T
nu = 1;
y = tpdf(sampleBins,nu);
%y = p(sampleBins);
hold on;
plot(sampleBins, y/sum(y) , 'r-', 'LineWidth', 2);
legend('Samples',sprintf('Theoretic\nStudent''s t'))
axis tight
xlim([-10 10]);
采樣過程如下面的GIF圖所示，圖片來自原文
圖中展示了前30次抽樣過程，黑線表示目標分布函數 \(p(x)\)，在 x 軸上移動的紅線表示建議分布函數 \(q(x|x^{(t-1)})\) 。藍色豎線（它在紅色分布函數的軸線上）表示 \(p(x^{(t-1)})\) 的值，紅色豎線表示 \(p(x^*)\)，\(x^*\) 是根據紅色建議分布函數生成的建議采樣。每一次迭代，如果紅色豎線比藍色豎線長，那么采樣 \(x^*\) 將會被接受，建議分布函數（紅色曲線）的中心將會移動到新的采樣點的位置。如果藍色的豎線更長，那么建議采樣將會隨機地被接受或拒絕。
但是為什么要隨機地接受“壞的”（\(p(x^*)<p(x^{(t-1)})\) 的點 \(x^*\)）采樣呢？ 因為這樣做可以使得馬爾可夫鏈時不時地訪問目標分布中概率小的狀態點。這是一個合理的性質，如果我們想要對目標分布的整體（包括尾巴的部分）進行隨機采樣。
Metropolis 算法一個很有吸引力的地方是，它不需要目標分布 \(p(x)\) 歸一化為概率分布。這是因為，接受概率是基於兩個目標分布函數值的比值的，分子分母的歸一化因子 \(Z=\int p(x)\mathrm{dx}\) 可以消掉：
如果 \(p(x)\) 是一個未歸一化的分布函數,那么：
\(p^*(x) = \frac{p(x)}{Z}\)
是其歸一化的分布函數，歸一化常數為 Z，於是:
\(p(x) = Zp^*(x)\)
然后，當我們計算接受概率 \(\alpha\) 時，我們需要計算：
\(\frac{p(a)}{p(b)} = \frac{Zp^*(a)}{Zp^*(b)} = \frac{p^*(a)}{p^*(b)}\)
歸一化常數 Z 被消掉了！這個好性質在貝葉斯方法中十分有用，因為一個分布的歸一化常數不容易直接計算。事實上，我們上面例子中采樣的“神秘”分布是一個沒有歸一化的 1 自由度 Student's-t 分布。我們比較 \(p(x)\) 和 Student's-t 分布的定義：
\(Student(x,\nu) = \frac{\Gamma(\frac{\nu+1}{2})}{\sqrt{\nu\pi}\Gamma(\frac{\nu}{2})}\left(1+\frac{x^2}{\nu}\right)^{-\frac{\nu+1}{2}}\)
自由度為1時，\(\nu=1\)，帶入上式，令常數項為 \(\frac{1}{Z}\)，有：
\(Student(x,\nu) = \frac{(1+x^2)^{-1}}{Z} = \frac{p(x)}{Z}\)
我們看到 \(p(x)\) 是一個自由度為1的 Student's-t 分布，但是少了歸一化常數：
\(Z = \left(\frac{\Gamma(\frac{\nu +1}{2})}{\sqrt{\nu\pi}\Gamma(\frac{\nu}{2})}\right)^{-1}\)
下面是上方代碼的另一個輸出，展示了 Metropolis 采樣得到了符合歸一化了的 Student's-t 分布的隨機數，盡管 \(p(x)\) 沒有歸一化。
上面的輸出展示了馬爾可夫鏈從 \(x^{(0)}\) （上方）到 \(x^{(5000)}\) （下方）更新的過程。這個鏈的 burn in 區間定為500次迭代，它被藍色的虛線表示出來了。
下面的圖像用黑色表示了馬爾可夫鏈采集到的樣例。紅色曲線表示 1 自由度 Student's-t 分布的理論值。我們看到 Metropolis 采樣得到的的樣例符合 Student's-t 分布，盡管轉移算子中用到的函數 \(p(x)\) 並沒有被歸一化成概率分布函數。
轉移算子的可逆性（Reversibility of the transition operator）
事實證明為了使馬爾可夫鏈穩定在一個固定分布（也就是我們所關心的目標分布），狀態轉移算子有一個理論上的約束條件。這個約束條件說的是：\(x^{(t)}\rightarrow x^{(t+1)}\) 的轉移概率與反向 \(x^{(t+1)}\rightarrow x^{(t)}\) 轉移概率相等。這個可逆性通常被稱為細致平衡（detailed balance）。對於 Metropolis 算法的轉移算子，如果推薦概率分布 \(q(x|x^{(t-1)})\) 是對稱的，那么其可逆性將被得到保證。也就是說：\(q(x^{(t+1)}|x^{(t)}) = q(x^{(t)}|x^{(t+1)})\)，對於離散的情況，這個條件的意思是：推薦轉移矩陣是對稱的。連續的情況，有很多概率分布函數也是對稱的，例如正態分布，柯西分布，Student's-t 分布，均勻分布。
上面的例子中，我們的推薦分布函數為方差為1的正態分布，正態分布的均值為上一個函數取值，也就是：
\(q(x^{(t)}|x^{(t-1)})\thicksim\mathcal{N}(x^{(t-1)},1) = \frac{1}{\sqrt{2\pi}}\mathrm{exp}\left(\frac{(x^t-x^{(t-1)})^2}{2}\right)\)
顯然，
\(q(x^{(t)}|x^{(t-1)}) = q(x^{(t-1)}|x^{(t)}) \)
滿足了上面提到的約束條件，也就是：推薦函數是對稱的。
然而，使用對稱的推薦函數或許不能夠准確、高效地對所有目標分布函數取樣。例如，如果一個目標分布定義在 x 軸正半軸，我們希望推薦分布也定義在同樣的區間，這就導致了不對稱性的出現。這就是為什么我們接下來會引入 Metropolis-Hastings 采樣算法。我們接下來的文章將會介紹 Metropolis-Hastings 算法是如何通過簡單地改變接受概率的計算方式，來允許我們使用非對稱的推薦分布的。
譯者按：
這篇文章介紹了 Metropolis 算法，並給出了限制條件。算法可以總結為：
1. 選一個推薦分布（離散有限維空間下為一個轉移矩陣，連續空間下為一系列函數，例如上例中方差為1，均值為當前采樣點的正態分布函數）
2. 根據推薦分布生成采樣 \(x^*\)
3. 計算 \(\alpha = \mathrm{min}(\frac{p(x^*)}{p(x^{(t)})},1)\)
4. 生成0到1之間的均勻分布隨機數\(u\)，如果\(u<\alpha\) 那么就接受 \(x^*\) 作為第 t+1 個采樣，否則拒絕 \(x^*\)，第 t+1 個采樣仍然與 \(x^{(t)}\) 相同
重復上面到足夠多的次數，然后扔掉 burn in 的部分，剩下的就是符合目標分布的隨機采樣點。
使用這個算法，需要注意的一個重要條件是：推薦分布（狀態轉移函數）必須是對稱的，我們的例子中的正態分布正是對稱的。
但是，這篇文章也有不足的地方，那就是沒有告訴我們，為什么采用這種算法就可以得到我們的目標分布，而推薦分布（轉移算子）又為什么必須是對稱的。下面本人將嘗試給出一種不嚴格的證明，或者是解釋，來說明這個算法的確可以得到我們想要的分布。
假如我們希望馬爾可夫鏈可以收斂到目標分布函數 \(\pi(x)\)，轉移算子為 \(T\)，\(T(a,b)\) 表示從狀態 a 轉移到狀態 b 的概率。推薦轉移算子為\(Q\)，\(Q(a,b)\) 表示按照推薦分布從狀態 a 轉移到狀態 b 的概率，根據前面的要求，Q 必須是對稱的：\(Q(a,b)=Q(b,a)\)。那么根據我們的算法，從狀態 a 轉移到狀態 b 的概率等於推薦函數的轉移概率乘以接受率，即：
\(T(a,b) = Q(a,b)*\mathrm{min}\left(1,\frac{\pi(b)}{\pi(a)}\right)\)
那么，如果馬爾可夫鏈達到穩態分布 \(\pi(x)\)（也就是細致平穩條件）：
\(\pi(a)\cdot T(a,b) = \pi(a)\cdot Q(a,b)\mathrm{min}\left(1,\frac{\pi(b)}{\pi(a)}\right)\)
\(= Q(a,b)\mathrm{min}(\pi(a), \pi(b))\)
\(= \pi(b)\cdot Q(a,b)\mathrm{min}\left(\frac{\pi(a)}{\pi(b)}, 1\right)\)
\(= \pi(b)\cdot Q(b,a)\mathrm{min}\left(\frac{\pi(a)}{\pi(b)}, 1\right)\)
\(= \pi(b)\cdot T(b,a)\)
果然，在穩態分布時各個狀態之間的轉移達到平衡。我們同時也知道了，為什么一定要求推薦分布是對稱的。
Metropolis 算法給出了一個可以對任意已知分布進行采樣的方法，這里我們再回顧一下為什么我們要對已知分布進行采樣：
我們要計算:
\(E(x) = \int f(x)p(x)\mathrm{dx}\)
只要產生了足夠多的按照 \(p(x)\) 分布的點：\(x_1,x_2,...x_N\)，就可以將上面的積分近似為：
\(E(x)=\frac{1}{N}\sum\limits_{n=1}^Nf(x_n)\)
應用實例
寫到這里，再舉一個經典的應用 Matropolis 算法的例子："Hard disk in a box model"。
統計物理學家想研究盒子中的氣體的一些物理性質，這些性質是所有氣體分子的位置的函數，我們假設二維的邊長為1的盒子中有 N 個氣體分子，它們的位置由各自的橫、縱坐標唯一確定。那么，每一個位置分布都代表了盒內氣體的一種狀態，狀態空間為：
\(X=(x_1,y_1,x_2,y_2,...,x_N,y_N)\)
玻爾茲曼方程可以計算任何一種分布的概率：
對於任何一個狀態 \(x\in X\)，其概率為：
\(\pi(x) = \frac{1}{Z}\mathrm{exp}(-\frac{\mathcal{E}(x)}{kT})\)
其中 Z 為歸一化常數，k 是玻爾茲曼常數，T 是溫度。\(\mathcal{E}(x)\) 是氣體的能量，將每個氣體分子的位置帶入公式可以算出來。歸一化常數 Z 是所有可能的位置帶入到 \(\pi(x)\) 的分子中，算出來的數再累加的結果，顯然 Z 是不可解的，因為我們無法窮盡氣體分子分布的所有可能情況。
我們如果想計算該系統的某個宏觀參數，事實上，氣體系統的一切宏觀參數都是某種函數的“均值”（或者說期望）。例如壓強 p，可以把它看成某個關於氣體分布狀態 \(x\) 的函數 \(f(x)\) 的期望：
\(p = E(f(x)) = \int f(x)\pi(x)\mathrm{dx}\)
該如何計算呢？顯然根據上面的討論 \(\pi(x)\) 因為歸一化常數 Z 太復雜，無法直接算出來，更別提對它進行積分了。但是我們可以采用蒙特卡洛方法計算這個 \(f(x)\) 的期望:
首先生成符合 \(\pi(x)\) 分布的一系列采樣，\(x_1,x_2,x_3,...x_N\)
然后帶入用公式：
\(E(f(x))\approx\frac{1}{N}\sum_{n=1}^Nf(x_n)\)
求出 \(f(x)\) 的期望。
但是，如何生成關於 \(\pi(x)\) 的采樣呢？我們就要用到本文前面介紹的 Metropolis 算法了：
雖然我們不知道 \(pi(x)\) 的歸一化常數 Z，但是因為它是常數，我們可以在 Metropolis 算法中計算接受概率時把它扔掉，只考慮分子部分。
也就是說，在算法中，我們使用 \(\tilde{\pi}=\mathrm{exp}(-\frac{\mathcal{E}(x)}{kT})\)
下面我們需要確定的是，推薦轉移概率分布（proposal transition distribution），這里我們做如下分析：
1. 只要有一個分子的位置發生改變，那么這個氣體系統的狀態就發生了改變，因此，每次我們從 N 個分子中挑出一個分子，改變它的位置，從而獲得一個新的狀態。
2. 一個氣體分子雖然可以隨機移動，但是基本上只會在周圍一定的空間內隨機移動，而不會發生從盒子左上角一下子移動到右下角的情況，因此，我們認為每次單個分子的移動范圍在以它為中心，\(2\alpha\) 為邊長的正方形盒子里。
3. 單個分子在以它為中心 \(2\alpha\) 為半徑的盒子內是均勻隨機移動的。
由於每次隨機選取的點是均勻分布的，因此這里的推薦轉移分布是對稱的，那么，我們可以按這樣的步驟采樣：
1. 隨機從 N 個分子中抽取一個分子
2. 將這個分子按均勻分布隨機移動到以它為中心，\(2\alpha\) 為邊長的盒子中的某個點，得到新的狀態 \(x^*\)
3. 計算接受概率\(p(accept) = \mathrm{min}(1,\frac{\tilde{\pi}(x^*)}{\tilde{\pi}(x^{(t)})})\)
4. 隨機接受或者拒絕這個取樣
5. 重復上面1至4步，得到足夠多的取樣
通過上面的步驟，我們就可以采集到符合 \(\pi(x)\) 分布的樣本了，接下來將它們代入蒙特卡洛方法的公式中，計算出系統的壓強 p 。當然，這里的壓強也可以換成別的什么宏觀參數。