2.2　概率圖模型

概率圖模型是通過(guò)以圖為工具表達(dá)變量關(guān)系的概率模型，一般一個(gè)節(jié)點(diǎn)表示一個(gè)或一組隨機(jī)變量，變量之間的概率關(guān)系通過(guò)節(jié)點(diǎn)間的邊表示。概率圖模型可以通過(guò)有向無(wú)環(huán)圖和無(wú)向圖表示變量間的依賴關(guān)系，前者稱為有貝葉斯網(wǎng)或有向圖模型，后者稱為馬爾可夫網(wǎng)或無(wú)向圖模型。

2.2.1　隱馬爾可夫模型

馬爾可夫過(guò)程表示數(shù)學(xué)中具有馬爾可夫性質(zhì)的離散隨機(jī)過(guò)程，在該過(guò)程中，下一時(shí)刻的狀態(tài)只與當(dāng)前狀態(tài)有關(guān)，與上一時(shí)刻狀態(tài)無(wú)關(guān)。最簡(jiǎn)單的馬爾可夫過(guò)程就是一階過(guò)程，每一個(gè)狀態(tài)的轉(zhuǎn)移只依賴其之前的那一個(gè)狀態(tài)。

假定一天的天氣狀況有晴（高溫狀態(tài)）、陰（低溫狀態(tài)）、雨（涼爽狀態(tài)）三種，這三件事按圖2.1所示的方向轉(zhuǎn)移：

這個(gè)過(guò)程就是一個(gè)簡(jiǎn)單的馬爾可夫過(guò)程，但馬爾可夫過(guò)程和確定性系統(tǒng)不同，狀態(tài)之間的轉(zhuǎn)移是有概率的，天氣的狀態(tài)是經(jīng)常變化的，而且會(huì)在兩個(gè)狀態(tài)間任意切換，如圖2.2所示。

圖2.2中箭頭表示天氣從一個(gè)狀態(tài)切換到另一個(gè)狀態(tài)的概率，保持晴天的概率是0.5，由晴天轉(zhuǎn)為陰天的概率為0.4，由晴天轉(zhuǎn)為下雨的概率為0.1。

當(dāng)前狀態(tài)的轉(zhuǎn)移依賴之前的n個(gè)狀態(tài)，當(dāng)n為1時(shí)即為馬爾可夫假設(shè)；并由此可以得出馬爾可夫鏈的如下定義：

馬爾可夫鏈?zhǔn)且唤M具有馬爾可夫性質(zhì)的離散隨機(jī)變量的集合S=｛St∶t＞0｝，隨機(jī)變量的所有可能取值的集合被稱為狀態(tài)空間，St的值則是在時(shí)間t的狀態(tài)。若St+1對(duì)于過(guò)去狀態(tài)的條件概率分布僅是St的一個(gè)函數(shù)，即稱S為馬爾可夫鏈。

這里定義的馬爾可夫鏈?zhǔn)请x散時(shí)間馬爾可夫鏈，具有連續(xù)指數(shù)集的情形稱為馬爾可夫過(guò)程。

隱馬爾可夫模型是描述含有隱含未知參數(shù)的馬爾可夫過(guò)程，從可觀察的參數(shù)中確定隱含參數(shù)，然后利用這些參數(shù)來(lái)做進(jìn)一步的統(tǒng)計(jì)與預(yù)測(cè)，是最簡(jiǎn)單的貝葉斯網(wǎng)，多用于動(dòng)態(tài)時(shí)序數(shù)據(jù)建模，被建模的系統(tǒng)被認(rèn)為是一個(gè)馬爾可夫過(guò)程與未觀測(cè)到的（隱藏的）狀態(tài)統(tǒng)計(jì)。

隱馬爾可夫模型包含兩個(gè)狀態(tài)集合，分別為可觀測(cè)的狀態(tài)集和隱藏的狀態(tài)集，通過(guò)可觀測(cè)狀態(tài)預(yù)測(cè)隱藏狀態(tài)，如圖2.3所示。

由此可以通過(guò)觀測(cè)狀態(tài)溫度、大氣壓和降水等觀測(cè)狀態(tài)O=｛OT，OA，OR｝來(lái)判斷隱藏狀態(tài)S=｛SF，SC，SR｝。在晴天狀態(tài)下，觀測(cè)到溫度高的概率高，觀測(cè)到氣壓低的概率低，觀測(cè)到降水的概率低；如果觀測(cè)到溫度低的概率高，氣壓低的概率高，降水的概率較低的時(shí)候，則大概率是陰天的狀態(tài)；如果觀測(cè)到溫度低的概率低，氣壓低的概率低，降水的概率高的時(shí)候，則大概率是下雨的狀態(tài)。由此可知，天氣可觀測(cè)的狀態(tài)序列和隱藏的狀態(tài)序列是概率相關(guān)的，則可以將此類型的過(guò)程建模為一個(gè)隱藏的馬爾可夫過(guò)程和一個(gè)與之概率相關(guān)的可觀測(cè)狀態(tài)集合，即隱馬爾可夫模型。

隱馬爾可夫模型可以通過(guò)混淆矩陣表示，矩陣的行表示隱藏狀態(tài)，每一行的概率值之和為1，列表示可觀測(cè)狀態(tài)，圖2.3的混淆矩陣為：

上述問(wèn)題的狀態(tài)轉(zhuǎn)移矩陣為：

隱馬爾可夫模型可以用5元組｛N，M，π，A，B｝表示，N表示隱藏狀態(tài)的數(shù)量，M表示可觀測(cè)狀態(tài)的數(shù)量，π=｛πi｝表示初始隱藏狀態(tài)的概率，A表示隱藏狀態(tài)的轉(zhuǎn)移矩陣，B表示混淆矩陣。

圖2.3中隱藏狀態(tài)N=3，可觀測(cè)狀態(tài)M=3。

在實(shí)際應(yīng)用中，給定狀態(tài)轉(zhuǎn)移矩陣、混淆矩陣和初始狀態(tài)概率可確定隱馬爾可夫模型，可按照如下過(guò)程產(chǎn)生相應(yīng)的觀測(cè)序列：

（1）設(shè)置時(shí)間步為1，根據(jù)初始狀態(tài)的概率選擇初始狀態(tài)值；

（2）根據(jù)當(dāng)前狀態(tài)值和觀測(cè)概率選擇觀測(cè)變量；

（3）根據(jù)當(dāng)前狀態(tài)值和狀態(tài)轉(zhuǎn)移矩陣選擇下一時(shí)間步的狀態(tài)值；

（4）直至滿足確定的時(shí)間步要求。

常見(jiàn)的和隱馬爾可夫模型相關(guān)問(wèn)題如下：

（1）給定隱馬爾可夫模型的初始狀態(tài)、狀態(tài)轉(zhuǎn)移矩陣和混淆矩陣，計(jì)算產(chǎn)生給定觀測(cè)序列的概率。即根據(jù)已有的觀測(cè)序列推測(cè)當(dāng)前時(shí)刻最可能出現(xiàn)的觀測(cè)值。

（2）給定隱馬爾可夫模型的初始狀態(tài)、狀態(tài)轉(zhuǎn)移矩陣、混淆矩陣和觀測(cè)序列。計(jì)算與之對(duì)應(yīng)的隱藏狀態(tài)序列，即根據(jù)觀測(cè)信號(hào)推斷最大可能的隱藏狀態(tài)序列。

（3）根據(jù)給定的觀測(cè)序列，計(jì)算隱馬爾可夫模型的參數(shù)，且使得出現(xiàn)給定序列的概率最大，常用于通過(guò)學(xué)習(xí)樣本來(lái)得到模型參數(shù)，以便更好地描述觀測(cè)數(shù)據(jù)。

2.2.2　貝葉斯網(wǎng)絡(luò)

貝葉斯網(wǎng)絡(luò)也稱信念網(wǎng)絡(luò)或有向無(wú)環(huán)圖模型，是一種概率圖模型，于1985年由Judea Pearl首先提出，是一種模擬人類推理過(guò)程中因果關(guān)系的不確定性處理模型。有向無(wú)環(huán)圖可表示其網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)，圖中的節(jié)點(diǎn)表示隨機(jī)變量，可以表示可觀察到的變量、隱變量、未知參數(shù)等。有因果關(guān)系或非條件獨(dú)立的變量和命題在有向無(wú)環(huán)圖中通過(guò)箭頭對(duì)表示“因”和“果”的兩個(gè)節(jié)點(diǎn)進(jìn)行連接，連接在一起的兩個(gè)節(jié)點(diǎn)會(huì)產(chǎn)生一個(gè)條件概率值。例如，假定節(jié)點(diǎn)A直接影響到節(jié)點(diǎn)B，通過(guò)A→B來(lái)表示，可用從節(jié)點(diǎn)A指向節(jié)點(diǎn)B的箭頭來(lái)構(gòu)建A到B的有向弧（A，B），權(quán)值可用條件概率來(lái)表示，用圓圈表示隨機(jī)變量，用箭頭表示變量之間的條件依賴，如圖2.4所示。

將某系統(tǒng)中涉及的全部隨機(jī)變量，根據(jù)變量之間是否條件獨(dú)立表示在一個(gè)有向圖中，用來(lái)表示隨機(jī)變量之間的條件依賴，即構(gòu)成了貝葉斯網(wǎng)絡(luò)。G=（I，E）表示有向無(wú)環(huán)圖，式中I表示有向無(wú)環(huán)圖中所有的節(jié)點(diǎn)的集合，E表示有向連接弧的集合，假設(shè)X=（Xj），j∈I表示圖2.4中的某節(jié)點(diǎn)j表示的隨機(jī)變量，則節(jié)點(diǎn)X的聯(lián)合概率可以表示為

X即為相對(duì)于有向無(wú)環(huán)圖G的貝葉斯網(wǎng)絡(luò)。

概率的先驗(yàn)獨(dú)立性和條件獨(dú)立性如下：

先驗(yàn)獨(dú)立性：P（X，Y）=P（X）×P（Y）

條件獨(dú)立性：

采用多個(gè)條件概率的乘積來(lái)計(jì)算聯(lián)合概率的方法叫作鏈?zhǔn)椒▌t，每一個(gè)貝葉斯網(wǎng)絡(luò)都對(duì)應(yīng)了一套鏈?zhǔn)椒▌t來(lái)計(jì)算聯(lián)合概率，同時(shí)每一個(gè)鏈?zhǔn)椒▌t在忽略弱依賴關(guān)系時(shí)，對(duì)于任意隨機(jī)變量，聯(lián)合概率可由各局部條件概率分布相乘求得：

簡(jiǎn)單貝葉斯網(wǎng)絡(luò)如圖2.5所示：

因?yàn)?span id="k2jhmkl" class="emphasis_italic">A導(dǎo)致B，A和B導(dǎo)致C，所以有下式：

貝葉斯網(wǎng)絡(luò)共有三種結(jié)構(gòu)形式，分別為head-to-head、tail-to-tail和head-to-tail。

貝葉斯網(wǎng)絡(luò)如圖2.6所示。

從圖2.6可以比較直觀地看出：

x1，x2，…，x7的聯(lián)合分布為：

（1）；

（2）x1、x2獨(dú)立；

（3）x4和x5在x6給定的條件下獨(dú)立。

拓?fù)浣Y(jié)構(gòu)1：匯連結(jié)構(gòu)（head-to-head）。匯連結(jié)構(gòu)形式如圖2.7所示：

則有成立，化簡(jiǎn)可得：

可得：

在x7未知的條件下，x1、x2是被阻斷的，稱為head-to-head條件獨(dú)立。

拓?fù)浣Y(jié)構(gòu)2：分連結(jié)構(gòu)（tail-to-tail）。分連結(jié)構(gòu)形式如圖2.8所示：

考慮C未知和C已知這兩種情況：

在x6未知時(shí)，有，此時(shí)無(wú)法得到p（x4，x5）=p（x4）p（x5），即x6未知時(shí)，x4、x5不獨(dú)立。

在x6已知時(shí)，有，將p（x4，x5，x6）=代入式中，可以得到：即x6已知時(shí)，x4、x5獨(dú)立。因此，當(dāng)x6已知時(shí)，x4、x5是被阻斷的，稱為tail-to-tail條件獨(dú)立，即x4和x5在x6給定的條件下獨(dú)立。

拓?fù)浣Y(jié)構(gòu)3：直連結(jié)構(gòu)（head-to-tail）。直連結(jié)構(gòu)形式如圖2.9所示：

當(dāng)x6未知時(shí)，有，無(wú)法推出p（x1，x4）=p（x1）p（x4），即x6未知時(shí)，x1、x4不獨(dú)立。

當(dāng)x6已知時(shí)，有，根據(jù)p（x1，，可得：

即在x6已知時(shí)，x1、x4被阻斷（blocked），稱為head-to-tail條件獨(dú)立。

顯然，head-to-tail其實(shí)就是一個(gè)鏈?zhǔn)骄W(wǎng)絡(luò)，在x6給定的條件下，x4的分布和x1條件獨(dú)立，也就是x4的分布狀態(tài)只和x6有關(guān)，和其他變量條件獨(dú)立。即當(dāng)前狀態(tài)只跟上一狀態(tài)有關(guān)，這一隨機(jī)過(guò)程叫作馬爾可夫鏈。

例2.1：煤炭是工業(yè)的食糧，是人類的生產(chǎn)生活必不可缺的能量來(lái)源之一，煤炭的供應(yīng)也關(guān)系到各國(guó)的工業(yè)乃至整個(gè)社會(huì)發(fā)展的穩(wěn)定。造成煤炭?jī)r(jià)格波動(dòng)的因素較為復(fù)雜，同時(shí)煤炭?jī)r(jià)格的波動(dòng)也會(huì)影響其他相關(guān)行業(yè)的健康發(fā)展。當(dāng)煤炭的價(jià)格上漲時(shí)，有80%的可能性會(huì)引起當(dāng)?shù)仉妰r(jià)的上漲，并且有40%的可能性會(huì)帶來(lái)失業(yè)率的上升；如果煤炭?jī)r(jià)格穩(wěn)定，引起電價(jià)上漲和失業(yè)率提高的可能性只有10%。若煤炭的產(chǎn)量持續(xù)增加且不出現(xiàn)長(zhǎng)期的高溫天氣，煤炭有70%的可能性保持價(jià)格穩(wěn)定；如果產(chǎn)量下降或出現(xiàn)持續(xù)高溫天氣，煤炭有60%的可能性保持價(jià)格穩(wěn)定。當(dāng)產(chǎn)量下降的同時(shí)又伴隨持續(xù)高溫，則煤炭只有30%的可能性保持價(jià)格穩(wěn)定。根據(jù)統(tǒng)計(jì)數(shù)據(jù)得知，煤炭產(chǎn)量持續(xù)增加的可能性為50%，同時(shí)有20%的可能性出現(xiàn)持續(xù)高溫天氣。

上述問(wèn)題可以通過(guò)圖2.10有向無(wú)環(huán)圖表示，橢圓節(jié)點(diǎn)表示變量，有向邊表示變量之間的依賴關(guān)系，起始節(jié)點(diǎn)表示父節(jié)點(diǎn)，終止節(jié)點(diǎn)為子節(jié)點(diǎn)，比如圖2.10中Y節(jié)點(diǎn)是C節(jié)點(diǎn)的父節(jié)點(diǎn)，C節(jié)點(diǎn)為Y節(jié)點(diǎn)的子節(jié)點(diǎn)。每個(gè)節(jié)點(diǎn)的概率值以條件概率表的形式給出。條件概率表是在父節(jié)點(diǎn)約束下的子節(jié)點(diǎn)條件分布，變量以二值的形式取值，比如Y節(jié)點(diǎn)取值T時(shí)表示煤炭產(chǎn)量穩(wěn)定，取值為F時(shí)表示產(chǎn)量不穩(wěn)定，則該問(wèn)題各變量的聯(lián)合概率分布為P（Y，H，C，E，R）。根據(jù)條件概率公式和變量的獨(dú)立性，可以得到如下的等量變換：

2.2.3　貝葉斯網(wǎng)絡(luò)推理

貝葉斯網(wǎng)絡(luò)推理算法一般分為精確推理算法和近似推理算法。其中精確推理算法希望能計(jì)算出目標(biāo)變量的邊際分布或條件分布的精確值，但此類算法的計(jì)算復(fù)雜度隨著極大團(tuán)規(guī)模的增長(zhǎng)呈指數(shù)增長(zhǎng)，不適用于規(guī)模較大的貝葉斯網(wǎng)絡(luò)，因此當(dāng)貝葉斯網(wǎng)絡(luò)的規(guī)模較大時(shí)，一般采用近似推理，以便在較低時(shí)間復(fù)雜度下獲得問(wèn)題的近似解。常見(jiàn)的貝葉斯精確推理算法有：多樹(shù)傳播（Polytree Propagation）推理算法；團(tuán)樹(shù)傳播（Clique Tree Propagation）推理算法，比如聯(lián)結(jié)樹(shù)（Junction Tree Propagation）推理算法；基于組合優(yōu)化的求解算法，如符號(hào)推理（Symbolic ProbabilisticInference）和桶消元推理（Bucket Elemination Inference）算法。近似推理算法主要有：基于搜索的（Search-Based）方法、Monte Carlo算法等。

2.2.3.1　精確推理

貝葉斯網(wǎng)絡(luò)精確推理的常用方法有變量消元法、枚舉推理法等，本質(zhì)上是通過(guò)動(dòng)態(tài)規(guī)劃算法來(lái)計(jì)算目標(biāo)概率值，通過(guò)給定的證據(jù)變量計(jì)算查詢變量的后驗(yàn)概率分布。變量消元法是一種基于變量的條件獨(dú)立性的精確推理算法，通過(guò)計(jì)算邊際概率消去非計(jì)算目標(biāo)的概率值，如圖2.11所示。

變量消元法就是對(duì)聯(lián)合概率求和消去其他變量得到邊緣分布，如圖2.11首先把z消元，得到僅含x和y的表，然后再對(duì)y進(jìn)行求和得到只含x的概率表。

以圖2.12中的貝葉斯網(wǎng)為例，計(jì)算P（D2），則有：

為了降低推理計(jì)算的復(fù)雜度，可以分解聯(lián)合分布的因子。與變量A有關(guān)的因子是P（A）和，與變量B有關(guān)的因子是和，與變量C有關(guān)的因子是和，因此，有：

消元變量A只涉及自身和變量B，消元變量B只涉及自身和變量C，消元變量D1只涉及自身和變量C，消元變量C只涉及自身和變量D2，因此可以通過(guò)局部運(yùn)算有效降低推理的復(fù)雜度，可以如下描述消元法的算法。

（1）設(shè)貝葉斯網(wǎng)中概率分布的集合為Ω，有：

（2）將集合Ω中含變量A的因子P（A）和合并求和得到新因子：

并將新因子f1（B）加入集合Ω，同時(shí)消去因子P（A）和，達(dá)到消元變量A的目的，得到更新的集合Ω：

（3）將集合Ω中含變量B的因子f1（B）和合并求和得到新因子：

并將新因子f2（C）加入集合Ω，同時(shí)消去因子f1（B）和，達(dá)到消元變量B的目的，得到更新的集合Ω：

（4）將集合Ω中含變量D1的因子計(jì)算得到新因子：

并將新因子f（C）加入集合Ω，同時(shí)消去因子f1（B）和，達(dá)到消元變量D1的目的，得到更新的集合Ω：

（5）將集合Ω中含變量C的因子f2（C）、f3（C）和合并求和得到新因子：

達(dá)到消元變量C的目的。

（6）返回f4（D2），即問(wèn)題的目標(biāo)P（D2），如圖2.12所示。

變量消元法計(jì)算的復(fù)雜度和消元的順序有關(guān)系，而最優(yōu)的消元順序是NP-hard問(wèn)題，所以實(shí)際計(jì)算中一般采取啟發(fā)式的規(guī)則進(jìn)行計(jì)算，常用的算法有最大勢(shì)搜索和最小缺邊搜索。

圖2.13是包含8個(gè)節(jié)點(diǎn)的無(wú)向圖，最大勢(shì)搜索是對(duì)圖中所有節(jié)點(diǎn)按照如下規(guī)則進(jìn)行編號(hào)，在t步中選擇具有最多已編號(hào)相鄰節(jié)點(diǎn)的未編號(hào)節(jié)點(diǎn)，若有多個(gè)則選擇其中一個(gè)并將編號(hào)定為8-t+1。當(dāng)所有節(jié)點(diǎn)被編號(hào)后，按照從小到大的順序排序節(jié)點(diǎn)，即可得到變量消元順序。

如圖2.13，任意選擇節(jié)點(diǎn)H，并編號(hào)為8。下一可選節(jié)點(diǎn)為有1個(gè)已編號(hào)鄰接節(jié)點(diǎn)H的E、F或G，任意選擇節(jié)點(diǎn)F，并編號(hào)為7。下一節(jié)點(diǎn)則只能選擇有2個(gè)已編號(hào)鄰接節(jié)點(diǎn)H和F的E，編號(hào)為6。下一節(jié)點(diǎn)選擇有2個(gè)已編號(hào)鄰接節(jié)點(diǎn)E和H的G，編號(hào)為5。下一可選節(jié)點(diǎn)有1個(gè)已編號(hào)鄰接節(jié)點(diǎn)的C、D或A，任意選擇節(jié)點(diǎn)A，并編號(hào)為4。下一節(jié)點(diǎn)選擇有2個(gè)已編號(hào)鄰接節(jié)點(diǎn)A和E的D，編號(hào)為3。下一節(jié)點(diǎn)選擇有2個(gè)已編號(hào)鄰接節(jié)點(diǎn)D和E的C，編號(hào)為2。最后將節(jié)點(diǎn)B編號(hào)為1，得到最終的消元順序?yàn)椋?span id="vx7rgwc" class="emphasis_italic">B，C，D，A，G，E，F，H＞。

變量消元算法的缺點(diǎn)主要表現(xiàn)在需要計(jì)算多個(gè)邊緣分布，多次執(zhí)行變量消元算法，中間會(huì)有很多結(jié)果可以復(fù)用，但是變量消元算法需計(jì)算多次，導(dǎo)致效率低下。

2.2.3.2　近似推理

精確推理需要大量的計(jì)算開(kāi)銷，在實(shí)際應(yīng)用中面對(duì)較大規(guī)模的問(wèn)題時(shí)近似推理更為常用。近似推理的核心思想是通過(guò)隨機(jī)采樣得到相應(yīng)的樣本集，進(jìn)而通過(guò)樣本的頻率來(lái)逼近查詢概率。蒙特卡洛方法也被稱為統(tǒng)計(jì)模擬方法、統(tǒng)計(jì)試驗(yàn)法、隨機(jī)抽樣技術(shù)等，該方法于20世紀(jì)40年代由美國(guó)在第二次世界大戰(zhàn)中研制原子彈的“曼哈頓計(jì)劃”的成員S.M.烏拉姆和J.馮·諾伊曼首先提出，馮·諾伊曼用賭城摩納哥的Monte Carlo來(lái)命名。廣泛應(yīng)用于人工智能、金融工程學(xué)、宏觀經(jīng)濟(jì)學(xué)、生物醫(yī)學(xué)、計(jì)算物理學(xué)，諸如粒子輸運(yùn)計(jì)算、量子熱力學(xué)計(jì)算、空氣動(dòng)力學(xué)計(jì)算、核工程等領(lǐng)域。常見(jiàn)的采樣方法有直接采樣方法、接受—拒絕采樣方法、重要性采樣方法、MCMC采樣方法、Metropolis-Hastings采樣方法等。

（1）直接采樣方法。

蒙特卡洛方法是通過(guò)隨機(jī)模擬的方式求解不太容易求解的和或者積分問(wèn)題，使用采樣+平均核心思想估計(jì)無(wú)法計(jì)算的值。

直接采樣的方法是根據(jù)概率分布進(jìn)行采樣。對(duì)一個(gè)已知概率密度函數(shù)與累積概率密度函數(shù)的概率分布，我們可以直接從累積分布函數(shù)進(jìn)行采樣。

在例2.1中，考慮如下推理問(wèn)題：假設(shè)觀察到失業(yè)率提高，計(jì)算該結(jié)果由于煤炭產(chǎn)量下降引起的概率是多少。針對(duì)此問(wèn)題，首先可以根據(jù)煤炭產(chǎn)量和高溫天氣的先驗(yàn)概率生成一個(gè)樣本，并記錄樣本的取值；其次以樣本的取值為條件，按照煤炭?jī)r(jià)格上漲的條件概率生成樣本，并記錄樣本的取值；最后根據(jù)獨(dú)立向分別計(jì)算在煤炭?jī)r(jià)格為條件時(shí)生成電價(jià)和失業(yè)率的樣本。產(chǎn)生大量樣本后，統(tǒng)計(jì)取值＜*，r，*，y，*＞和＜*，r，*，*，*＞的樣本個(gè)數(shù)，分別記為N（*，r，*，y，*）和N（*，r，*，*，*），根據(jù)概率的定義，當(dāng)采集的樣本量較大時(shí)，則N（*，r，*，y，*）/N（*，r，*，*，*）即可作為近似結(jié)果輸出，當(dāng)樣本總量趨于無(wú)窮時(shí)，結(jié)果即為問(wèn)題的概率。在實(shí)際應(yīng)用中，由于很多分布無(wú)法寫出概率密度函數(shù)與累積分布函數(shù)，因此直接采樣方法的適用范圍較窄。

例2.2：應(yīng)用蒙特卡洛方法求圓周率，如圖2.14所示。圓的半徑為r=1，正方形的邊長(zhǎng)為1，根據(jù)圓面積的計(jì)算公式可知圓的面積為πr2=π，同時(shí)正方形內(nèi)部的相切圓的面積為整個(gè)圓面積的1/4，也就是。正方形的面積為1。然后向正方形中隨機(jī)打點(diǎn)，則點(diǎn)就會(huì)以一定的概率落在1/4圓中，而點(diǎn)落在1/4圓中的概率為：

圓的面積/正方形面積=，

然后即可推出圓周率的計(jì)算公式為：

（2）接受—拒絕采樣方法。

接受—拒絕采樣方法用于累積分布函數(shù)未知的問(wèn)題。如圖2.15所示，p（z）是希望采樣的分布，很多實(shí)際問(wèn)題中，p（z）因太復(fù)雜在程序中沒(méi)法直接采樣，需要借助其他的手段來(lái)采樣，因此設(shè)置提議的程序可抽樣分布q（z），比如高斯分布，且kq（z）＞p（z）。在kq（z）中按照直接采樣的方法采樣，然后判斷這個(gè)樣本的取值是否符合給定觀測(cè)，不符合給定觀測(cè)的樣本予以拒絕，符合給定觀測(cè)的樣本予以接受，最終得到符合p（z）的N個(gè)樣本。樣本量越大，得到的查詢概率越接近真實(shí)值，當(dāng)樣本量區(qū)域無(wú)窮大時(shí)得到的查詢概率和真實(shí)值嚴(yán)格相等。接受—拒絕采樣方法因?yàn)樵诓蓸舆^(guò)程中要丟棄大量不符合觀測(cè)值的樣本，所以計(jì)算效率不高。

（3）重要性采樣方法。

接受—拒絕采樣解決了累積分布函數(shù)未知的采樣問(wèn)題，但采樣的效果取決于提議分布的選擇，如果提議分布選擇得不好，則會(huì)出現(xiàn)計(jì)算效率低下且難以得到符合觀測(cè)值的樣本。與直接采樣與接受—拒絕采樣將每個(gè)樣本默認(rèn)相同的權(quán)重不同，重要性采樣給予每個(gè)樣本不同的權(quán)重，使用加權(quán)平均的方法來(lái)計(jì)算期望。

通過(guò)從提議分布q（z）中采樣樣本z1，z2，…，zn，每個(gè)樣本的權(quán)重是，通過(guò)加權(quán)平均的方式計(jì)算期望，如：

（4）MCMC抽樣方法。

馬爾可夫鏈?zhǔn)歉怕收摵蛿?shù)理統(tǒng)計(jì)中具有馬爾可夫性質(zhì)且存在于離散的指數(shù)集和狀態(tài)空間內(nèi)的隨機(jī)過(guò)程，也就是通過(guò)時(shí)間串聯(lián)的一系列隨機(jī)變量，可通過(guò)轉(zhuǎn)移矩陣和轉(zhuǎn)移圖定義。馬爾可夫鏈可被應(yīng)用于蒙特卡洛方法中，形成馬爾可夫鏈蒙特卡洛采樣方法，即MCMC方法。

舉例說(shuō)明馬爾可夫鏈。按照十年河?xùn)|十年河西的說(shuō)法，窮人和富人之間會(huì)在一定概率下相互轉(zhuǎn)變。為便于計(jì)算，假定窮人人口和富人人口每年會(huì)發(fā)生一次轉(zhuǎn)變，且轉(zhuǎn)變概率是固定的。假定每年富人轉(zhuǎn)變?yōu)楦F人的概率是4%，窮人轉(zhuǎn)變?yōu)楦蝗说母怕适?%。如果某一年，某區(qū)域的富人和窮人的人口分別是5000和35000，那么該區(qū)域下一年窮人和富人的分布如何呢？

首先構(gòu)造轉(zhuǎn)移矩陣：

則富人和窮人的分布為：

因此，每個(gè)人作為一個(gè)個(gè)體，其個(gè)人的窮富狀態(tài)其實(shí)根據(jù)個(gè)人的奮斗情況在轉(zhuǎn)變：比如t=1時(shí)，是窮人；t=2時(shí)，經(jīng)過(guò)奮斗則會(huì)變成富人。馬爾可夫鏈就是生成這樣一段狀態(tài)序列的隨機(jī)過(guò)程，其中城市和農(nóng)村互相流動(dòng)的矩陣，叫作遷移矩陣。馬爾可夫鏈的這個(gè)隨機(jī)過(guò)程滿足馬爾可夫性質(zhì)，也就是某一個(gè)狀態(tài)的值，只跟前一個(gè)狀態(tài)相關(guān)，即一階齊次馬爾可夫鏈。用公式表示就是：

在遷移矩陣中，每一個(gè)元素對(duì)應(yīng)一個(gè)條件概率值，經(jīng)過(guò)一定階段的迭代之后，最終窮人和富人的人數(shù)相對(duì)固定了，富人是26667，窮人是13333。對(duì)于該區(qū)域而言，只要總?cè)丝谑?0000，無(wú)論初始的窮人和富人比例如何，最終都會(huì)收斂到富人是26667，窮人是13333。簡(jiǎn)單分析原因，窮人經(jīng)過(guò)奮斗會(huì)變得富有，而富人如果墮落則會(huì)變得窮困。隨著時(shí)間的推移，窮人逐漸減少，而富人逐漸增加，但最終窮人變富和富人變窮的數(shù)量是相等的，即兩類人群的流入和流出是相等的，達(dá)到一種相對(duì)平穩(wěn)的分布，這個(gè)分布就是馬爾可夫鏈的平穩(wěn)分布。

馬爾可夫鏈?zhǔn)諗坑谄椒€(wěn)分布π，π是方程πP=π唯一非負(fù)解。把這個(gè)解用向量方式表示，即：

且，。

因此，對(duì)于一個(gè)非周期的馬爾可夫鏈轉(zhuǎn)移矩陣P，其任何兩個(gè)狀態(tài)是連通的，存在，并且與i獨(dú)立，記為，則有：

π是方程πP=π唯一非負(fù)解：

且，π=［π（1），π（2），…，π（j），…］，

則π成為馬爾可夫鏈的平穩(wěn)分布。

MCMC方法的核心是對(duì)于任意給定的目標(biāo)分布，找到以該目標(biāo)分布為穩(wěn)態(tài)分布的馬爾可夫鏈，并基于馬爾可夫鏈采樣方法對(duì)其近似采樣。

（5）Metropolis-Hastings采樣方法。

若要基于給定的概率分布高效生成對(duì)應(yīng)的樣本，由于馬爾可夫鏈可以收斂到平穩(wěn)分布，最直觀的需要是構(gòu)造一個(gè)轉(zhuǎn)移矩陣為P的馬爾可夫鏈，使得該馬爾可夫鏈的平穩(wěn)分布是π（x），則從任何一個(gè)初始狀態(tài)均可以得到相應(yīng)的轉(zhuǎn)移序列，收集馬爾可夫鏈在收斂之后的樣本。

所以，基于馬爾可夫鏈采樣的關(guān)鍵問(wèn)題是構(gòu)造轉(zhuǎn)移矩陣，使得轉(zhuǎn)移矩陣滿足一定的條件，使平穩(wěn)分布恰好是目標(biāo)分布P（x）。

轉(zhuǎn)移矩陣需要滿足的這個(gè)條件就是細(xì)致平穩(wěn)條件，細(xì)致平穩(wěn)條件是指對(duì)于給定的馬爾可夫鏈的轉(zhuǎn)移矩陣P和分布π，以及馬爾可夫鏈中的任意兩個(gè)狀態(tài)x和x*，如果有：

則分布π即為該馬爾可夫鏈的平穩(wěn)分布。

簡(jiǎn)單分析細(xì)致平穩(wěn)條件，對(duì)于條件中等式兩側(cè)關(guān)于狀態(tài)x求和，則有：

則有，可得到平穩(wěn)分布的條件πP=π，即滿足細(xì)致平穩(wěn)條件的目標(biāo)分布是基于轉(zhuǎn)移矩陣P的馬爾可夫鏈平穩(wěn)分布。

Metropolis-Hastings采樣方法提供了找到滿足細(xì)致平穩(wěn)條件的轉(zhuǎn)移矩陣P的思路。

在一般情況下，對(duì)于給定的目標(biāo)概率π，任意轉(zhuǎn)移矩陣P不滿足π（x），即此時(shí)的轉(zhuǎn)移矩陣不滿足細(xì)致平穩(wěn)條件，Metropolis-Hastings方法試圖通過(guò)引入α（x，x*）和α（x*，x）因子來(lái)滿足細(xì)致平穩(wěn)條件，即：

為方便表示，將記為p（x，x*），將記為p（x*，x）。為了使上式成立，按照對(duì)稱性，只需：

于是具備了設(shè)置因子α（x，x*）和α（x*，x）的可能性，從而等式π（x）p（x，x*）α（x，x*）=π（x*）p（x*，x）α（x*，x）可以得到滿足，對(duì)其進(jìn)行如下改寫：

則有：

于是就將原轉(zhuǎn)移矩陣改造成了轉(zhuǎn)移矩陣Q，且Q滿足細(xì)致平穩(wěn)條件，引入的因子被稱為接受率，物理意義是以此接受率為概率接受轉(zhuǎn)移。以接受率α（x，x*）為例，當(dāng)從狀態(tài)x轉(zhuǎn)移到狀態(tài)x*時(shí)，以α（x，x*）為概率接受這個(gè)轉(zhuǎn)移，則新的馬爾可夫鏈的轉(zhuǎn)移概率由原來(lái)的p（x，x*）變?yōu)?span id="bkm7app" class="emphasis_italic">p（x，x*）α（x，x*）。如圖2.16（a）所示。

狀態(tài)2以0.4的概率向狀態(tài)1轉(zhuǎn)變，以0.5的概率向狀態(tài)3轉(zhuǎn)變，以0.1的概率保持在狀態(tài)2。疊加上接受率α（2，1）=0.6和α（2，3）=0.8之后，如圖2.16（b）所示。

疊加上接受率之后圖示的轉(zhuǎn)移概率發(fā)生了改變，狀態(tài)2到狀態(tài)1的轉(zhuǎn)移概率原為0.4，因?yàn)榀B加了接受率α（2，1）=0.6，則狀態(tài)2到狀態(tài)1的轉(zhuǎn)移概率變?yōu)?.4×0.6=0.24，有0.4×0.4=0.16的概率仍保留在狀態(tài)2；狀態(tài)2到狀態(tài)3的轉(zhuǎn)移概率原為0.5，因?yàn)榀B加了接受率α（2，3）=0.8，則狀態(tài)2到狀態(tài)3的轉(zhuǎn)移概率變?yōu)?.5×0.8=0.4，有0.5×0.2=0.1的概率仍保留在狀態(tài)2；而保留狀態(tài)2的概率則由原來(lái)的概率0.1增大到0.36（0.1+0.16+0.1=0.36）。

在馬爾可夫鏈蒙特卡洛方法中的Metropolis-Hastings方法中進(jìn)一步明確了接受概率α的表達(dá)式：

以圖2.16為例，假定接受率α（2，1）=0.1和α（2，3）=0.2，在實(shí)際采樣過(guò)程中效率會(huì)很低，這是因?yàn)樵诘仁?span id="b3zvhez" class="emphasis_italic">π（x）p（x，x*）×0.1=π（x*）p（x*，x）×0.2中等號(hào)兩端的接受率過(guò)小，此時(shí)如果對(duì)等號(hào)兩端的接受率等比例放大，比如：

顯然仍滿足細(xì)致平穩(wěn)條件，而采樣的效率可以得到大幅提高。