申探社：不確定度建模（中）-蒙特卡洛采樣

3050 閱讀 0 評論 8 點贊

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

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

第一步：假設先驗??，對似然率建模?
第二步：計算后驗?
展開后有?
第三步：計算預測值的分布?
第四步：計算y的期望值??；
如果需要，計算y的方差??作為預估值的不確定度

通常我們很容易的可以把后驗分布的表達式??寫出來，但是到第三步，是一個積分，這里就比較麻煩了。如果我們選擇的似然率的模型和先驗模型不是剛好有些特殊性質可以利用，很難直接把積分求出來。另外第四步里也需要用積分來求期望值（或者方差）。

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

一、蒙特卡洛采樣

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

如果我們是蒙特卡洛賭城的老板，我們需要設計一臺老虎機，怎么做可以知道這臺老虎機是不是穩賺不賠的？一個最簡單的方法就是找一批測試人員輪流模擬玩家在這臺機器上試玩1000次，看看最終機器是不是賺錢了，賺了多少錢，然后再除以1000，得到的就是玩家平均每玩一次賭場賺的錢。這個叫做老虎機的ev（expected value）。賭場的機器，對于玩家來說都是負ev的。所以長期玩下去，賭場肯定是賺的。

這種求單次賭場平均盈利的方法，就是蒙特卡洛采樣法。更通用地來說，蒙特卡洛采樣法是指用從一個分布中采樣很多樣本，然后用這些樣本來代替這個分布進行相關計算（例如求期望值）的方法。比如說我們要下面這個積分

回顧一下期望值計算公式，上式其實就是當x的分布為??時，的期望值。

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

我們先回到貝葉斯推斷里關鍵的第三步和第四步，把第三步的分布帶入到第四步的期望值公式里，則有??。因為y與??無關，可以放到第二個積分內部，則有。再調整下積分的順序，并把與y無關的??挪到外面，則有?

上式中，我們會發現括號內的積分，其實就是?的均值。在大多數時候，我們在對似然率建模的時候，會間接通過對均值建模來實現。例如在回歸問題里，神經網絡的輸出作為高斯分布的均值來建模似然率，即??，其中??表示神經網絡的輸入和參數，??表示神經網絡的輸出值，也是高斯分布的均值。?一般是個超參數。如果是分類問題，則神經網絡的輸出會作為伯努利分布里的y=1的概率，也是伯努利分布的均值。所以上面的式子就變成?

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

那么現在的問題就變成：我們如何根據一個分布的表達式（在這里主要是指），從這個分布中采樣出很多樣本來。（采樣問題）

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

二、基本分布的采樣方法

對于均勻分布，我們可以用隨機數發生器很容易的采樣出來。借助隨機數發生器，我們可以實現另外一些簡單的分布，例如伯努利分布，用簡單的if else就可以實現；也可以實現稍微復雜一些的分布，比如高斯分布，可以借助Box-Muller方法來實現：

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

因為不是本文重點，不展開證明。

但是對于更復雜的分布，采樣就變得并不容易了。有一系列采樣算法來幫助我們通過上面介紹的這些基本分布的采樣方法，來實現從我們想要的分布中去采樣樣本。

三、復雜分布的采樣算法

拒絕采樣

圖片來源《Pattern Recognition and Machine Learning》

假設如果我們開了一家飯店，開業大酬賓搞了一個促銷活動，來吃飯的人只要扔一次骰子，如果扔6就能享受免單。那么免單概率是1/6。如果我們發現這樣下去要虧本，想把免單概率調整為1/12怎么辦？那么我們可以修改規則為對扔到6的人，隨機一半概率才能免單，就可以把概率變為1/12了。具體實施方式為，要求扔到6的人，再扔一次硬幣，如果是反面，則我們拒絕免單，如果是正面，我們接受免單。這樣就達到了以1/12概率采樣的效果。

拒絕采樣也是類似的思想。假設我們已經能從分布q中采樣樣本（這個分布也叫Proposal Distribution，例如一個高斯分布），而我們的目標是要對目標分布 p去采樣。

在一開始，我們先對兩個分布都做一個變換。先把p分布乘以一個歸一化分母??（是個常數），得到??。回顧上面說的??有一個歸一化分母不好計算的問題，而做拒絕采樣只需要用到??，因此就不需要考慮歸一化分母計算的問題了。然后我們再對分布q的概率密度函數乘以一個常數k，使得新的概率密度函數始終在的上方，也就是讓藍線始終在紅線上方或者和紅線重合。

然后我們從分布q中進行采樣，假設我們采樣出的樣本是??，它對應的放大后的q分布的概率密度值是??，它對應的目標分布的概率密度是??，那么拒絕概率就是??，或者說接受概率是，這樣最后被接受的樣本，就會符合??的分布。（PS：嚴格來說因為被放大了k倍，放大了倍，都已經不是概率函數了，因為積分不為1）

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

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

圖片來源為電視劇《老公們的私房錢》

2. 重要性采樣

重要性采樣其實并不是真的采樣，或者說可以認為是一種軟采樣。我們并沒有真的拒絕掉或接受一些樣本，而是通過把這個樣本的權重調低/調高來表達對樣本的偏好（拒絕采樣中被拒絕的樣本可以認為樣本權重被置0了）。

現在我們要用從p(x)里采樣的樣本來計算f(x)的均值，而我們已經實現的是從q(x)中采樣（例如q(x)是個高斯分布）。那么我們拆解下f(x)均值計算的公式如下

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

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

樂隊的夏天中，超級樂迷有10票，相當于每個人一票，但是權重是10; 專業樂迷，每人兩票，相當于權重是2；大眾樂迷每人1票，相當于權重是1。這是因為組委會認為，超級樂迷和專業樂迷水平更高，他們的評審水平，更接近于作品的真實水平。

《樂隊的夏天》投票是280票制，四個超級樂迷40票（1人10票），20個專業樂迷40票（1人2票），200個大眾樂迷200票（1人1票）

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

解決了采樣的問題，我們來看看怎么在重要性采樣里解決的歸一化分母不好計算的問題。我們還是用表示沒有歸一化后的目標“概率密度”函數，用??表示沒有歸一化后的proposal distribution的“概率密度”函數。把這兩個歸一化分母提到前面，則有

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

上式的倒數，就是要求的。即對于采樣樣本??，計算??的均值，再求一個倒數就可以了。至此，歸一化問題也被解決了。

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

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

但是我們需要一點耐心，先從馬爾可夫鏈慢慢聊起。一階馬爾可夫鏈是指一個序列的隨機變量??，他們滿足下面的式子

也就是說，??只和序列里上一個隨機變量??有關系，而和之前的其他隨機變量（例如??)都沒有關系。

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

（a) 齊次馬爾科夫鏈

如果對于所有的m取值，轉移概率??對于任意的取值和的取值，都相同，則我們說這個馬爾可夫鏈是齊次馬爾可夫鏈。也就是說對于任意??，??(注意這里用加括號表示馬爾可夫鏈中的隨機變量，不加表示普通的隨機變量)，齊次馬爾可夫鏈都有?

可以理解為在根據每個隨機變量的值產生序列里的下一個隨機變量的值時，遵循著同樣的一個轉移概率矩陣，和在序列里的位置無關。

如果我們把轉移概率表示為??，齊次馬爾可夫鏈的轉移概率因為對所有m都相同，也就是和m的取值無關，可以計為?，即?

（b) 概率矩陣，平穩分布與細致平穩條件

如果我們把對應的轉移概率放到矩陣??里，使得矩陣里每個元素，則矩陣?就是齊次馬爾可夫鏈的轉移概率矩陣。其中表示的是從狀態??轉移到狀態??的概率。

還是剛才的例子，假設隨機變量的取值范圍是離散的集合??，我們把當前隨機變量x等于不同取值的概率寫到一個長度為n向量里。

對于一個以向量表示的分布乘以一個矩陣，即?，物理意義就是該分布經過一次由矩陣T定義的轉移之后，產生的序列里下一個隨機變量的概率分布。如果有

也就說下一個隨機變量的概率分布和當前一致，沒有發生變化。我們說這個分布是關于這個馬爾可夫鏈的平穩分布。回憶一下我們的目標，是要讓從馬爾可夫鏈中產生出來的樣本服從我們的目標分布，那至少要讓這些樣本都是從同一個分布中出來的。拆解向量和矩陣乘法的計算規則，會發現下面的式子是p為平穩分布的充要條件

?（式1）

一個充分但是不必要的條件是下面這個式子，這個式子也叫做細致平穩條件（detailed balance）

?（式2）

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

（c) 轉移矩陣的構造

所以我們現在要做的，就是對于我們想要采樣的分布p作為初始分布，構造一個滿足上述細致平穩條件的轉移矩陣??，利用這個轉移矩陣來生成一系列樣本，那么每個樣本也會服從這個分布p。

假設我們有一個轉移矩陣，但是這個矩陣不符合細致平穩條件，也就是

我們可以通過一個改造來讓細致平穩條件得到滿足，我們在兩邊分別乘以??和?

?（式3）

其中

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

也就是說我們新的轉移矩陣?需要修改為

我們就得到了一個讓細致平穩條件得到滿足的轉移矩陣。所以我們就可以選擇一個我們已經能采樣的分布（例如高斯分布）作為??，然后再以??的概率接受這個采樣的結果。如此反復，就能得到一個序列的樣本，且樣本符合??的分布。看起來勝利就在不遠處。

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

另外一個問題是這里還是有一個的拒絕采樣，如果這個值比較小，采樣效率還是會比較低。不過有個小trick，可以讓這個概率大一些。我們在細致平穩條件的等式兩邊都乘以同一個數C，等式仍然成立

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

?且?

所以C最大可以取到

這樣，新的接受率

回憶一下，我們還需要考慮的一個問題是分布的歸一化分母難計算的問題，這里也剛好就解決了，因為分子分母中都含有歸一化分母，自然就約掉了

(d) Metropolis-Hastings算法

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

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

如果我們選擇的T函數，滿足對稱性（因為T是我們自己選的，很容易滿足），即

則有?

這種特殊的M-H算法，也叫Metropolis Sampling算法。（這個算法在1953年就已經誕生了，而M-H是在1970年在這個基礎上進行了通用化）

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

雖然這里??已經被我們專門放大過了，但是在緯度比較高的時候，還是可能會出現接受率比較低的情況。有沒有辦法把這個接受率始終變成1呢？有的，比如吉布斯采樣（Gibbs Sampling）方法。

在后面的推導中，我們需要把隨機變量的每一維區分出來，有如下表示方法

其中??表示隨機變量?的第1個緯度，表示的除了第k維的其他緯度。

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

同理有

接下來，我們來證明按上面這個轉移矩陣采樣樣本后，接受率為1。回顧一下M-H中的接受率公式，再根據概率公式展開下有

因為需要被判斷接受與否的樣本，都是從上述的特殊的??中采樣出來的，都會滿足?（否則轉移概率是0，不會被采樣出來），所以有??，??，代入則有

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

所以只要按上述的轉移矩陣??來采樣，接受率是100%。這里又還有一個小問題沒解決，即k值是怎么得到的？假設緯度為n，則k是從1到n中隨機得到的。嚴謹一些的話，需要把k值產生的概率，也加到??里，在上面接受率的證明里，很快就會被消元，不影響結論。不過在實際應用的時候，會讓k從1到n輪流取得，方便實現。

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

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

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

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

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

M-H采樣訓練貝葉斯神經網絡的代碼樣例：

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

1. 構造一些樣本，用來后面訓練和展現效果。（此處參考了這篇文章中的樣本產生和部分代碼，鏈接：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 # 樣本數
noise = 1.0 # 噪音規模

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) # 樣本的真實值

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. 畫出的圖中，實線是y值的真實分布，“+”號是加上一定噪音后，我們觀測得到的樣本，也是后面訓練用的樣本。

3. 先驗分布選擇了一個由兩個高斯分布組成的混合高斯分布

?，

轉移矩陣T選擇的是高斯分布??，可以發現這個分布對于交換i和j后，值不變，也就是對稱的。消去接受率公式分子和分母中的 T項，所以有接受率??，即為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):
# 先驗分布假設的各種參數
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

# 模型中的所有參數，即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):
# 簡單的3層神經網絡結構
    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):
# 高斯隨機游走，即從T（x^i -> x^j）= Gaussian(x^i - x^j|0,1.0)中根據已知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

# 存儲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):
# 先驗分布假設
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):   

# 計算先驗分布的比率
    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))

# 計算似然率的比率
    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))

# 接受率等于先驗分布的比率*似然率的比率
    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 = [] # 保存接受率，監控用
    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):
# 用所有存儲的參數，分別都預估下模型的輸出值
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. 畫出來的曲線是每一輪的接受率

5. 預估階段，我們把預估的均值和不確定度都畫出來


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輸出值的均值，也就作為最終的預估值。而紅色的區域寬窄，則反應了所有BNN輸出值的不確定度（為了方便可視化這里用分位數）?？梢钥吹?，對于沒有樣本的區域，不確定度較大，而對于樣本比較密集的地方，不確定度小。另外，在樣本有覆蓋的領域，轉彎的地方，因為要偏離原來的路線，不確定度大；而直線的地方，則不確定度更小。