[數據分析] Markov Chain Monte Carlo

01-26

Markov Chain Monte Carlo簡稱MCMC，是一個抽樣方法，用於解決難以直接抽樣的分布的隨機抽樣模擬問題。

在基礎概率課我們有學過，已知一個概率分布函數F(X)，那麼用電腦產生服從Uniform分布的隨機數U，代入 $F^{-1}(X)$ ，那麼 $X=F^{-1}(U)$ 就是服從F(X)的隨機變數。這個方法在金融領域使用很廣，即Monte Carlo Simulation：通過模擬股票價格來計算期權價格。當然，這個方法很重要的一個假設是，已知概率分布函數。但很多情況下我們並沒有方法得到概率分布和其反函數，因此需要其他手段。

本文翻譯學習自 CHRISTOPHE ANDRIEU等人所寫的An Introduction to MCMC for Machine Learning。

文章主要介紹了MCMC的重要理念，實現MCMC的流行演算法，以及前沿研究。

動因

MCMC通常用於解決高維度的積分和最優化問題，這兩種問題也是機器學習、物理、統計、計量等領域的基礎。比如：

貝葉斯推論和學習，給定未知變數x和y，貝葉斯統計主要需計算以下積分：

標準化（normalization），即貝葉斯估計的分母 $p(x|y)=frac{p(y|x)p(x)}{int_{cal X}p(y|x)p(x)dx}$
邊緣分布（marginalisation），給定聯合後驗分布f（x，z），計算 $p(x|y)=int_{cal Z}p(x,z|y)dz$
期望（Expectation）， $E_{p(x|y)}(f(x))=int_{cal X}f(x)p(x|y)dx$ ，方差也是期望的一種

統計力學，用概率理論來總結一個力學系統的平均行為。。。在計算中有運用類似統計推論里計算標準化常數的方法。
最優化，最小化或最大化目標函數
帶有懲罰函數的似然模型選擇（Penalized likelihood model selection），分為兩個步驟，一，對每一個模型計算極大似然估計參數。然後用懲罰項（如AIC，BIC）來選擇複雜度和擬合度都較好的模型。

Monte Carlo原理

依據給定的目標概率分布（高維度的），產生獨立同分布的樣本，利用產生的樣本可以對該分布進行離散式近似 $p_N(x)=frac{1}{N}sum^N_{i=1}delta_{x^{(i)}}(x)$ ，其中 $delta_{x^{(i)}}(x)$ 是xi的狄拉克δ函數。然後就可以用這個近似的分布函數計算各種積分 $I(f(x))=frac{1}{N}sum_{i=1}^Nf(x_i)$ 。但是，很多分布是難以直接取樣的，需要藉助一些方法，就是用可以直接抽樣的較為簡單的分布作為依據（proposal distribution）。如rejection sampling, importance sampling 和 MCMC.

狄拉克δ函數delta-Dirac，它在除x以外的點上都等於零，且其在整個定義域上的積分等於1。是一種連續函數的表達式，若用離散函數來表達的話，就是 $p_N(x)=frac{1}{N}sum^N_{i=1}mathbb I{x=x_i}$

接受-拒絕抽樣 (Rejection sampling)

也叫Accept-Reject Method，如一個未知分布p(x)，不過能夠知道他的上限。那麼選取一個很好抽樣的分布函數q(x)滿足 $p(x)le Mq(x),M<infty$ 。然後，執行如下演算法：

for i = 1 to N: 產生一個服從均勻分布的隨機量u，和服從q(x)的隨機量x if u<p(x)/M*q(x)，那麼x就可以選位一個樣本（accept） i+=1 else：reject

可以證明，如此這般，accepted x服從p(x)。

因為 $Pr(X|Accept) = frac{Pr(Accept|X)Pr(X)}{Pr(Accept)}$
$Pr(Accept|X)=Pr(u<frac{p(x)}{Mq(x)})=frac{p(x)}{Mq(x)}$
$Pr(X)=q(X)$
$Pr(Accept)=int_X Pr(Accept,X)dX=int_X frac{p(x)}{Mq(x)}q(x)dx=frac{1}{M}$
所以 $Pr(X|Accept) =p(X)$

但是這個演算法的限制在於：不一定經常能有 $p(x)le Mq(x),M<infty$ 這樣的關係存在，以及如果M太大的話，抽樣的速度會很慢，因為 $p(accept)propto 1/M$

用python對p(x)為標準正態分布作了個實驗（code在文末）

重要性抽樣 (Importance Sampling)

目標分布p(x)，參考分布q(x)，且p(x)=w(x)q(x)，w(x)稱為重要性權重。比如要計算 $I(f)=int f(x)w(x)q(x)dx$ ，可以通過產生N個獨立且服從q(x)的隨機量，再根據w(x)，就可以得到I（f）的估計 $I_N (f)=sum_{i=1}^Nf(x)w(x)$ 。（這裡不需要乘以q(x)!，因為是sample的匯總，已經把概率分布融在裡面了）。離散式近似 $p_N(x)=sum^N_{i=1}w(x)delta_{x^{(i)}}(x)$ 。選取q(x)令f(x)w(x)的方差最小，即最小化 $E_{q(x)}(f^2(x)w^2(x))$ ，可以得到最優參考分布 $q(x)propto |f(x)|p(x)$ ，說明選取q(x)的時候若多考慮 $|f(x)|p(x)$ 值很大的部分（重要區域important region）就可以得到很高的效率。這意味著，根據重要區域選取q(x)來抽樣得到的結果甚至比直接使用p(x)還好！

咋看之下很奇怪，但是其實是用於應對極端情況問題的抽樣，因此很適合用在金融領域。比如我們已知某風險因子（risk factor）服從p(X)，然後根據風險因子計算的風險值f(X)只在X很大的時候才有值，其他時候均為0，即 $f(X)>0 for Xin A, and =0 otherwise$ ，且 $Pr(Xin A)$ 很小。那麼，如果單純使用p(X)來抽樣可能每抽樣上千次才有一個使得f(x)>0。但是用q(x)=c|f(x)|p(x)，就可以只關注重要的部分。

就是關注的問題不一樣，一般是希望找到p(x)的估計，但這裡我們關注的是基於p(x)的積分的估計。

下面是兩個拓展：

隨著問題的維數增大，參考分布也很難找，所以需要採用參數化參考分布 $q(x, heta)$ ，並在實際模擬中調整參數。如此就引出一種稱為Adaptive importance sampling的方法。一種很潮的方法是對前述方差 $E_{q(x)}(f^2(x)w^2(x))$ 對 $heta$ 求導 $D( heta)=E_{q(x, heta)}(f^2(x)w(x, heta)frac{partial w(x, heta)}{partial heta}$ 。然後用Newton method更新 $heta$ ， $heta_{t+1}= heta_t-alphafrac{1}{N}sum^N_{i=1}f^2(x)w(x, heta_t)frac{partial w(x, heta)}{partial heta}$ 。
若p(x)的標準化常數無法計算，還是可以計算I(f)。 $I(f)=frac{int f(x)w(x)q(x)dx}{int w(x)q(x)dx}$ ，然後w(x)這時候就是真實權重的常數倍。I(f)的Monte Carlo估計 $I(f)=frac{1/Nsum f(x)w(x)}{1/Nsum w(x)}=sum_{i=1}^Nf(x_i)hat w(x_i)$ ， $hat w(x)$ 就是標準化後的重要性權重。

標準化常數（Normalizing constant）：就是使概率函數變為概率密度函數的常數，比如 $int p(x)dx=c$ ，那麼令 $p^*(x)=frac{1}{c}p(x)$ ， $int p^*(x)dx=1$ 。

如果，不僅是想算積分，還想得到服從p(x)的樣本，那就需要再抽樣。因為前面我們只需要權重w(x)來計算積分，而得到的樣本還是服從參考分布的。為了得到真正服從目標分布的樣本，需要藉助一個技術稱為Sampling Importance Resampling(SIR)。比如想要獲得M個樣本，就需要重複M次對 $p_N(x)$ 抽樣，產生M個服從 $p_N(x)$ 的樣本，然後預測的目標分布就是 $ilde p_M(x)=frac{1}{M}sum^M_{i=1}delta_{x_i}(x)$ 。

根據前述 $p_N(x)=sum^N_{i=1}w(x)delta_{x^{(i)}}(x)$ ，其實就是N個sample xi每一個被選中的概率約為w(xi)，然後依此來抽樣。用p(x)等於標準正態，q(x)等於均勻分布，在python中實驗結果如下（code在文末）。

不過即使有上述幾個拓展，有的時候還是很難得到目標分布的良好估計。因此，我們需要更加英霸的技術，基於馬爾科夫鏈。

馬爾科夫鏈蒙地卡羅（MCMC）

是我覺得目前來看縮寫最好記的方法，值得一學。

簡而言之，就是用馬爾科夫鏈方法來抽樣。本文一貫的假設是目標分布p(x)難以直接抽樣，但是我們對他的常數倍了如指掌（原文是「can evaluate p(x) up to a normalising constant」）

馬爾科夫鏈

首先，xi是離散的隨機序列（stochastic process），而且在每個時刻只有s個可能的值（state）。那麼，這個隨機序列若滿足以下條件就可以稱為馬爾科夫鏈： $p(x_i|x_{i-1},...x_1)=T(x_i|x_{i-1})$ 。即每一時刻的條件概率分布只由上一期和一個轉換函數/矩陣決定。

在隨機微積分里對馬爾科夫序列的定義是：對於任意時刻i，和任意函數f，都存在函數g，使得 $E_n[f(X_{n+1})]=g(X_n)$ ，也即未來的序列的期望值完全由現在的實現值的某種函數決定，與過去的實現值無關。

令T表示轉換概率矩陣，即第i行第j列表示x的值在當期等於i的條件下在下一期轉換到j的概率。若在1時刻，x的概率分布為 $mu(x_1)$ ，那麼在2時刻，x的概率分布就是 $mu(x_1)T$ 。馬爾科夫鏈神奇的地方就是，只要T滿足一定的條件，那麼 $mu(x_1)T^N$ 會收斂，而且不論一開始的x的概率分布是咋樣，最後收斂的值都一樣，為x的非條件概率分布（unconditional distribution）。

比如說下圖的關係就可以用 $x_i=[x_1,x_2,x_3]$ ， $T=egin{bmatrix} 0 & 1 & 0 \ 0 & 0.1 & 0.9 \ 0.6 & 0.4 & 0 end{bmatrix}$ 來表示

T需要滿足的條件是：

不可化簡，即對於x任何可能的取值，都有機會（概率大於0）到達其他的取值（不一定要下一期，只要是在有限的時刻內），即上面的圖是閉環的。
不能有死循環

另外，還有一個p(x)是可收斂分布的充分條件（非必要）是 $p(x_i)T(x_{i-1}|x_i)=p(x_{i-1})T(x_i|x_{i-1})$

用互聯網來具象化上述兩個條件的重要意義：第一個條件說明我們從任意一個網頁出發都可以訪問其他任何網頁，而且在任何一個網頁都必須可以到其他網頁，而不是死循環。例如谷歌的PageRank演算法就是定義T=L+E，其中L矩陣的各元素表示從一個網頁到另一個網頁的鏈接數，E是一個均勻隨機數矩陣，用於滿足不可化簡性和排除死循環（再次的網站也要有正的概率可以擺在搜索排行前面，雖然概率很低）。

MCMC的目的就是用Markov Chain收斂於p(x)這一性質來抽樣，即不管一開始抽取的樣本服從什麼分布，只要用了恰當的轉換矩陣，經過多次轉換後之後抽取的眼本就會服從目標分布。概念說起來很簡單，但是轉換矩陣的不同選擇就形成了不同的演算法。

在X是連續變數的情況下，轉換矩陣T變為核函數積分形式： $p(x_{i+1})=int p(x_i)K(x_{i+1}|x_i)dx_i$ ， $K(x_{i+1}|x_i)=Pr(x_{i+1}|x_i)$

Metropolis Hasting Algorithm

最出名的演算法叫做Metropolis Hasting Algorithm （MH）：很多其他實用的MCMC演算法究其根本都是MH演算法的延伸。

演算法介紹：對於目標分布p(x)，首先給定參考條件概率分布q(x*|x)，然後基於當前x的值模擬出一個新樣本x*，然後馬爾科夫鏈將依據一定的概率移動到x*（否則還留在x），這個概率是 $A(x,x^*)=minleft{ 1,frac{p(x^*)q(x|x^*)}{p(x)q(x^*|x)} ight}$ 。具體演算法偽代碼如下：

1. 隨便整一個x0，比如都為02. for i = 0 to N-1: 產生一個服從均勻分布的變數u 產生一個服從q(x*|x)的變數x* 如果u<A(x_i,x*): x_i+1 = x* 否則： x_i+1 = x_i

對偽代碼的解釋就是，如果 $frac{p(x^*)q(x|x^*)}{p(x)q(x^*|x)}>1$ （粗看就是需要p(x*)>p(x)），那麼鏈一定會移動到x*，但如果<1，那麼有一定的概率會移動到x*，意思就是，如果下一個建議的值計算得出的概率密度函數值更大，那麼這個值被保留在樣本里的可能性更高。基本就是直方圖的構建過程。

個人在stackexchange上看到一個很形象的比喻：假設你目前在一片高低不平的山地上，你此行的目的是在海拔越高的地方停留越久（p(x)大的時候，樣本里x就多）。你的方法是隨便指一個新的地方，如果這個地方的海拔更高，那麼就移動過去；但如果這個地方的海拔比當前低，你就拋一個不均勻的硬幣決定是否過去，而硬幣的不均勻程度相當於新海拔和當前海拔的比例。也即新海拔若是當前海拔的一半，你就只有1/2的概率會過去。MH演算法中q(x*|x)就是按照一定的規律指出一個新方向，A(x,x*)就是計算相對高度。

用python實現結果如下：

MH演算法很簡單，但是對q(x*|x)的選擇非常關鍵。而且其他不同的演算法也主要是對q的不同選擇。

MH演算法的轉換核(transition kernel)是 $K_{MH}(x_{i+1}|x_i)=q(x_{i+1}|x_i)A(x_i|x_{i+1})+delta_{x_i}(x_{i+1})r(x_i)$ ，其中 $r(x_i)=int_{cal X}q(x^*|x_i)(1-A(x_i,x^*))dx^*$ 。這個轉換核滿足平衡條件 $p(x_i)K_{MH}(x_{i+1}|x_i)=p(x_{i+1})K_{MH}(x_i|x_{i+1})$ 。

轉換核就是前面的T在連續情況下

可以證明，通過MH演算法得到的樣本就是近似服從目標分布p(x)的樣本：

證明：
首先，一個轉換核只要滿足不可簡化性和無死循環，那麼這個鏈就一定會收斂於一個概率分布。而這樣的轉換核對應唯一的收斂概率，即 $p(x_{i+1})=int p(x_i)K(x_{i+1}|x_i)dx_i$ 。(即滿足這個式子的分布就是收斂分布)
如果轉換核滿足平衡條件 $p(x_i)K_{MH}(x_{i+1}|x_i)=p(x_{i+1})K_{MH}(x_i|x_{i+1})$ ，p(x)是這個馬爾科夫鏈的收斂概率分布，即長遠來看，xi將會逐漸服從p(x)。因為 $int p(x_{i})K(x_{i+1}|x_i)dx_i=int p(x_{i+1})K(x_i|x_{i+1})dx_i=p(x_{i+1})int K(x_i|x_{i+1})dx_i=p(x_{i+1})$ 。因此，MH的轉換核滿足平衡條件，所以抽樣得到的樣本將會趨向於p(x)。

因為演算法內部有拒絕的設定，所以不存在死循環（即每一次產生一個新的proposal x*，都有機會移動到該值上）。為了證明MH演算法不存在可簡化性，需要證明在有限步內x的每一個取值都有正的概率。

MH演算法有兩個簡易的特例：Independent Sampler和Metropolis Algorithm。

Independent Sampler使用的q(x*|x)=q(x*)，即獨立於當前實現值。所以 $A(x_i,x^*)=minleft{1,frac{p(x^*)q(x_i)}{p(x_i)q(x^*)} ight}=minleft{1,frac{w(x^*)}{w(x_i)} ight}$ ，w(x)就是前面重要性抽樣里的重要性權重。這個方法類似重要性抽樣，但是這個方法的樣本是相關的。
Metropolis Algorithm：使用對稱的隨機遊走（random walk）參考分布q，q(x*|x)=q(x|x*)，所以 $A(x_i,x^*)=minleft{1,frac{p(x^*)}{p(x_i)} ight}$

另外，MH演算法有幾個重要性質：

目標分布p(x)的標準化常數不需要知道；
可以改進前面的偽代碼實現同時製作多條鏈；
對q的選擇決定了MH是否可以成功。

------------------接下去就是MCMC如何解決具體問題-----------------

尋找p(x)最大值（simulated annealing）

即通過MCMC尋找p(x)的眾數。單純地使用 $hat x = arg max p(x_i)$ 效率很低，因為在遠離眾數的區間浪費太多計算資源。這裡介紹一種模擬退火(simulated annealing)方法，即對前述MCMC作出調整，使得其不在收斂於p(x)而是 $p_i(x)propto p^{1/T_i}(x)$ ，其中Ti是一個冷卻因子，隨著i不斷增大，其趨向於0。依據是 $p^{infty}(x)$ 的圖像將會集中在p(x)最大值處。而方法的成敗受到q和T的影響。偽代碼如下： $A(x_i,x^*)=minleft{1,frac{p^{1/T_i}(x^*)q(x|x^*)}{p^{1/T_i}(x_i)q(x^*|x_i)} ight}$

1. 隨便整一個x0，比如都為0，另T0=12. for i = 0 to N-1: 產生一個服從均勻分布的變數u 產生一個服從q(x*|x)的變數x* 如果u<A(x_i,x*): x_i+1 = x* 否則： x_i+1 = x_i 更新Ti

效果如圖：

轉換核的混搭和循環

MCMC還有一個優點就是可以結合多種抽樣器於一個混合的（mixture）或循環的（cycle）抽樣器。即若K1和K2是具有相同收斂分布p的核，那麼循環核K1K2和混搭核vK1+(1-v)K2, $0le vle1$ ，都是具有收斂分布p的核。

混搭核可以通過一個總體參考分布先探索可能的取值範圍，然後再針對不同的抽樣器選用不同的參考函數來探索目標分布的性質。這對於有很多極值的目標分布很有效。

循環核可以很好地分離多元變數形成不同組分並分別在演算法中更新。循環核應用中有一個很流行的叫做Gibbs抽樣，通過使用全條件分布 $p(x_j|x_{-j})=p(x_j|x_1,...,x_{j-1},x_{j+1},...,x_n)$ （xi表示第i維的input，即input矩陣的第i列）作為參考概率分布。

Gibbs 抽樣

數據有n個維度，對於各維度，建議使用參考分布 $q(x^*|x_i)=egin{cases}p(x^*_j|x_{-j,i}) ext{ if }x^*_{-j}=x_{-j,i} \0 ext{ Otherwise}end{cases}$ ，這般設置後，A=1。其最有名的就是應用在有向非循環圖（directed acyclic graphs)

Monte Carlo EM

EM（期望-最大化）演算法是極大似然點估計的一個標準演算法。若X包含顯性和隱性變數（不可直接觀察值）那麼似然率 $p(x_v| heta)$ 的一個局部最大值可以通過如下兩個步驟循環得到：

E步：基於隱變數計算似然函數的均值 $Q( heta)=int_{{cal X}h}log(p(x_h,x_v| heta))p(x_h|x_v, heta^{old})dx_h$
M步：選取 $heta$ 最大化 $Q( heta)$

其中，E步就可以用MCMC來抽取 $p(x_h|x_v, heta^{old})$ 的樣本。

輔助變數抽樣（Auxiliary variable samplers)

有時聯合分布比邊緣分布好抽樣，所以可以先通過聯合分布得到樣本，然後忽略掉作為輔助的變數。下面介紹兩種。

混血蒙地卡羅（Hybrid Monte Carlo 無腦直譯）：通過包含目標分布的梯度信息來提高高維度抽樣的效率。下面介紹「青蛙跳」HMC演算法：首先假設目標分布p(x)可導。聯合分布p(x,u)定義為 $p(x){cal N}(u;0,I_{nx})$ ，再定義梯度向量 $Delta (x)=partial log p(x)/partial x$ ，以及一個固定的步長 $ho$ 。在抽樣過程中，每一步中我們都重複L次（自行決定）「青蛙跳」，就是u和x用對方的值進行更新。這般重複後得到的x和u就是提議的樣本，進入下一個環節，即「接受-拒絕」考驗。L和 $ho$ 的選擇決定了演算法的速度。大的L可以得到與初始值差距很大的提議值，但是算的滿。另外，對不同維度的x設置不同的 $ho$ 也會提高效率。
切片抽樣（slice sampler）：是Gibbs抽樣的一個概括。輔助變數和目標變數組成的聯合分布是 $p^*(x,u)=egin{cases}1 if 0le ule p(x) \0 otherwiseend{cases}$ ，所以 $int p^*(x,u)du=int ^{p(x)}_0 du=p(x)$ 。條件概率分布為 $p(u|x)=U_{[0,p(x)]}(u), p(x|u)=U_A(x), where A={x:p(x)ge u}$ 。演算法就是基於當前步得到的sample xi，先在0到p(xi)之間均勻抽樣一個u，然後在p(x)>u的區域A內均勻抽樣一個x。

可逆跳躍MCMC（Reversible jump MCMC）（最新補充）

涉及模型選擇問題，如選擇神經網路的神經元個數、多元自適應Spline回歸的splines個數、或者時間序列自回歸模型的延後期（lag）個數、以及最優的輸入數據集。

給定一類模型集合M，我們將構建遍歷馬爾科夫鏈(ergodic markov china)，使得其收斂分布為p(m,xm)。可以實現模擬以及對比不同維度的分布。

遍歷馬爾科夫鏈：若馬爾科夫鏈上的每一個狀態（state）都是可遍歷的，那麼這個MC就是一個EMC，其實就是具有不可化簡性，可以從任何狀態到達任何狀態。

根據Yanan Fan and Scott A. Sisson的Paper Reversible Jump Markov chain MonteCarlo，MH演算法是Reversible jump MCMC是的一個特殊形式，後者可以應用於更多狀態空間。

給定一個模型集合M={M1,M2,...,Mk}，k表示模型參數的一個代號，具體對應於一個模型M，代號k的參數是 $heta_k$ 是一個多維參數向量。所以基於觀測值x，（k, $heta_k$ )的後驗分布可以通過貝葉斯方法得到 $pi(k, heta_k|x)=frac{L(x|k, heta_k)p( heta_k|k)p(k)}{sum_{kin K}int_{R^{n_{k}}}L(x|k, heta_{k})p( heta_{k}|k)p(k)d heta_{k}}$ ，分母就是觀測值的非條件概率（基於所有可能的模型和所有可能的參數）。RJ-MCMC就是用這個聯合後驗分布作為MCMC的收斂分布，且該馬爾科夫鏈的狀態空間是 $Theta = igcup_{kin K}(left{ k ight} imes R^{n_k})$ ，即每一個狀態長（k, $heta_k$ ）這樣。因此，通過單個馬爾科夫鏈，可以得到基於觀測值的完整後驗概率，以及單個模型的後驗分布。

具體的實現如下：

給定初始的k和theta_k，for i = 1 to N: 模型內移動：在模型k內，用參考分布抽樣一個新的參數向量theta_k 模型間跳動：根據預設的接受概率，在不同模型間移動，即同時更新k和theta_k

MCMC的前沿研究

收斂速度

希望找到一個界限，使得在該步以後MC基本收斂。

自適應MCMC

在抽樣過程中根據抽取的樣本調整參考分布，不過目前的研究都沒有顯著效果

Sequential Monte Carlo and particle filters

適用於連續型數據分析，例如信號等。

機器學習領域

主要應用在高維度模型，大數據。

----------------------------------以下是附錄----------------------------------

Code 1 Rejection Sampling

import numpy as npimport matplotlib.pyplot as pltimport mathdef p(x): #standard normal mu=0 sigma=1 return 1/(math.pi*2)**0.5/sigma*np.exp(-(x-mu)**2/2/sigma**2)#uniform proposal distribution on [-4,4]def q(x): #uniform return np.array([0.125 for i in range(len(x))])x = np.linspace(-4,4,500)M = 3.5N=1000 #number of samples neededi = 1count = 0X = np.array([])while i < N: u = np.random.rand(10) #evaluate 10 each loop x = (np.random.rand(10)-0.5)*8 res = u < p(x)/q(x)/M if any(res): X = np.hstack((X,x[res])) i+=len(x[res]) count+=10count -= len(X) - 1000X=X[:1000]x = np.linspace(-4,4,500)plt.plot(x,p(x))plt.hist(X,bins=100,normed=True)plt.title(Rejection Sampling)print N/count #proportion of raw sample used

Code 2: Sampling Importance Resampling

import numpy as npimport matplotlib.pyplot as pltdef p(x): #standard normal mu=0 sigma=1 return 1/(math.pi*2)**0.5/sigma*np.exp(-(x-mu)**2/2/sigma**2)#uniform proposal distribution on [-4,4]def q(x): #uniform return np.array([0.125 for i in range(len(x))])#draw N samples that conform to q(x), and then draw M from then that approximately conform to p(x)N=100000M=1000x = (np.random.rand(N)-0.5)*8w_x = p(x)/q(x)w_x = w_x/sum(w_x)w_xc = np.cumsum(w_x) #used for uniform quantile inverse# resample from x with replacement with probability of w_xX=np.array([])for i in range(M): u = np.random.rand() X = np.hstack((X,x[w_xc>u][0]))x = np.linspace(-4,4,500)plt.plot(x,p(x))plt.hist(X,bins=100,normed=True)plt.title(Sampling Importance Resampling)

Code 3: Metropolis_Hastings for MCMC

def q(x_star,x): #normal distribution mu=x sigma=10 return 1/(math.pi*2)**0.5/sigma*np.exp(-(x_star-mu)**2/2/sigma**2) def p(x): #target distribution return 0.3*np.exp(-0.2*x**2)+0.7*np.exp(-0.2*(x-10)**2) N = [100,500,1000,5000]fig = plt.figure()for i in range(4): X = np.array([]) x = 0.1 #initialize x0 to be 0.1 for j in range(N[i]): u = np.random.rand() x_star = np.random.normal(x,10) A = min(1,p(x_star)*q(x,x_star)/p(x)/q(x_star,x)) if u < A: x = x_star X=np.hstack((X,x)) ax = fig.add_subplot(2,2,i+1) ax.hist(X,bins=100,normed=True) x = np.linspace(-10,20,5000) ax.plot(x,p(x)/2.7) #2.7 is just a number that approximates the normalizing constant ax.set_ylim(0,0.35) ax.text(-9,0.25,I=%d%N[i])fig.suptitle(Metropolis_Hastings for MCMC)

Reference

CHRISTOPHE, NANDO, ARNAUD, MICHAEL(2003) An Introduction to MCMC for Machine Learning
Yanan Fan and Scott A. Sisson (2010) Reversible Jump Markov chain MonteCarlo