3.1 隨機模擬

隨機模擬(或者統(tǒng)計模擬)方法有一個很酷的別名是蒙特卡羅方法(Monte Carlo Simulation)。這個方法的發(fā)展始于20世紀40年代，和原子彈制造的曼哈頓計劃密切相關(guān)，當時的幾個大牛，包括烏拉姆、馮.諾依曼、費米、費曼、Nicholas Metropolis， 在美國洛斯阿拉莫斯國家實驗室研究裂變物質(zhì)的中子連鎖反應(yīng)的時候，開始使用統(tǒng)計模擬的方法,并在最早的計算機上進行編程實現(xiàn)。

隨機模擬與計算機

現(xiàn)代的統(tǒng)計模擬方法最早由數(shù)學家烏拉姆提出，被Metropolis命名為蒙特卡羅方法，蒙特卡羅是著名的賭場，賭博總是和統(tǒng)計密切關(guān)聯(lián)的，所以這個命名風趣而貼切，很快被大家廣泛接受。被不過據(jù)說費米之前就已經(jīng)在實驗中使用了，但是沒有發(fā)表。說起蒙特卡羅方法的源頭，可以追溯到18世紀，布豐當年用于計算[Math Processing Error]π的著名的投針實驗就是蒙特卡羅模擬實驗。統(tǒng)計采樣的方法其實數(shù)學家們很早就知道，但是在計算機出現(xiàn)以前，隨機數(shù)生成的成本很高，所以該方法也沒有實用價值。隨著計算機技術(shù)在二十世紀后半葉的迅猛發(fā)展，隨機模擬技術(shù)很快進入實用階段。對那些用確定算法不可行或不可能解決的問題，蒙特卡羅方法常常為人們帶來希望。

蒙特卡羅方法

統(tǒng)計模擬中有一個重要的問題就是給定一個概率分布[Math Processing Error]p(x)，我們?nèi)绾卧谟嬎銠C中生成它的樣本。一般而言均勻分布?[Math Processing Error]Uniform(0,1)的樣本是相對容易生成的。 通過線性同余發(fā)生器可以生成偽隨機數(shù)，我們用確定性算法生成[Math Processing Error][0,1]之間的偽隨機數(shù)序列后，這些序列的各種統(tǒng)計指標和均勻分布?[Math Processing Error]Uniform(0,1)?的理論計算結(jié)果非常接近。這樣的偽隨機序列就有比較好的統(tǒng)計性質(zhì)，可以被當成真實的隨機數(shù)使用。

生成一個概率分布的樣本

而我們常見的概率分布，無論是連續(xù)的還是離散的分布，都可以基于[Math Processing Error]Uniform(0,1)?的樣本生成。例如正態(tài)分布可以通過著名的 Box-Muller 變換得到

[Box-Muller 變換] ?如果隨機變量?[Math Processing Error]U1,U2?獨立且[Math Processing Error]U1,U2～Uniform[0,1]，
[Math Processing Error]Z0=?2ln?U1cos(2πU2)Z1=?2ln?U1sin(2πU2)
則?[Math Processing Error]Z0,Z1?獨立且服從標準正態(tài)分布。

其它幾個著名的連續(xù)分布，包括指數(shù)分布、Gamma 分布、t 分布、F 分布、Beta 分布、Dirichlet 分布等等,也都可以通過類似的數(shù)學變換得到；離散的分布通過均勻分布更加容易生成。更多的統(tǒng)計分布如何通過均勻分布的變換生成出來，大家可以參考統(tǒng)計計算的書，其中 Sheldon M. Ross 的《統(tǒng)計模擬》是寫得非常通俗易懂的一本。

不過我們并不是總是這么幸運的，當[Math Processing Error]p(x)的形式很復(fù)雜，或者[Math Processing Error]p(x)?是個高維的分布的時候，樣本的生成就可能很困難了。 譬如有如下的情況

• [Math Processing Error]p(x)=p~(x)∫p~(x)dx,而?[Math Processing Error]p~(x)?我們是可以計算的，但是底下的積分式無法顯式計算。
• [Math Processing Error]p(x,y)?是一個二維的分布函數(shù)，這個函數(shù)本身計算很困難，但是條件分布?[Math Processing Error]p(x|y),p(y|x)的計算相對簡單;如果?[Math Processing Error]p(x)?是高維的，這種情形就更加明顯。

此時就需要使用一些更加復(fù)雜的隨機模擬的方法來生成樣本。而本節(jié)中將要重點介紹的 MCMC(Markov Chain Monte Carlo) 和 Gibbs Sampling算法就是最常用的一種，這兩個方法在現(xiàn)代貝葉斯分析中被廣泛使用。要了解這兩個算法，我們首先要對馬氏鏈的平穩(wěn)分布的性質(zhì)有基本的認識。

3.2 馬氏鏈及其平穩(wěn)分布

馬氏鏈的數(shù)學定義很簡單
[Math Processing Error]P(Xt+1=x|Xt,Xt?1,?)=P(Xt+1=x|Xt)
也就是狀態(tài)轉(zhuǎn)移的概率只依賴于前一個狀態(tài)。

我們先來看馬氏鏈的一個具體的例子。社會學家經(jīng)常把人按其經(jīng)濟狀況分成3類：下層(lower-class)、中層(middle-class)、上層(upper-class)，我們用1,2,3 分別代表這三個階層。社會學家們發(fā)現(xiàn)決定一個人的收入階層的最重要的因素就是其父母的收入階層。如果一個人的收入屬于下層類別，那么他的孩子屬于下層收入的概率是 0.65, 屬于中層收入的概率是 0.28, 屬于上層收入的概率是 0.07。事實上，從父代到子代，收入階層的變化的轉(zhuǎn)移概率如下

使用矩陣的表示方式，轉(zhuǎn)移概率矩陣記為
[Math Processing Error]P=[0.650.280.070.150.670.180.120.360.52]

假設(shè)當前這一代人處在下層、中層、上層的人的比例是概率分布向量?[Math Processing Error]π0=[π0(1),π0(2),π0(3)]，那么他們的子女的分布比例將是?[Math Processing Error]π1=π0P, 他們的孫子代的分布比例將是?[Math Processing Error]π2=π1P=π0P2, ……, 第[Math Processing Error]n代子孫的收入分布比例將是?[Math Processing Error]πn=πn?1P=π0Pn。

假設(shè)初始概率分布為[Math Processing Error]π0=[0.21,0.68,0.11]，則我們可以計算前[Math Processing Error]n代人的分布狀況如下

我們發(fā)現(xiàn)從第7代人開始，這個分布就穩(wěn)定不變了，這個是偶然的嗎？我們換一個初始概率分布[Math Processing Error]π0=[0.75,0.15,0.1]?試試看，繼續(xù)計算前[Math Processing Error]n代人的分布狀況如下

我們發(fā)現(xiàn)，到第9代人的時候, 分布又收斂了。最為奇特的是，兩次給定不同的初始概率分布，最終都收斂到概率分布?[Math Processing Error]π=[0.286,0.489,0.225]，也就是說收斂的行為和初始概率分布?[Math Processing Error]π0?無關(guān)。這說明這個收斂行為主要是由概率轉(zhuǎn)移矩陣[Math Processing Error]P決定的。我們計算一下?[Math Processing Error]Pn

[Math Processing Error]P20=P21=?=P100=?=[0.2860.4890.2250.2860.4890.2250.2860.4890.225]

我們發(fā)現(xiàn)，當?[Math Processing Error]n?足夠大的時候，這個[Math Processing Error]Pn矩陣的每一行都是穩(wěn)定地收斂到[Math Processing Error]π=[0.286,0.489,0.225]?這個概率分布。自然的，這個收斂現(xiàn)象并非是我們這個馬氏鏈獨有的，而是絕大多數(shù)馬氏鏈的共同行為，關(guān)于馬氏鏈的收斂我們有如下漂亮的定理：

馬氏鏈定理：?如果一個非周期馬氏鏈具有轉(zhuǎn)移概率矩陣[Math Processing Error]P,且它的任何兩個狀態(tài)是連通的，那么?[Math Processing Error]limn→∞Pijn?存在且與[Math Processing Error]i無關(guān)，記?[Math Processing Error]limn→∞Pijn=π(j), 我們有

[Math Processing Error]limn→∞Pn=[π(1)π(2)?π(j)?π(1)π(2)?π(j)??????π(1)π(2)?π(j)??????]

?[Math Processing Error]π(j)=∑i=0∞π(i)Pij

?[Math Processing Error]π?是方程?[Math Processing Error]πP=π?的唯一非負解

其中,
[Math Processing Error]π=[π(1),π(2),?,π(j),?],∑i=0∞πi=1
[Math Processing Error]π稱為馬氏鏈的平穩(wěn)分布。

這個馬氏鏈的收斂定理非常重要，所有的 MCMC(Markov Chain Monte Carlo) 方法都是以這個定理作為理論基礎(chǔ)的。?定理的證明相對復(fù)雜，一般的隨機過程課本中也不給證明，所以我們就不用糾結(jié)它的證明了，直接用這個定理的結(jié)論就好了。我們對這個定理的內(nèi)容做一些解釋說明：

該定理中馬氏鏈的狀態(tài)不要求有限，可以是有無窮多個的；

定理中的“非周期“這個概念我們不打算解釋了，因為我們遇到的絕大多數(shù)馬氏鏈都是非周期的；

兩個狀態(tài)[Math Processing Error]i,j是連通并非指[Math Processing Error]i?可以直接一步轉(zhuǎn)移到[Math Processing Error]j([Math Processing Error]Pij>0),而是指?[Math Processing Error]i?可以通過有限的[Math Processing Error]n步轉(zhuǎn)移到達[Math Processing Error]j([Math Processing Error]Pijn>0)。馬氏鏈的任何兩個狀態(tài)是連通的含義是指存在一個[Math Processing Error]n, 使得矩陣[Math Processing Error]Pn?中的任何一個元素的數(shù)值都大于零。

我們用?[Math Processing Error]Xi?表示在馬氏鏈上跳轉(zhuǎn)第[Math Processing Error]i步后所處的狀態(tài)，如果?[Math Processing Error]limn→∞Pijn=π(j)?存在，很容易證明以上定理的第二個結(jié)論。由于
[Math Processing Error]P(Xn+1=j)=∑i=0∞P(Xn=i)P(Xn+1=j|Xn=i)=∑i=0∞P(Xn=i)Pij
上式兩邊取極限就得到?[Math Processing Error]π(j)=∑i=0∞π(i)Pij

從初始概率分布?[Math Processing Error]π0?出發(fā)，我們在馬氏鏈上做狀態(tài)轉(zhuǎn)移，記[Math Processing Error]Xi的概率分布為[Math Processing Error]πi, 則有
[Math Processing Error]X0～π0(x)Xi～πi(x),πi(x)=πi?1(x)P=π0(x)Pn
由馬氏鏈收斂的定理, 概率分布[Math Processing Error]πi(x)將收斂到平穩(wěn)分布?[Math Processing Error]π(x)。假設(shè)到第[Math Processing Error]n步的時候馬氏鏈收斂，則有
[Math Processing Error]X0～π0(x)X1～π1(x)?Xn～πn(x)=π(x)Xn+1～π(x)Xn+2～π(x)?
所以?[Math Processing Error]Xn,Xn+1,Xn+2,?～π(x)?都是同分布的隨機變量，當然他們并不獨立。如果我們從一個具體的初始狀態(tài)?[Math Processing Error]x0?開始,沿著馬氏鏈按照概率轉(zhuǎn)移矩陣做跳轉(zhuǎn)，那么我們得到一個轉(zhuǎn)移序列?[Math Processing Error]x0,x1,x2,?xn,xn+1?,?由于馬氏鏈的收斂行為，[Math Processing Error]xn,xn+1,??都將是平穩(wěn)分布?[Math Processing Error]π(x)?的樣本。

3.3 Markov Chain Monte Carlo

對于給定的概率分布[Math Processing Error]p(x),我們希望能有便捷的方式生成它對應(yīng)的樣本。由于馬氏鏈能收斂到平穩(wěn)分布， 于是一個很的漂亮想法是：如果我們能構(gòu)造一個轉(zhuǎn)移矩陣為[Math Processing Error]P的馬氏鏈，使得該馬氏鏈的平穩(wěn)分布恰好是[Math Processing Error]p(x), 那么我們從任何一個初始狀態(tài)[Math Processing Error]x0出發(fā)沿著馬氏鏈轉(zhuǎn)移, 得到一個轉(zhuǎn)移序列?[Math Processing Error]x0,x1,x2,?xn,xn+1?,， 如果馬氏鏈在第[Math Processing Error]n步已經(jīng)收斂了，于是我們就得到了?[Math Processing Error]π(x)?的樣本[Math Processing Error]xn,xn+1?。

這個絕妙的想法在1953年被 Metropolis想到了，為了研究粒子系統(tǒng)的平穩(wěn)性質(zhì)， Metropolis 考慮了物理學中常見的波爾茲曼分布的采樣問題，首次提出了基于馬氏鏈的蒙特卡羅方法，即Metropolis算法，并在最早的計算機上編程實現(xiàn)。Metropolis 算法是首個普適的采樣方法，并啟發(fā)了一系列 MCMC方法，所以人們把它視為隨機模擬技術(shù)騰飛的起點。 Metropolis的這篇論文被收錄在《統(tǒng)計學中的重大突破》中， Metropolis算法也被遴選為二十世紀的十個最重要的算法之一。

我們接下來介紹的MCMC 算法是 Metropolis 算法的一個改進變種，即常用的 Metropolis-Hastings 算法。由上一節(jié)的例子和定理我們看到了，馬氏鏈的收斂性質(zhì)主要由轉(zhuǎn)移矩陣[Math Processing Error]P?決定, 所以基于馬氏鏈做采樣的關(guān)鍵問題是如何構(gòu)造轉(zhuǎn)移矩陣[Math Processing Error]P,使得平穩(wěn)分布恰好是我們要的分布[Math Processing Error]p(x)。如何能做到這一點呢？我們主要使用如下的定理。

定理：[細致平穩(wěn)條件]?如果非周期馬氏鏈的轉(zhuǎn)移矩陣[Math Processing Error]P和分布[Math Processing Error]π(x)?滿足
[Math Processing Error](1)π(i)Pij=π(j)Pjifor alli,j
則?[Math Processing Error]π(x)?是馬氏鏈的平穩(wěn)分布，上式被稱為細致平穩(wěn)條件(detailed balance condition)。

其實這個定理是顯而易見的，因為細致平穩(wěn)條件的物理含義就是對于任何兩個狀態(tài)[Math Processing Error]i,j, 從?[Math Processing Error]i?轉(zhuǎn)移出去到[Math Processing Error]j?而丟失的概率質(zhì)量，恰好會被從?[Math Processing Error]j?轉(zhuǎn)移回[Math Processing Error]i?的概率質(zhì)量補充回來，所以狀態(tài)[Math Processing Error]i上的概率質(zhì)量[Math Processing Error]π(i)是穩(wěn)定的，從而[Math Processing Error]π(x)是馬氏鏈的平穩(wěn)分布。數(shù)學上的證明也很簡單，由細致平穩(wěn)條件可得
[Math Processing Error]∑i=1∞π(i)Pij=∑i=1∞π(j)Pji=π(j)∑i=1∞Pji=π(j)?πP=π
由于[Math Processing Error]π?是方程?[Math Processing Error]πP=π的解，所以[Math Processing Error]π是平穩(wěn)分布。

假設(shè)我們已經(jīng)有一個轉(zhuǎn)移矩陣為[Math Processing Error]Q馬氏鏈([Math Processing Error]q(i,j)表示從狀態(tài)?[Math Processing Error]i轉(zhuǎn)移到狀態(tài)[Math Processing Error]j的概率，也可以寫為?[Math Processing Error]q(j|i)或者[Math Processing Error]q(i→j)), 顯然，通常情況下?[Math Processing Error]p(i)q(i,j)≠p(j)q(j,i)?也就是細致平穩(wěn)條件不成立，所以?[Math Processing Error]p(x)?不太可能是這個馬氏鏈的平穩(wěn)分布。我們可否對馬氏鏈做一個改造，使得細致平穩(wěn)條件成立呢？譬如，我們引入一個?[Math Processing Error]α(i,j), 我們希望
[Math Processing Error](2)p(i)q(i,j)α(i,j)=p(j)q(j,i)α(j,i)(?)
取什么樣的?[Math Processing Error]α(i,j)?以上等式能成立呢？最簡單的，按照對稱性，我們可以取
[Math Processing Error]α(i,j)=p(j)q(j,i)，α(j,i)=p(i)q(i,j)
于是(*)式就成立了。所以有
[Math Processing Error](3)p(i)q(i,j)α(i,j)?Q′(i,j)=p(j)q(j,i)α(j,i)?Q′(j,i)(??)
于是我們把原來具有轉(zhuǎn)移矩陣[Math Processing Error]Q的一個很普通的馬氏鏈，改造為了具有轉(zhuǎn)移矩陣[Math Processing Error]Q′的馬氏鏈，而?[Math Processing Error]Q′恰好滿足細致平穩(wěn)條件，由此馬氏鏈[Math Processing Error]Q′的平穩(wěn)分布就是[Math Processing Error]p(x)！

在改造?[Math Processing Error]Q?的過程中引入的?[Math Processing Error]α(i,j)稱為接受率，物理意義可以理解為在原來的馬氏鏈上，從狀態(tài)?[Math Processing Error]i?以[Math Processing Error]q(i,j)?的概率轉(zhuǎn)跳轉(zhuǎn)到狀態(tài)[Math Processing Error]j?的時候，我們以[Math Processing Error]α(i,j)的概率接受這個轉(zhuǎn)移，于是得到新的馬氏鏈[Math Processing Error]Q′的轉(zhuǎn)移概率為[Math Processing Error]q(i,j)α(i,j)。

馬氏鏈轉(zhuǎn)移和接受概率

假設(shè)我們已經(jīng)有一個轉(zhuǎn)移矩陣Q(對應(yīng)元素為[Math Processing Error]q(i,j)), 把以上的過程整理一下，我們就得到了如下的用于采樣概率分布[Math Processing Error]p(x)的算法。

上述過程中?[Math Processing Error]p(x),q(x|y)?說的都是離散的情形，事實上即便這兩個分布是連續(xù)的，以上算法仍然是有效，于是就得到更一般的連續(xù)概率分布?[Math Processing Error]p(x)的采樣算法，而?[Math Processing Error]q(x|y)?就是任意一個連續(xù)二元概率分布對應(yīng)的條件分布。

以上的 MCMC 采樣算法已經(jīng)能很漂亮的工作了，不過它有一個小的問題：馬氏鏈[Math Processing Error]Q在轉(zhuǎn)移的過程中的接受率?[Math Processing Error]α(i,j)?可能偏小，這樣采樣過程中馬氏鏈容易原地踏步，拒絕大量的跳轉(zhuǎn)，這使得馬氏鏈遍歷所有的狀態(tài)空間要花費太長的時間，收斂到平穩(wěn)分布[Math Processing Error]p(x)的速度太慢。有沒有辦法提升一些接受率呢?

假設(shè)?[Math Processing Error]α(i,j)=0.1,α(j,i)=0.2, 此時滿足細致平穩(wěn)條件，于是
[Math Processing Error]p(i)q(i,j)×0.1=p(j)q(j,i)×0.2
上式兩邊擴大5倍，我們改寫為
[Math Processing Error]p(i)q(i,j)×0.5=p(j)q(j,i)×1
看，我們提高了接受率，而細致平穩(wěn)條件并沒有打破！這啟發(fā)我們可以把細致平穩(wěn)條件(**) 式中的[Math Processing Error]α(i,j),α(j,i)?同比例放大，使得兩數(shù)中最大的一個放大到1，這樣我們就提高了采樣中的跳轉(zhuǎn)接受率。所以我們可以取
[Math Processing Error]α(i,j)=min{p(j)q(j,i)p(i)q(i,j),1}
于是，經(jīng)過對上述MCMC 采樣算法中接受率的微小改造，我們就得到了如下教科書中最常見的 Metropolis-Hastings 算法。

對于分布?[Math Processing Error]p(x),我們構(gòu)造轉(zhuǎn)移矩陣?[Math Processing Error]Q′?使其滿足細致平穩(wěn)條件
[Math Processing Error]p(x)Q′(x→y)=p(y)Q′(y→x)
此處?[Math Processing Error]x?并不要求是一維的，對于高維空間的?[Math Processing Error]p(x)，如果滿足細致平穩(wěn)條件
[Math Processing Error]p(x)Q′(x→y)=p(y)Q′(y→x)
那么以上的 Metropolis-Hastings 算法一樣有效。

3.2 Gibbs Sampling

對于高維的情形，由于接受率?[Math Processing Error]α的存在(通常?[Math Processing Error]α<1), 以上 Metropolis-Hastings 算法的效率不夠高。能否找到一個轉(zhuǎn)移矩陣Q使得接受率?[Math Processing Error]α=1?呢？我們先看看二維的情形，假設(shè)有一個概率分布[Math Processing Error]p(x,y), 考察[Math Processing Error]x坐標相同的兩個點[Math Processing Error]A(x1,y1),B(x1,y2)，我們發(fā)現(xiàn)
[Math Processing Error]p(x1,y1)p(y2|x1)=p(x1)p(y1|x1)p(y2|x1)p(x1,y2)p(y1|x1)=p(x1)p(y2|x1)p(y1|x1)
所以得到
[Math Processing Error](4)p(x1,y1)p(y2|x1)=p(x1,y2)p(y1|x1)(???)
即
[Math Processing Error]p(A)p(y2|x1)=p(B)p(y1|x1)
基于以上等式，我們發(fā)現(xiàn)，在?[Math Processing Error]x=x1?這條平行于?[Math Processing Error]y軸的直線上，如果使用條件分布?[Math Processing Error]p(y|x1)做為任何兩個點之間的轉(zhuǎn)移概率，那么任何兩個點之間的轉(zhuǎn)移滿足細致平穩(wěn)條件。同樣的，如果我們在[Math Processing Error]y=y1?這條直線上任意取兩個點?[Math Processing Error]A(x1,y1),C(x2,y1),也有如下等式
[Math Processing Error]p(A)p(x2|y1)=p(C)p(x1|y1).

平面上馬氏鏈轉(zhuǎn)移矩陣的構(gòu)造

于是我們可以如下構(gòu)造平面上任意兩點之間的轉(zhuǎn)移概率矩陣Q
[Math Processing Error]Q(A→B)=p(yB|x1)如果xA=xB=x1Q(A→C)=p(xC|y1)如果yA=yC=y1Q(A→D)=0其它

有了如上的轉(zhuǎn)移矩陣 Q, 我們很容易驗證對平面上任意兩點?[Math Processing Error]X,Y, 滿足細致平穩(wěn)條件
[Math Processing Error]p(X)Q(X→Y)=p(Y)Q(Y→X)
于是這個二維空間上的馬氏鏈將收斂到平穩(wěn)分布?[Math Processing Error]p(x,y)。而這個算法就稱為 Gibbs Sampling 算法,是 Stuart Geman 和Donald Geman 這兩兄弟于1984年提出來的，之所以叫做Gibbs Sampling 是因為他們研究了Gibbs random field, 這個算法在現(xiàn)代貝葉斯分析中占據(jù)重要位置。

Gibbs Sampling 算法中的馬氏鏈轉(zhuǎn)移

以上采樣過程中，如圖所示，馬氏鏈的轉(zhuǎn)移只是輪換的沿著坐標軸?[Math Processing Error]x軸和[Math Processing Error]y軸做轉(zhuǎn)移，于是得到樣本?[Math Processing Error](x0,y0),(x0,y1),(x1,y1),(x1,y2),(x2,y2),?馬氏鏈收斂后，最終得到的樣本就是?[Math Processing Error]p(x,y)?的樣本，而收斂之前的階段稱為 burn-in period。額外說明一下，我們看到教科書上的 Gibbs Sampling 算法大都是坐標軸輪換采樣的，但是這其實是不強制要求的。最一般的情形可以是，在[Math Processing Error]t時刻，可以在[Math Processing Error]x軸和[Math Processing Error]y軸之間隨機的選一個坐標軸，然后按條件概率做轉(zhuǎn)移，馬氏鏈也是一樣收斂的。輪換兩個坐標軸只是一種方便的形式。

以上的過程我們很容易推廣到高維的情形，對于(***) 式，如果[Math Processing Error]x1變?yōu)槎嗑S情形[Math Processing Error]x1，可以看出推導(dǎo)過程不變，所以細致平穩(wěn)條件同樣是成立的
[Math Processing Error](5)p(x1,y1)p(y2|x1)=p(x1,y2)p(y1|x1)
此時轉(zhuǎn)移矩陣 Q 由條件分布?[Math Processing Error]p(y|x1)?定義。上式只是說明了一根坐標軸的情形，和二維情形類似，很容易驗證對所有坐標軸都有類似的結(jié)論。所以[Math Processing Error]n維空間中對于概率分布?[Math Processing Error]p(x1,x2,?,xn)?可以如下定義轉(zhuǎn)移矩陣

如果當前狀態(tài)為[Math Processing Error](x1,x2,?,xn)，馬氏鏈轉(zhuǎn)移的過程中，只能沿著坐標軸做轉(zhuǎn)移。沿著?[Math Processing Error]xi?這根坐標軸做轉(zhuǎn)移的時候，轉(zhuǎn)移概率由條件概率?[Math Processing Error]p(xi|x1,?,xi?1,xi+1,?,xn)?定義；

?其它無法沿著單根坐標軸進行的跳轉(zhuǎn)，轉(zhuǎn)移概率都設(shè)置為 0。

于是我們可以把Gibbs Smapling 算法從采樣二維的?[Math Processing Error]p(x,y)?推廣到采樣[Math Processing Error]n?維的?[Math Processing Error]p(x1,x2,?,xn)

以上算法收斂后，得到的就是概率分布[Math Processing Error]p(x1,x2,?,xn)的樣本，當然這些樣本并不獨立，但是我們此處要求的是采樣得到的樣本符合給定的概率分布，并不要求獨立。同樣的，在以上算法中，坐標軸輪換采樣不是必須的，可以在坐標軸輪換中引入隨機性，這時候轉(zhuǎn)移矩陣?[Math Processing Error]Q?中任何兩個點的轉(zhuǎn)移概率中就會包含坐標軸選擇的概率，而在通常的 Gibbs Sampling 算法中，坐標軸輪換是一個確定性的過程，也就是在給定時刻[Math Processing Error]t，在一根固定的坐標軸上轉(zhuǎn)移的概率是1。

轉(zhuǎn)載地址：https://cos.name/2013/01/lda-math-mcmc-and-gibbs-sampling/

這篇鏈接http://www.myexception.cn/cloud/2100158.html也值得學習。

總結(jié)

以上是生活随笔為你收集整理的(转载)机器学习知识点(十三)吉布斯采样法(Gibbs Sampling)的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。