Ⅰ MCMC基本原理與應用(一)
貝葉斯統計推斷是一個非常 藝術 的東西,他的先驗結果(prior knowledge)是一種專家的經驗,而 貝葉斯公式 給出的後驗分布是貝葉斯推斷的基本工具。每個大學生都應該學過數理統計,其中提及的貝葉斯統計推斷都是一些基於簡單先驗和簡單後驗的結果。
但當我們一旦碰到復雜的先驗(比如高維參數、復雜先驗分布),我們用貝葉斯公式得出的後驗分布將變得非常復雜,在計算上會非常困難。為此,先人提出了 MCMC演算法 方便我們可以對任何後驗分布進行計算或推斷。其思想就是其名字:兩個MC。
第一個MC: Monte Carlo(蒙特卡洛)。這個簡單來說是讓我們使用隨機數(隨機抽樣)來解決計算問題。在MCMC中意味著:後驗分布作為一個隨機樣本生成器,我們利用它來生成樣本(simulation),然後通過這些樣本對一些感興趣的計算問題(特徵數,預測)進行估計。
第二個MC:Markov Chain(馬爾科夫鏈)。第二個MC是這個方法的關鍵,因為我們在第一個MC中看到,我們需要利用 後驗分布 生成隨機樣本,但後驗分布太復雜,一些package中根本沒有相應的隨機數生成函數(如 rnorm(),rbinom() )怎麼辦?答案是我們可以利用Markov Chain的 平穩分布 這個概念實現對復雜後驗分布的抽樣。
為了能順利闡述MCMC演算法,這里就簡單地講一下所涉及到的馬爾科夫鏈概念。因為都是我個人的理解,這里所講不敢說都是正確的,希望各位童鞋能夠獨立思考。
隨機過程${X_n},X_n$的狀態空間為有限集或可列集,比如當$X_n=i$,就稱過程在n處於狀態i。
定義:如果對任何一列狀態$i_0,i_1,...,i_{n-1},i,j$,及對任何n≥0,隨機過程${X_n}$滿足Markov性質:
$$P{X_{n+1}=j|X_0=i_0,...X_{n-1}=i_{n-1},X_n=i}=P{X_{n+1}=j|X_n=i}$$
轉移概率
$P{x_{n+1}=j|X_n=i}$成為Markov鏈的一步轉移概率$P_{ij}^{n,n+1}$,當這個概率與n無關,則記之為$P_{ij}$
轉移概率矩陣P
這個矩陣的元素就是一步轉移概率$P_{ij}$
結論
一個Markov鏈可以由它的初始狀態以及轉移概率矩陣P完全確定。(證明略,自行網路或翻書)
所謂n步轉移概率就是從狀態i走n步正好到狀態j的概率,我們記為$P_{ij}^{(n)}$。
利用概率分割的思想,由基礎概率論中 全概率公式 可以得到$$P_{ij} {(n)}=Sigma_{k=0} {infty}P_{ik}P_{kj}^{(n-1)}$$ 寫成矩陣形式就是$$P^{(n)}=P imes P^{(n-1)}$$
進一步推廣,我們就推出了著名的Chapman-Kolmogorov方程:$$P_{ij} {(n+m)}=Sigma_{k=0} {infty}P_{ik} {(n)}P_{kj} {(m)}$$寫成矩陣形式就是$$P {(m+n)}=P {(m)} imes P^{(n)}$$
現在我們已經有了n步轉移概率的概念,一個很簡單的想法就是如果n趨向無窮會怎麼樣?這個問題也是後面 極限分布 以及 平穩分布 的來源,更是MCMC演算法中第二個MC的關鍵。
要回答這個問題首先要掌握幾個關鍵概念,我先列出來,如果不熟悉的可以自行網路或翻書:
幾個重要結論(證明略,自行網路,或者call me)
是時候回答上面那個問題了!(摘自網路共享PPT)
另外對於一個正常返、非周期的狀態(也稱為遍歷ergodic),我們有結論:$$對所有ileftrightarrow j,有limP_{ji} {(n)}=limP_{ii} {(n)}=frac{1}{mu_i}$$
什麼是平穩分布?它和求極限概率分布有什麼關系呢?
定義:Markov鏈有轉移概率矩陣P,如果有一個概率分布${pi_i }滿足pi_j =Sigma_{i=0}^{infty}pi_i P_{ij}$,則稱為這個Markov鏈的平穩分布。這個定義用矩陣形式寫出來就是π*P=π.
重要定理
若一個markov鏈中所有的狀態都是遍歷的,則對所有i,j有$limP_{ij} {(n)}=pi_j存在,且pi={pi_j,j≥0}就是平穩分布$。反之,拖一個不可約Markov鏈只存在一個平穩分布,且這個Markov鏈的所有狀態都是遍歷的,則這個平穩分布就是該Markov鏈的極限概率分布。$$limP_{ij} {(n)}=pi_j$$