在（上）篇中，我們討論了什么是不確定度，為什么需要關(guān)注不確定度建模，以及不確定度可以怎么用。也從最大似然估計(jì)（MLE）到最大后驗(yàn)概率（MAP），講到了貝葉斯推斷(Bayesian Inference)。而我們希望用來(lái)建模不確定的目標(biāo)模型是貝葉斯神經(jīng)網(wǎng)絡(luò)（BNN），它是一種用神經(jīng)網(wǎng)絡(luò)來(lái)建模似然率，然后進(jìn)行貝葉斯推斷的方法。（中）篇會(huì)介紹如何用蒙特卡洛采樣的方法來(lái)進(jìn)行貝葉斯推斷。主要是概念圖中的這一部分：

回顧一下，貝葉斯推斷通常遵循以下的四步：

第一步：假設(shè)先驗(yàn)??，對(duì)似然率建模?
第二步：計(jì)算后驗(yàn)?
展開(kāi)后有?
第三步：計(jì)算預(yù)測(cè)值的分布?
第四步：計(jì)算y的期望值??；
如果需要，計(jì)算y的方差??作為預(yù)估值的不確定度

通常我們很容易的可以把后驗(yàn)分布的表達(dá)式??寫(xiě)出來(lái)，但是到第三步，是一個(gè)積分，這里就比較麻煩了。如果我們選擇的似然率的模型和先驗(yàn)?zāi)Ｐ筒皇莿偤糜行┨厥庑再|(zhì)可以利用，很難直接把積分求出來(lái)。另外第四步里也需要用積分來(lái)求期望值（或者方差）。

如果不追求精確的解析解，而可以接受近似解的話，有兩個(gè)方法可以解決這個(gè)問(wèn)題：一個(gè)是蒙特卡洛采樣，一個(gè)是變分推斷。

一、蒙特卡洛采樣

我們先來(lái)看蒙特卡洛采樣的方法。拉斯維加斯，澳門(mén)，都是世界著名賭城，還有另外一座城市，和它們并稱“世界三大賭城”。這個(gè)城市的名字就叫做蒙特卡洛，這個(gè)算法的名字也就起源于這個(gè)城市。

如果我們是蒙特卡洛賭城的老板，我們需要設(shè)計(jì)一臺(tái)老虎機(jī)，怎么做可以知道這臺(tái)老虎機(jī)是不是穩(wěn)賺不賠的？一個(gè)最簡(jiǎn)單的方法就是找一批測(cè)試人員輪流模擬玩家在這臺(tái)機(jī)器上試玩1000次，看看最終機(jī)器是不是賺錢(qián)了，賺了多少錢(qián)，然后再除以1000，得到的就是玩家平均每玩一次賭場(chǎng)賺的錢(qián)。這個(gè)叫做老虎機(jī)的ev（expected value）。賭場(chǎng)的機(jī)器，對(duì)于玩家來(lái)說(shuō)都是負(fù)ev的。所以長(zhǎng)期玩下去，賭場(chǎng)肯定是賺的。

這種求單次賭場(chǎng)平均盈利的方法，就是蒙特卡洛采樣法。更通用地來(lái)說(shuō)，蒙特卡洛采樣法是指用從一個(gè)分布中采樣很多樣本，然后用這些樣本來(lái)代替這個(gè)分布進(jìn)行相關(guān)計(jì)算（例如求期望值）的方法。比如說(shuō)我們要下面這個(gè)積分

回顧一下期望值計(jì)算公式，上式其實(shí)就是當(dāng)x的分布為??時(shí)，的期望值。

如果我們不能直接把積分的解析解算出來(lái)，就可以用蒙特卡洛采樣的方法，從?中采樣很多樣本??（類比賭場(chǎng)找來(lái)的測(cè)試人員）然后把這些樣本帶入??（類比測(cè)試人員去玩一次老虎機(jī)），得到（類比測(cè)試人員玩老虎機(jī)的盈虧）再把他們的均值算出來(lái)就得到了這個(gè)積分的一個(gè)近似解（類比老虎機(jī)的ev）。采樣的樣本越多，得到的解越接近解析解（大數(shù)定律）

我們先回到貝葉斯推斷里關(guān)鍵的第三步和第四步，把第三步的分布帶入到第四步的期望值公式里，則有??。因?yàn)閥與??無(wú)關(guān)，可以放到第二個(gè)積分內(nèi)部，則有。再調(diào)整下積分的順序，并把與y無(wú)關(guān)的??挪到外面，則有?

上式中，我們會(huì)發(fā)現(xiàn)括號(hào)內(nèi)的積分，其實(shí)就是?的均值。在大多數(shù)時(shí)候，我們?cè)趯?duì)似然率建模的時(shí)候，會(huì)間接通過(guò)對(duì)均值建模來(lái)實(shí)現(xiàn)。例如在回歸問(wèn)題里，神經(jīng)網(wǎng)絡(luò)的輸出作為高斯分布的均值來(lái)建模似然率，即??，其中??表示神經(jīng)網(wǎng)絡(luò)的輸入和參數(shù)，??表示神經(jīng)網(wǎng)絡(luò)的輸出值，也是高斯分布的均值。?一般是個(gè)超參數(shù)。如果是分類問(wèn)題，則神經(jīng)網(wǎng)絡(luò)的輸出會(huì)作為伯努利分布里的y=1的概率，也是伯努利分布的均值。所以上面的式子就變成?

所以，利用蒙特卡洛采樣的方法，如果我們能從??分布里采樣很多??，然后分別帶入到里算出神經(jīng)網(wǎng)絡(luò)的一系列輸出值，再求均值，也就是最終我們要的??了。

那么現(xiàn)在的問(wèn)題就變成：我們?nèi)绾胃鶕?jù)一個(gè)分布的表達(dá)式（在這里主要是指），從這個(gè)分布中采樣出很多樣本來(lái)。（采樣問(wèn)題）

另外一個(gè)問(wèn)題是，我們要采樣的分布??里，有一個(gè)積分，這個(gè)積分也很不好算（歸一化問(wèn)題）。不過(guò)它是作為一個(gè)歸一化分母來(lái)使用的，如果我們能在采樣的時(shí)候，有辦法忽略這個(gè)歸一化分母就好了。大家記住這個(gè)問(wèn)題，后面我們會(huì)多次提到怎么來(lái)化解。

二、基本分布的采樣方法

對(duì)于均勻分布，我們可以用隨機(jī)數(shù)發(fā)生器很容易的采樣出來(lái)。借助隨機(jī)數(shù)發(fā)生器，我們可以實(shí)現(xiàn)另外一些簡(jiǎn)單的分布，例如伯努利分布，用簡(jiǎn)單的if else就可以實(shí)現(xiàn)；也可以實(shí)現(xiàn)稍微復(fù)雜一些的分布，比如高斯分布，可以借助Box-Muller方法來(lái)實(shí)現(xiàn)：

Box-Muller方法：選取兩個(gè)服從[0,1]上均勻分布的隨機(jī)變量U1，U2，則X，Y為均值為0，方差為1的高斯分布，且X，Y獨(dú)立。

因?yàn)椴皇潜疚闹攸c(diǎn)，不展開(kāi)證明。

但是對(duì)于更復(fù)雜的分布，采樣就變得并不容易了。有一系列采樣算法來(lái)幫助我們通過(guò)上面介紹的這些基本分布的采樣方法，來(lái)實(shí)現(xiàn)從我們想要的分布中去采樣樣本。

三、復(fù)雜分布的采樣算法

1. 拒絕采樣

假設(shè)如果我們開(kāi)了一家飯店，開(kāi)業(yè)大酬賓搞了一個(gè)促銷活動(dòng)，來(lái)吃飯的人只要扔一次骰子，如果扔6就能享受免單。那么免單概率是1/6。如果我們發(fā)現(xiàn)這樣下去要虧本，想把免單概率調(diào)整為1/12怎么辦？那么我們可以修改規(guī)則為對(duì)扔到6的人，隨機(jī)一半概率才能免單，就可以把概率變?yōu)?/12了。具體實(shí)施方式為，要求扔到6的人，再扔一次硬幣，如果是反面，則我們拒絕免單，如果是正面，我們接受免單。這樣就達(dá)到了以1/12概率采樣的效果。

拒絕采樣也是類似的思想。假設(shè)我們已經(jīng)能從分布q中采樣樣本（這個(gè)分布也叫Proposal Distribution，例如一個(gè)高斯分布），而我們的目標(biāo)是要對(duì)目標(biāo)分布 p去采樣。

在一開(kāi)始，我們先對(duì)兩個(gè)分布都做一個(gè)變換。先把p分布乘以一個(gè)歸一化分母??（是個(gè)常數(shù)），得到???；仡櫳厦嬲f(shuō)的??有一個(gè)歸一化分母不好計(jì)算的問(wèn)題，而做拒絕采樣只需要用到??，因此就不需要考慮歸一化分母計(jì)算的問(wèn)題了。然后我們?cè)賹?duì)分布q的概率密度函數(shù)乘以一個(gè)常數(shù)k，使得新的概率密度函數(shù)始終在的上方，也就是讓藍(lán)線始終在紅線上方或者和紅線重合。

然后我們從分布q中進(jìn)行采樣，假設(shè)我們采樣出的樣本是??，它對(duì)應(yīng)的放大后的q分布的概率密度值是??，它對(duì)應(yīng)的目標(biāo)分布的概率密度是??，那么拒絕概率就是??，或者說(shuō)接受概率是，這樣最后被接受的樣本，就會(huì)符合??的分布。（PS：嚴(yán)格來(lái)說(shuō)因?yàn)楸环糯罅薻倍，放大了倍，都已經(jīng)不是概率函數(shù)了，因?yàn)榉e分不為1）

拒絕采樣的一個(gè)問(wèn)題是，如果拒絕率很高，則采樣的效率很低。因此proposal distribution與目標(biāo)采樣分布越接近越好，這樣被拒絕的樣本更少。另外如果在高維空間進(jìn)行采樣，則每個(gè)緯度上都要被接受，最終樣本才能被接受。緯度災(zāi)難會(huì)讓拒絕采樣在高維空間的接受率低到不能接受。我們還需要看看有沒(méi)有更好的算法。

部分現(xiàn)代丈母娘要求的維度越來(lái)越多，從文憑到車(chē)子房子，從戶口到車(chē)牌。維度越多，候選碼農(nóng)被拒絕的概率就越大。為了讓拒絕率降低一些，碼農(nóng)們只能讓自己在這些考核維度的分布，與丈母娘要求的分布越接近越好。但是丈母娘要求的分布也會(huì)水漲船高，總是高于碼農(nóng)實(shí)現(xiàn)的分布。

2. 重要性采樣

重要性采樣其實(shí)并不是真的采樣，或者說(shuō)可以認(rèn)為是一種軟采樣。我們并沒(méi)有真的拒絕掉或接受一些樣本，而是通過(guò)把這個(gè)樣本的權(quán)重調(diào)低/調(diào)高來(lái)表達(dá)對(duì)樣本的偏好（拒絕采樣中被拒絕的樣本可以認(rèn)為樣本權(quán)重被置0了）。

現(xiàn)在我們要用從p(x)里采樣的樣本來(lái)計(jì)算f(x)的均值，而我們已經(jīng)實(shí)現(xiàn)的是從q(x)中采樣（例如q(x)是個(gè)高斯分布）。那么我們拆解下f(x)均值計(jì)算的公式如下

可以發(fā)現(xiàn)，我們只要把從q(x)中采樣出來(lái)的樣本??，帶入到f(x)中，然后用??作為權(quán)重，計(jì)算加權(quán)平均就可以了。擴(kuò)展一下，如果是希望訓(xùn)練p(x)分布下的模型，但是收集的樣本是遵循q(x)分布的時(shí)候，也可以用作為樣本權(quán)重，去訓(xùn)練模型，就一定程度上抹平分布的差異。這也是重要性采樣經(jīng)常被使用的場(chǎng)景。

重要性采樣的一個(gè)明顯的弱點(diǎn)在于，在q(x)=0的地方，因?yàn)闊o(wú)法采樣到樣本，所以即使p(x)在這個(gè)地方不為0，也不會(huì)采樣到樣本。對(duì)于q(x)很接近0但是略大于0的地方雖然稍微好點(diǎn)，但是由于被采樣的概率很低，也需要采樣很多的樣本才可能采樣到這一區(qū)間，樣本效率很低。因此，和拒絕采樣類似，重要性采樣也需要proposal distribution與目標(biāo)分布越接近越好，這樣越少的樣本就可以得到一個(gè)越準(zhǔn)確的值。也同樣會(huì)有在高緯度空間效率指數(shù)級(jí)下降的問(wèn)題。不過(guò)比拒絕采樣好的是，不需要再去尋找讓始終在上方的k值了。

樂(lè)隊(duì)的夏天中，超級(jí)樂(lè)迷有10票，相當(dāng)于每個(gè)人一票，但是權(quán)重是10; 專業(yè)樂(lè)迷，每人兩票，相當(dāng)于權(quán)重是2；大眾樂(lè)迷每人1票，相當(dāng)于權(quán)重是1。這是因?yàn)榻M委會(huì)認(rèn)為，超級(jí)樂(lè)迷和專業(yè)樂(lè)迷水平更高，他們的評(píng)審水平，更接近于作品的真實(shí)水平。

重要性采樣的應(yīng)用領(lǐng)域很廣，可以用來(lái)解在建?；蚪y(tǒng)計(jì)的分布和數(shù)據(jù)來(lái)源分布不一致的問(wèn)題。例如在強(qiáng)化學(xué)習(xí)里（例如Youtube那篇年內(nèi)最大收益的《Top-K Off-Policy Correction for a REINFORCE Recommender System》），因果推斷里（用來(lái)達(dá)到CIA條件做的Propensity Score Weighting，即PSW，其實(shí)是目標(biāo)分布為均勻分布的重要性采樣的特例）都經(jīng)常用到。又是一門(mén)算法工程師居家旅行必備技術(shù)。

解決了采樣的問(wèn)題，我們來(lái)看看怎么在重要性采樣里解決的歸一化分母不好計(jì)算的問(wèn)題。我們還是用表示沒(méi)有歸一化后的目標(biāo)“概率密度”函數(shù)，用??表示沒(méi)有歸一化后的proposal distribution的“概率密度”函數(shù)。把這兩個(gè)歸一化分母提到前面，則有

這里就把計(jì)算所需要的分布從p(x)變成了沒(méi)有歸一化的??，不過(guò)我們還需要求出??才行。不過(guò)這個(gè)值并不難求。

上式的倒數(shù)，就是要求的。即對(duì)于采樣樣本??，計(jì)算??的均值，再求一個(gè)倒數(shù)就可以了。至此，歸一化問(wèn)題也被解決了。

3. MCMC（馬爾可夫鏈蒙特卡洛采樣）

接下來(lái)要介紹的MCMC的方法，在一些特殊的情況下，可以緩解或解決拒絕采樣和重要性采樣遇到的高維空間采樣效率低的問(wèn)題，并且同樣可以規(guī)避的歸一化因子難計(jì)算的問(wèn)題。

但是我們需要一點(diǎn)耐心，先從馬爾可夫鏈慢慢聊起。一階馬爾可夫鏈?zhǔn)侵敢粋€(gè)序列的隨機(jī)變量??，他們滿足下面的式子

也就是說(shuō)，??只和序列里上一個(gè)隨機(jī)變量??有關(guān)系，而和之前的其他隨機(jī)變量（例如??)都沒(méi)有關(guān)系。

馬爾可夫鏈蒙特卡洛采樣的基本想法是構(gòu)造一個(gè)特殊的馬爾可夫鏈，使得從這個(gè)序列里依次生成出來(lái)的值的集合，會(huì)遵循我們想要采樣的分布。接下來(lái)，我們看看如何構(gòu)造這種特殊的馬爾可夫鏈。在此之前，還有幾個(gè)概念要先了解下。

（a) 齊次馬爾科夫鏈

如果對(duì)于所有的m取值，轉(zhuǎn)移概率??對(duì)于任意的取值和的取值，都相同，則我們說(shuō)這個(gè)馬爾可夫鏈?zhǔn)驱R次馬爾可夫鏈。也就是說(shuō)對(duì)于任意??，??(注意這里用加括號(hào)表示馬爾可夫鏈中的隨機(jī)變量，不加表示普通的隨機(jī)變量)，齊次馬爾可夫鏈都有?

可以理解為在根據(jù)每個(gè)隨機(jī)變量的值產(chǎn)生序列里的下一個(gè)隨機(jī)變量的值時(shí)，遵循著同樣的一個(gè)轉(zhuǎn)移概率矩陣，和在序列里的位置無(wú)關(guān)。

如果我們把轉(zhuǎn)移概率表示為??，齊次馬爾可夫鏈的轉(zhuǎn)移概率因?yàn)閷?duì)所有m都相同，也就是和m的取值無(wú)關(guān)，可以計(jì)為?，即?

（b) 概率矩陣，平穩(wěn)分布與細(xì)致平穩(wěn)條件

如果我們把對(duì)應(yīng)的轉(zhuǎn)移概率放到矩陣??里，使得矩陣?yán)锩總€(gè)元素，則矩陣?就是齊次馬爾可夫鏈的轉(zhuǎn)移概率矩陣。其中表示的是從狀態(tài)??轉(zhuǎn)移到狀態(tài)??的概率。

還是剛才的例子，假設(shè)隨機(jī)變量的取值范圍是離散的集合??，我們把當(dāng)前隨機(jī)變量x等于不同取值的概率寫(xiě)到一個(gè)長(zhǎng)度為n向量里。

對(duì)于一個(gè)以向量表示的分布乘以一個(gè)矩陣，即?，物理意義就是該分布經(jīng)過(guò)一次由矩陣T定義的轉(zhuǎn)移之后，產(chǎn)生的序列里下一個(gè)隨機(jī)變量的概率分布。如果有

也就說(shuō)下一個(gè)隨機(jī)變量的概率分布和當(dāng)前一致，沒(méi)有發(fā)生變化。我們說(shuō)這個(gè)分布是關(guān)于這個(gè)馬爾可夫鏈的平穩(wěn)分布。回憶一下我們的目標(biāo)，是要讓從馬爾可夫鏈中產(chǎn)生出來(lái)的樣本服從我們的目標(biāo)分布，那至少要讓這些樣本都是從同一個(gè)分布中出來(lái)的。拆解向量和矩陣乘法的計(jì)算規(guī)則，會(huì)發(fā)現(xiàn)下面的式子是p為平穩(wěn)分布的充要條件

?（式1）

一個(gè)充分但是不必要的條件是下面這個(gè)式子，這個(gè)式子也叫做細(xì)致平穩(wěn)條件（detailed balance）

?（式2）

我們可以證明滿足細(xì)致平穩(wěn)條件（式2）的分布，都滿足（式1），因此都是平穩(wěn)分布，證明如下：

（c) 轉(zhuǎn)移矩陣的構(gòu)造

所以我們現(xiàn)在要做的，就是對(duì)于我們想要采樣的分布p作為初始分布，構(gòu)造一個(gè)滿足上述細(xì)致平穩(wěn)條件的轉(zhuǎn)移矩陣??，利用這個(gè)轉(zhuǎn)移矩陣來(lái)生成一系列樣本，那么每個(gè)樣本也會(huì)服從這個(gè)分布p。

假設(shè)我們有一個(gè)轉(zhuǎn)移矩陣，但是這個(gè)矩陣不符合細(xì)致平穩(wěn)條件，也就是

我們可以通過(guò)一個(gè)改造來(lái)讓細(xì)致平穩(wěn)條件得到滿足，我們?cè)趦蛇叿謩e乘以??和?

?（式3）

其中

把上式代入（式3）這樣等式的兩邊就完全一樣了，自然細(xì)致平穩(wěn)條件也就滿足了

也就是說(shuō)我們新的轉(zhuǎn)移矩陣?需要修改為

我們就得到了一個(gè)讓細(xì)致平穩(wěn)條件得到滿足的轉(zhuǎn)移矩陣。所以我們就可以選擇一個(gè)我們已經(jīng)能采樣的分布（例如高斯分布）作為??，然后再以??的概率接受這個(gè)采樣的結(jié)果。如此反復(fù)，就能得到一個(gè)序列的樣本，且樣本符合??的分布。看起來(lái)勝利就在不遠(yuǎn)處。

這里還有一個(gè)小問(wèn)題，我們的初始分布一開(kāi)始并不是符合??的話，那么經(jīng)過(guò)概率轉(zhuǎn)移以后生成的下一個(gè)隨機(jī)變量的分布也不符合。還好如果是齊次馬爾可夫鏈，滿足一些很簡(jiǎn)單的條件(不是本文重點(diǎn)，不增加理解復(fù)雜度不展開(kāi)了，感興趣的讀者google之)，就能成為一個(gè)“各態(tài)歷經(jīng)”的馬爾可夫鏈。不管初始分布是怎么樣的，只要對(duì)應(yīng)的轉(zhuǎn)移矩陣??與某一個(gè)分布??滿足了細(xì)致平穩(wěn)條件，在經(jīng)過(guò)一定次數(shù)的概率轉(zhuǎn)移后（即m趨向于無(wú)窮大之后），產(chǎn)生的分布會(huì)逼近唯一的平穩(wěn)分布。

另外一個(gè)問(wèn)題是這里還是有一個(gè)的拒絕采樣，如果這個(gè)值比較小，采樣效率還是會(huì)比較低。不過(guò)有個(gè)小trick，可以讓這個(gè)概率大一些。我們?cè)诩?xì)致平穩(wěn)條件的等式兩邊都乘以同一個(gè)數(shù)C，等式仍然成立

也就是說(shuō)我們可以用??代替原來(lái)的作為拒絕采樣的接受率，我們想辦法讓盡可能大就好了，但是作為概率值，也不能超過(guò)1。也就是說(shuō)

?且?

所以C最大可以取到

這樣，新的接受率

回憶一下，我們還需要考慮的一個(gè)問(wèn)題是分布的歸一化分母難計(jì)算的問(wèn)題，這里也剛好就解決了，因?yàn)榉肿臃帜钢卸己袣w一化分母，自然就約掉了

(d) Metropolis-Hastings算法

到目前為止，我們就得到了一種效率相對(duì)比較高的MCMC采樣方式，叫做Metropolis-Hastings算法（簡(jiǎn)寫(xiě)M-H算法）?？偨Y(jié)一下M-H算法如下：

1. 選擇一個(gè)我們已經(jīng)能采樣的分布（例如高斯分布）作為?。(例如??）
2. 隨機(jī)獲得一個(gè)初始樣本??，令?
3. 根據(jù)??采樣出新的樣本??，以??的概率接受采樣出的新樣本，如果被拒絕（??的概率），則?。令?
4. 重復(fù)3的操作，得到??，選擇??之后的樣本（分布收斂到平穩(wěn)分布之后）作為采樣的最終結(jié)果。其中M是一個(gè)超參數(shù)。

如果我們選擇的T函數(shù)，滿足對(duì)稱性（因?yàn)門(mén)是我們自己選的，很容易滿足），即

則有?

這種特殊的M-H算法，也叫Metropolis Sampling算法。（這個(gè)算法在1953年就已經(jīng)誕生了，而M-H是在1970年在這個(gè)基礎(chǔ)上進(jìn)行了通用化）

（e) 吉布斯采樣（Gibbs Sampling)

雖然這里??已經(jīng)被我們專門(mén)放大過(guò)了，但是在緯度比較高的時(shí)候，還是可能會(huì)出現(xiàn)接受率比較低的情況。有沒(méi)有辦法把這個(gè)接受率始終變成1呢？有的，比如吉布斯采樣（Gibbs Sampling）方法。

在后面的推導(dǎo)中，我們需要把隨機(jī)變量的每一維區(qū)分出來(lái)，有如下表示方法

其中??表示隨機(jī)變量?的第1個(gè)緯度，表示的除了第k維的其他緯度。

吉布斯采樣是M-H采樣的一個(gè)特例。這里我們先給出吉布斯采樣的轉(zhuǎn)移矩陣??，再證明利用該轉(zhuǎn)移矩陣，會(huì)讓接受率為1。如果當(dāng)前隨機(jī)變量是??，馬爾可夫鏈下一個(gè)隨機(jī)變量??的時(shí)候，如果隨機(jī)變量和隨機(jī)變量??除了第k維，其他緯度都一樣，也就是說(shuō)有?，則轉(zhuǎn)移概率為目標(biāo)分布??的條件概率分布??。對(duì)于不滿足這個(gè)條件的??，轉(zhuǎn)移矩陣為0。

同理有

接下來(lái)，我們來(lái)證明按上面這個(gè)轉(zhuǎn)移矩陣采樣樣本后，接受率為1?；仡櫼幌翸-H中的接受率公式，再根據(jù)概率公式展開(kāi)下有

因?yàn)樾枰慌袛嘟邮芘c否的樣本，都是從上述的特殊的??中采樣出來(lái)的，都會(huì)滿足?（否則轉(zhuǎn)移概率是0，不會(huì)被采樣出來(lái)），所以有??，??，代入則有

同理，采樣出的樣本都會(huì)滿足??，代入得到

所以只要按上述的轉(zhuǎn)移矩陣??來(lái)采樣，接受率是100%。這里又還有一個(gè)小問(wèn)題沒(méi)解決，即k值是怎么得到的？假設(shè)緯度為n，則k是從1到n中隨機(jī)得到的。嚴(yán)謹(jǐn)一些的話，需要把k值產(chǎn)生的概率，也加到??里，在上面接受率的證明里，很快就會(huì)被消元，不影響結(jié)論。不過(guò)在實(shí)際應(yīng)用的時(shí)候，會(huì)讓k從1到n輪流取得，方便實(shí)現(xiàn)。

美中不足的是，在普通的M-H里，我們可以任意選擇可采樣的分布（例如高斯分布）作為??，在吉布斯采樣中不行了，采樣分布必須和目標(biāo)分布掛鉤。

總結(jié)一下吉布斯采樣（Gibbs Sampling)算法如下：

1. 根據(jù)我們的目標(biāo)分布??，得到每個(gè)緯度的條件分布?
2. 隨機(jī)獲得一個(gè)初始樣本??，初始化k=1，m=1
3. 根據(jù)每個(gè)緯度的條件分布??采樣得到??中第k個(gè)緯度的值??，對(duì)于其他緯度則直接復(fù)制??的值，即?，這樣得到?
4. k=k+1，如果k>n，則重置k=1；m=m+1
5. 重復(fù)3，4的操作，得到??，選擇??之后的樣本（分布收斂到平穩(wěn)分布之后）作為采樣的最終結(jié)果。其中M是一個(gè)超參數(shù)。

要從的每個(gè)緯度（??）的條件分布中采樣也很不容易。在一些具體應(yīng)用的時(shí)候，通常借助有向概率圖（也叫貝葉斯網(wǎng)絡(luò)，注意不是貝葉斯神經(jīng)網(wǎng)絡(luò)）或無(wú)向概率圖（也叫馬爾可夫隨機(jī)場(chǎng)）的工具來(lái)去掉一些條件依賴，以及巧妙地選擇共軛先后驗(yàn)分布，來(lái)解決這個(gè)問(wèn)題。例如經(jīng)典的用吉布斯采樣解LDA，吉布斯采樣解受限玻爾茲曼機(jī)（RBM）的例子。有興趣的同學(xué)可以自行再深入了解下。

所以，如果不對(duì)的分布做任何限制，M-H的采樣方法是更通用地從采樣的方法。不過(guò)缺點(diǎn)就是會(huì)有高維空間采樣效率較低的問(wèn)題。所以下面的代碼以M-H采樣為例（其實(shí)是M-H的特例Metropolis Sampling）。

M-H采樣訓(xùn)練貝葉斯神經(jīng)網(wǎng)絡(luò)的代碼樣例：

（說(shuō)明：樣例代碼的原則是能說(shuō)明算法原理，追求可讀性，所以不會(huì)考慮可擴(kuò)展性，性能等因素。為了兼容用pytorch丹爐的朋友，訓(xùn)練方式用和pytorch更類似的接口。運(yùn)行環(huán)境為python 3，tf 2.3）

1. 構(gòu)造一些樣本，用來(lái)后面訓(xùn)練和展現(xiàn)效果。（此處參考了這篇文章中的樣本產(chǎn)生和部分代碼，鏈接：krasserm.github.io/2019）

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

def f(x, sigma):
return 10 * np.sin(2 * np.pi * (x)) + np.random.randn(*x.shape) * sigma

num_of_samples = 64 # 樣本數(shù)
noise = 1.0 # 噪音規(guī)模

X = np.linspace(-0.5, 0.5, num_of_samples).reshape(-1, 1)
y_label = f(X, sigma=noise) # 樣本的label
y_truth = f(X, sigma=0.0) # 樣本的真實(shí)值

plt.scatter(X, y_label, marker='+', label='Training data')
plt.plot(X, y_truth, label='Ground Truth')
plt.title('Noisy training data and ground truth')
plt.legend();

2. 畫(huà)出的圖中，實(shí)線是y值的真實(shí)分布，“+”號(hào)是加上一定噪音后，我們觀測(cè)得到的樣本，也是后面訓(xùn)練用的樣本。

3. 先驗(yàn)分布選擇了一個(gè)由兩個(gè)高斯分布組成的混合高斯分布

?，

轉(zhuǎn)移矩陣T選擇的是高斯分布??，可以發(fā)現(xiàn)這個(gè)分布對(duì)于交換i和j后，值不變，也就是對(duì)稱的。消去接受率公式分子和分母中的 T項(xiàng)，所以有接受率??，即為M-H的特例：Metropolis Sampling。

from tensorflow.keras.activations import tanh
import tensorflow as tf

class BNN_MH():

def __init__(self, prior_sigma_1=1.5, prior_sigma_2=0.1, prior_pi=0.5):
# 先驗(yàn)分布假設(shè)的各種參數(shù)
self.prior_sigma_1 = prior_sigma_1
self.prior_sigma_2 = prior_sigma_2
self.prior_pi_1 = prior_pi
self.prior_pi_2 = 1.0 - prior_pi

# 模型中的所有參數(shù)，即theta=(w0,b0,w1,b1,w2,b2)
    w0, b0 = self.init_weights([1, 5])
    w1, b1 = self.init_weights([5, 10])
    w2, b2 = self.init_weights([10, 1])
self.w_list = [w0, w1, w2]
self.b_list = [w0, b1, b2]

def init_weights(self, shape):
    w = tf.Variable(tf.random.truncated_normal(shape, mean=0., stddev=1.))
    b = tf.Variable(tf.random.truncated_normal(shape[1:], mean=0., stddev=1.))
return w, b

def forward(self, X, w_list, b_list):
# 簡(jiǎn)單的3層神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)
    x = tanh(tf.matmul(X, w_list[0]) + b_list[0])
    x = tanh(tf.matmul(x, w_list[1]) + b_list[1])
return tf.matmul(x, w_list[2]) + b_list[2]

def generate_weights_by_gaussian_random_walk(self, step_size = 0.1):
# 高斯隨機(jī)游走，即從T（x^i -> x^j）= Gaussian(x^i - x^j|0,1.0)中根據(jù)已知x^i采樣x^j
    next_w0 = self.w_list[0] + tf.random.normal(self.w_list[0].shape)*step_size
    next_b0 = self.b_list[0] + tf.random.normal(self.b_list[0].shape)*step_size
    next_w1 = self.w_list[1] + tf.random.normal(self.w_list[1].shape)*step_size
    next_b1 = self.b_list[1] + tf.random.normal(self.b_list[1].shape)*step_size
    next_w2 = self.w_list[2] + tf.random.normal(self.w_list[2].shape)*step_size
    next_b2 = self.b_list[2] + tf.random.normal(self.b_list[2].shape)*step_size

# 存儲(chǔ)x^j
self.next_w_list = [next_w0, next_w1, next_w2]
self.next_b_list = [next_b0, next_b1, next_b2]   

def prior(self, w):
# 先驗(yàn)分布假設(shè)
return self.prior_pi_1 * self.gaussian_pdf(w, 0.0, self.prior_sigma_1) \
+ self.prior_pi_2 * self.gaussian_pdf(w, 0.0, self.prior_sigma_2)

def gaussian_pdf(self, x, mu, sigma):
return tf.compat.v1.distributions.Normal(mu,sigma).prob(x)    

def accept_rate(self, X, y_label):   

# 計(jì)算先驗(yàn)分布的比率
    prior_ratio = 1.0
for w, w_next, b, b_next in zip(self.w_list, self.next_w_list, \
self.b_list, self.next_b_list):                            
        p_theta = self.prior(w)*self.prior(b)
        p_theta_next = self.prior(w_next)*self.prior(b_next)
        prior_ratio *= tf.reduce_prod(tf.divide(p_theta_next,p_theta))

# 計(jì)算似然率的比率
    y_pred = self.forward(X, self.w_list, self.b_list)
    p_d_theta = self.gaussian_pdf(y_label, y_pred, 1.0)
    y_pred_next = self.forward(X, self.next_w_list, self.next_b_list)
    p_d_theta_next = self.gaussian_pdf(y_label, y_pred_next, 1.0) 

    likelihood_ratio = tf.reduce_prod(tf.divide(p_d_theta_next,p_d_theta))

# 接受率等于先驗(yàn)分布的比率*似然率的比率
    ac_rate = prior_ratio * likelihood_ratio
return min(ac_rate,1.0)

def train(self, X, y_label):
self.saved_weights = [] # 存所有w和b值的地方
self.ac_rates = [] # 保存接受率，監(jiān)控用
    num_of_rounds = 10000
    burn_in_rounds = 1000
for i in range(0,num_of_rounds):
self.generate_weights_by_gaussian_random_walk()
        ac_rate = self.accept_rate(X, y_label)
self.ac_rates.append(ac_rate)
if tf.random.uniform([1]) < ac_rate:
self.w_list = self.next_w_list
self.b_list = self.next_b_list
if i >= burn_in_rounds:
self.saved_weights.append({"w_list":self.w_list, "b_list":self.b_list})
return self.ac_rates

def predict(self, X):
# 用所有存儲(chǔ)的參數(shù)，分別都預(yù)估下模型的輸出值
return list(map(lambda saved_weight : self.forward(X, saved_weight["w_list"], \
                                    saved_weight["b_list"]), self.saved_weights))

X = X.astype('float32')
y_label = y_label.astype('float32')

model = BNN_MH()
ac_rates = model.train(X,y_label)
plt.plot(ac_rates)
plt.grid()
X_test = np.linspace(-1., 1., num_of_samples * 2).reshape(-1, 1).astype('float32')
y_preds = model.predict(X_test)
y_preds = np.concatenate(y_preds, axis=1)

plt.scatter(X, y_label, marker='+', label='Training data')
plt.plot(X_test, np.mean(y_preds, axis=1), 'r-', label='Predict Line')
plt.fill_between(X_test.reshape(-1), np.percentile(y_preds, 2.5, axis=1), np.percentile(y_preds, 97.5, axis=1), color='r', alpha=0.3, label='95% Confidence')
plt.grid()
plt.legend()

4. 畫(huà)出來(lái)的曲線是每一輪的接受率

5. 預(yù)估階段，我們把預(yù)估的均值和不確定度都畫(huà)出來(lái)

X_test = np.linspace(-1., 1., num_of_samples * 2).reshape(-1, 1).astype('float32')
y_preds = model.predict(X_test)
y_preds = np.concatenate(y_preds, axis=1)

plt.scatter(X, y_label, marker='+', label='Training data')
plt.plot(X_test, np.mean(y_preds, axis=1), 'r-', label='Predict Line')
plt.fill_between(X_test.reshape(-1), np.percentile(y_preds, 2.5, axis=1), np.percentile(y_preds, 97.5, axis=1), color='r', alpha=0.3, label='95% Confidence')
plt.grid()
plt.legend()

6. 圖中，紅色線條就是所有BNN輸出值的均值，也就作為最終的預(yù)估值。而紅色的區(qū)域?qū)捳?，則反應(yīng)了所有BNN輸出值的不確定度（為了方便可視化這里用分位數(shù)）。可以看到，對(duì)于沒(méi)有樣本的區(qū)域，不確定度較大，而對(duì)于樣本比較密集的地方，不確定度小。另外，在樣本有覆蓋的領(lǐng)域，轉(zhuǎn)彎的地方，因?yàn)橐x原來(lái)的路線，不確定度大；而直線的地方，則不確定度更小。

在（下）篇中，我們會(huì)介紹另外一種貝葉斯推斷的近似解法，變分推斷。也會(huì)介紹一種實(shí)用的BNN實(shí)現(xiàn)MC Dropout，以及討論實(shí)際應(yīng)用中的一些問(wèn)題。

還是廣告

又到了招聘旺季，我們正在大力尋找志同道合的朋友一起在某手商業(yè)化做最有趣最前沿的廣告算法，初中高級(jí)廣告算法職位均有HC（迅速上車(chē)，還能趕上上市）。感興趣的朋友歡迎加我的個(gè)人微信約飯約咖啡索要JD或發(fā)送簡(jiǎn)歷。（產(chǎn)品和運(yùn)營(yíng)也在招人，看機(jī)會(huì)的朋友我可以幫忙直接推薦給leader們。）

作者個(gè)人微信（添加注明申探社讀者及簡(jiǎn)單介紹）：

歡迎掃描下面二維碼關(guān)注本公眾號(hào)，也歡迎關(guān)注知乎“申探社”專欄