Metropolis-Hastings algorithm
1. 隨機模擬的基本思想
假設我們有一個矩形區域\(R\)，面積為\(S_0\)。在此區域中，有一個不規則區域\(M\)，其面積\(S\)待求。
方法1：把不規則區域\(M\)划分為多個小的規則區域，由這些規則區域的面積總和\(S'\)近似。
方法2：抓一把黃豆，均勻鋪在\(R\)中，再統計落在不規則區域\(M\)中的黃豆比例。該比例可近似\(\frac{S}{S_0}\)。
方法2采用的是抽樣（采樣）解決問題的思想，即隨機模擬。
因此，隨機模擬的關鍵，在於如何將實際復雜問題，轉化為抽樣可以解決的問題。當然，這非常靈活，也考驗我們的創造力。
我們接下來的核心問題，是利用計算機可直接實現的均勻抽樣，借助隨機模擬的方法，實現滿足特定概率分布\(p(x)\)的抽樣。
2. 拒絕抽樣
當我們的概率分布很簡單時，如\([0,1]\)上的均勻分布，抽樣是很好實現的。現成的隨機算法非常多。
但是，對於其他更復雜的分布，如果我們還想借助簡單的均勻隨機抽樣，就需要換換思路了。
例如上圖中的\(p(z)\)（忽略波浪號），是較為復雜的概率密度函數。
我們先構造比較簡單的、容易抽樣的高斯分布\(q(z)\)，乘一個常數\(k\)，使之滿足：
拒絕抽樣的方法是：
-
對高斯分布的樣本進行采樣，得到一個\(z_0\)，如圖；
-
在\([0,kq(z_0)]\)上均勻采樣，得到\(u_0\)；
-
如果\(u_0 > p(z)\)，就拒絕該采樣結果；反之接受；
-
重復直到接受足夠多的樣本。
我們理解一下為什么拒絕抽樣是可行的，這對我們之后運用隨機模擬思想非常關鍵！
-
該簡單分布與高斯分布的趨勢大致相同，中間大兩頭小，因此取到中間的\(z\)的概率比較大；
-
通過\(u_0\)進一步篩選\(z\)，\(p(z)\)較低時拒絕率較大，反之接受率較高，大致符合\(p(z)\)分布要求。
一句話，通過簡單抽樣和設置接受率，“強迫”結果趨於復雜分布。
這樣，我們就通過簡單的高斯分布采樣，以及計算機可直接實現的均勻采樣，間接實現了復雜的\(p(z)\)采樣。
然而，在高維情況下，該方法的劣勢很明顯：
-
合適的簡單分布\(q(z)\)很難找到。
-
合適的\(k\)值很難找到。
這兩個問題會導致拒絕率極高，計算效率很低。其根本缺陷在於：樣本之間是獨立無關的。
3. Metropolis-Hastings抽樣
3.1 引入思想
結合馬氏鏈知識（參考信息論相關書籍），我們知道：
對於遍歷的馬爾可夫鏈，當其轉移次數足夠多時，將會進入穩態分布\(\pi (x)\)，即各狀態的出現概率趨於穩定。
進一步，假設變量\(\boldsymbol X\)所有可能的取值，構成某遍歷馬氏鏈的狀態空間。
則無論從什么狀態出發，只要轉移步數足夠大，該馬氏鏈將趨於穩定，即逐漸開始依次輸出符合\(p(x)\)分布的狀態\(X_i\)。
這樣，通過收集馬氏鏈收斂后產生的\(X_i\)，我們就得到了符合分布\(p(x)\)的樣本。
此時穩態分布即\(\pi (x)=p(x)\)。
比如，我們希望抽樣樣本滿足指數分布。此時，整數0到正無窮都是狀態，但整數0的出現概率最大，隨着數增大，出現概率越來越小。我們構造一個穩定輸出服從指數分布狀態的馬氏鏈，則得到的樣本服從指數分布。
這個絕妙的想法在1953年由Metropolis提出。為了研究粒子系統的平穩性質，Metropolis 考慮了物理學中常見的波爾茲曼分布的采樣問題，首次提出了基於馬氏鏈的蒙特卡羅方法（Metropolis算法），並在最早的計算機上編程實現。
Metropolis 算法是首個普適的采樣方法，並啟發了一系列 MCMC方法，所以人們把它視為隨機模擬技術騰飛的起點。 Metropolis的這篇論文被收錄在《統計學中的重大突破》中， Metropolis算法也被選為二十世紀的十個最重要的算法之一。
我們接下來介紹的MH算法是Metropolis 算法的一個常用的改進變種。
由於馬氏鏈的收斂性質主要由轉移矩陣\(\boldsymbol P\)決定, 所以基於馬氏鏈做采樣的關鍵問題是：如何構造轉移矩陣\(\boldsymbol P\)，使平穩分布恰好是我們要的分布\(p(x)\)。
3.2 理論基礎：細致平穩條件
定理——細致平穩條件：
若非周期馬氏鏈的轉移矩陣\(\boldsymbol P\)和分布\(\pi (x)\)滿足：
則分布\(\pi (x)\)是馬氏鏈的平穩分布。
證明：
這說明：
因此\(\pi (x)\)是平穩分布。
物理意義：
對於任意兩個狀態\(i,j\)，從狀態\(i\)流失到狀態\(j\)的概率質量，等於反向流動的概率質量，因此是狀態概率是穩定的。
3.3 M-H算法實現
假設我們已經有了一個轉移矩陣\(\boldsymbol Q\)。一般情況下，該轉移矩陣不滿足細致平穩條件：
我們引入接受率\(\alpha(i,j)\)，滿足：
其中：
-
\(P(i)\)是狀態\(i\)出現的概率，是假設馬氏鏈滿足復雜分布時輸出狀態的概率。
-
\(\boldsymbol Q\)是初始轉移概率矩陣，滿足任意一個給定的簡單分布，便於抽樣即可。
顯然，接受率最簡單的構造方法是：
此時，細致平穩條件成立，該馬氏鏈輸出的狀態滿足穩態分布\(p(x)\)！
綜上，有幾個要點必須要理解，有助於我們編程實現：
-
我們引入接受率，使得該馬氏鏈以轉移概率\(\boldsymbol Q_{ij}\alpha(i,j)\)從狀態\(i\)轉移到狀態\(j\)。
-
該轉移概率的實現方式是：我們先按\(\boldsymbol Q_{ij}\)生成新狀態，再按\(\alpha(i,j)\)的概率接受轉移結果，乘積即為\(\boldsymbol Q_{ij}\alpha(i,j)\)轉移概率。
-
一定要理解這個接受率！在拒絕抽樣中，如果\(u_0 > p(z)\)，我們會放棄抽樣結果，不納入最終的統計。正是因為如此，抽樣才會逼近復雜分布。而這里的接受是“接受轉移”，如果\(u>\alpha\)，意味着\(t+1\)時刻狀態和\(t\)時刻相同，並且應該納入統計而不是放棄結果！！
-
由於狀態出現概率和轉移概率滿足細致平穩條件，因此狀態出現概率即穩態分布概率。長期轉移后，我們會得到想要的抽樣分布效果。
圖解轉移過程：
要注意的是，當馬氏鏈處在狀態\(i\)時，我們並不知道如何選擇下一個狀態\(j\)。
我們在這里采用抽樣的方法，並借助接受率，“強迫”轉移過程逼近理想的遍歷馬氏鏈轉移過程。
算法描述如下：
-
初始化馬氏鏈狀態\(X_0=x_0\)
-
從\(t=0\)開始，重復以下步驟：
-
假設該\(t\)時刻狀態為\(X_t=x_t\)；
-
對\(Q_{x_t,x}\)抽樣，隨機得到一個狀態\(x_{next}\)；
-
在\([0,1]\)上均勻抽樣，得到一個\(u\)；
-
若\(u<\alpha(x_t,x_{next})=P(x_{next})\boldsymbol Q_{x_{next},x_t}\)，則接受轉移\(X_{t+1}=x_{next}\)；
-
否則不轉移，即\(X_{t+1}=x_t\)。
-
3.4 算法升級
如果\(\alpha(x_t,x_{next})\)太小，會導致轉移很難發生，馬氏鏈收斂過慢。
我們回到細致平穩條件：
我們希望在保持條件成立的前提下，使接受率盡量增大。
由於接受率的本質是概率，因此\(\alpha(i,j)\)和\(\alpha(j,i)\)中的較大者不應大於1。那么我們作下述改進即可：
解釋：
-
如果\(\alpha(i,j)\)更大，那么應設為1。而第一項大於1，因此\(\alpha(i,j)=1\)，結束。
-
如果\(\alpha(j,i)\)更大，那么為了等式成立，\(\alpha(i,j)\)必須等於buf。而buf\(<1\)，得逞～
綜上，我們得到了升級版算法：
3.5 仿真實驗
仿真目標：
用標准正態分布模擬指數分布（\(\lambda=0.1, \mu=\frac{1}{\lambda}=10\)）。
簡化：
由於高斯分布是一個徑向基函數（取值僅僅依賴於離原點距離的實值函數），因此該矩陣是轉置對稱的，\(\boldsymbol Q_{ij}\)=\(\boldsymbol Q_{ji}\)
代碼實現：
function SamplefromExponential
clc;close all;clear;
% 隨意設
data_size=100000; % 越大越准
start=1; % 非負即可
% 初始化
data=start;
index=1;
% 直接生成多個服從均值為10的指數分布的樣本
data_dir=exprnd(10,[1,data_size]);
figure;
subplot(1,2,1)
histogram(data_dir,'BinWidth',1);
axis([0,100,0,data_size/10]);
% 隨機模擬
while(index<data_size)
u=rand; % 均勻抽樣
data_next=normrnd(data(index),1); % 簡單分布為標准正態分布
buf=exppdf(data_next,10)/exppdf(data(index),10);
alpha=min(1,buf);
if u<alpha
data=[data data_next];
else
data=[data data(index)];
end
index = index+1;
end
subplot(1,2,2)
histogram(data,'BinWidth',1);
axis([0,100,0,data_size/10]);
結果：
