天天看點

MCMC(一)蒙特卡羅方法1. MCMC概述2. 蒙特卡羅方法引入3. 機率分布采樣4. 接受-拒絕采樣5. 蒙特卡羅方法小結

    從名字我們可以看出,MCMC由兩個MC組成,即蒙特卡羅方法(Monte Carlo Simulation,簡稱MC)和馬爾科夫鍊(Markov Chain ,也簡稱MC)。要弄懂MCMC的原理我們首先得搞清楚蒙特卡羅方法和馬爾科夫鍊的原理。我們将用三篇來完整學習MCMC。在本篇,我們關注于蒙特卡羅方法。

    蒙特卡羅原來是一個賭場的名稱,用它作為名字大概是因為蒙特卡羅方法是一種随機模拟的方法,這很像賭博場裡面的扔骰子的過程。最早的蒙特卡羅方法都是為了求解一些不太好求解的求和或者積分問題。比如積分:

θ=∫baf(x)dxθ=∫abf(x)dx

    如果我們很難求解出f(x)f(x)的原函數,那麼這個積分比較難求解。當然我們可以通過蒙特卡羅方法來模拟求解近似值。如何模拟呢?假設我們函數圖像如下圖:

MCMC(一)蒙特卡羅方法1. MCMC概述2. 蒙特卡羅方法引入3. 機率分布采樣4. 接受-拒絕采樣5. 蒙特卡羅方法小結

    則一個簡單的近似求解方法是在[a,b]之間随機的采樣一個點。比如x0x0,然後用f(x0)f(x0)代表在[a,b]區間上所有的f(x)f(x)的值。那麼上面的定積分的近似求解為:

(b−a)f(x0)(b−a)f(x0)

    當然,用一個值代表[a,b]區間上所有的f(x)f(x)的值,這個假設太粗糙。那麼我們可以采樣[a,b]區間的n個值:x0,x1,...xn−1x0,x1,...xn−1,用它們的均值來代表[a,b]區間上所有的f(x)f(x)的值。這樣我們上面的定積分的近似求解為:

b−an∑i=0n−1f(xi)b−an∑i=0n−1f(xi)

    雖然上面的方法可以一定程度上求解出近似的解,但是它隐含了一個假定,即xx在[a,b]之間是均勻分布的,而絕大部分情況,xx在[a,b]之間不是均勻分布的。如果我們用上面的方法,則模拟求出的結果很可能和真實值相差甚遠。 

    怎麼解決這個問題呢? 如果我們可以得到xx在[a,b]的機率分布函數p(x)p(x),那麼我們的定積分求和可以這樣進行:

θ=∫baf(x)dx=∫baf(x)p(x)p(x)dx≈1n∑i=0n−1f(xi)p(xi)θ=∫abf(x)dx=∫abf(x)p(x)p(x)dx≈1n∑i=0n−1f(xi)p(xi)

    上式最右邊的這個形式就是蒙特卡羅方法的一般形式。當然這裡是連續函數形式的蒙特卡羅方法,但是在離散時一樣成立。

    可以看出,最上面我們假設xx在[a,b]之間是均勻分布的時候,p(xi)=1/(b−a)p(xi)=1/(b−a),帶入我們有機率分布的蒙特卡羅積分的上式,可以得到:

1n∑i=0n−1f(xi)1/(b−a)=b−an∑i=0n−1f(xi)1n∑i=0n−1f(xi)1/(b−a)=b−an∑i=0n−1f(xi)

    也就是說,我們最上面的均勻分布也可以作為一般機率分布函數p(x)p(x)在均勻分布時候的特例。那麼我們現在的問題轉到了如何求出xx的分布p(x)p(x)的若幹和樣本上來。

    上一節我們講到蒙特卡羅方法的關鍵是得到xx的機率分布。如果求出了xx的機率分布,我們可以基于機率分布去采樣基于這個機率分布的n個xx的樣本集,帶入蒙特卡羅求和的式子即可求解。但是還有一個關鍵的問題需要解決,即如何基于機率分布去采樣基于這個機率分布的n個xx的樣本集。 

    對于常見的均勻分布uniform(0,1)uniform(0,1)是非常容易采樣樣本的,一般通過線性同餘發生器可以很友善的生成(0,1)之間的僞随機數樣本。而其他常見的機率分布,無論是離散的分布還是連續的分布,它們的樣本都可以通過uniform(0,1)uniform(0,1)的樣本轉換而得。比如二維正态分布的樣本(Z1,Z2)(Z1,Z2)可以通過通過獨立采樣得到的uniform(0,1)uniform(0,1)樣本對(X1,X2)(X1,X2)通過如下的式子轉換而得:

Z1=−2lnX1−−−−−−−√cos(2πX2)Z1=−2lnX1cos(2πX2)

Z2=−2lnX1−−−−−−−√sin(2πX2)Z2=−2lnX1sin(2πX2)

    其他一些常見的連續分布,比如t分布,F分布,Beta分布,Gamma分布等,都可以通過類似的方式從uniform(0,1)uniform(0,1)得到的采樣樣本轉化得到。在python的numpy,scikit-learn等類庫中,都有生成這些常用分布樣本的函數可以使用。

    不過很多時候,我們的xx的機率分布不是常見的分布,這意味着我們沒法友善的得到這些非常見的機率分布的樣本集。那這個問題怎麼解決呢?

    對于機率分布不是常見的分布,一個可行的辦法是采用接受-拒絕采樣來得到該分布的樣本。既然 p(x)p(x) 太複雜在程式中沒法直接采樣,那麼我設定一個程式可采樣的分布 q(x)q(x) 比如高斯分布,然後按照一定的方法拒絕某些樣本,以達到接近 p(x)p(x) 分布的目的,其中q(x)q(x)叫做 proposal distribution。

MCMC(一)蒙特卡羅方法1. MCMC概述2. 蒙特卡羅方法引入3. 機率分布采樣4. 接受-拒絕采樣5. 蒙特卡羅方法小結

    具體采用過程如下,設定一個友善采樣的常用機率分布函數 q(x)q(x),以及一個常量 kk,使得 p(x)p(x) 總在 kq(x)kq(x) 的下方。如上圖。

    首先,采樣得到q(x)q(x)的一個樣本z0z0,采樣方法如第三節。然後,從均勻分布(0,kq(z0))(0,kq(z0))中采樣得到一個值uu。如果uu落在了上圖中的灰色區域,則拒絕這次抽樣,否則接受這個樣本z0z0。重複以上過程得到n個接受的樣本z0,z1,...zn−1z0,z1,...zn−1,則最後的蒙特卡羅方法求解結果為:

1n∑i=0n−1f(zi)p(zi)1n∑i=0n−1f(zi)p(zi)

    整個過程中,我們通過一系列的接受拒絕決策來達到用q(x)q(x)模拟p(x)p(x)機率分布的目的。

    使用接受-拒絕采樣,我們可以解決一些機率分布不是常見的分布的時候,得到其采樣集并用蒙特卡羅方法求和的目的。但是接受-拒絕采樣也隻能部分滿足我們的需求,在很多時候我們還是很難得到我們的機率分布的樣本集。比如:

    1)對于一些二維分布p(x,y)p(x,y),有時候我們隻能得到條件分布p(x|y)p(x|y)和p(y|x)p(y|x)和,卻很難得到二維分布p(x,y)p(x,y)一般形式,這時我們無法用接受-拒絕采樣得到其樣本集。

    2)對于一些高維的複雜非常見分布p(x1,x2,...,xn)p(x1,x2,...,xn),我們要找到一個合适的q(x)q(x)和kk非常困難。

    從上面可以看出,要想将蒙特卡羅方法作為一個通用的采樣模拟求和的方法,必須解決如何友善得到各種複雜機率分布的對應的采樣樣本集的問題。而我們下一篇要講到的馬爾科夫鍊就是幫助找到這些複雜機率分布的對應的采樣樣本集的白衣騎士。下一篇我們來總結馬爾科夫鍊的原理。

本文轉自劉建平Pinard部落格園部落格,原文連結:http://www.cnblogs.com/pinard/p/6625739.html,如需轉載請自行聯系原作者