快捷導(dǎo)航

python實(shí)現(xiàn)貝葉斯推斷的例子

更新時(shí)間：2021年09月10日 11:31:03 作者：笨牛慢耕

本文介紹一個(gè)貝葉斯推斷的python實(shí)現(xiàn)，并展現(xiàn)了基于標(biāo)量運(yùn)算的實(shí)現(xiàn)和基于numpy的矩陣運(yùn)算的實(shí)現(xiàn)之間的差別，感興趣的可以了解一下

1. 前言

本文介紹一個(gè)貝葉斯推斷的python實(shí)現(xiàn)例，并展現(xiàn)了基于標(biāo)量運(yùn)算的實(shí)現(xiàn)和基于numpy的矩陣運(yùn)算的實(shí)現(xiàn)之間的差別。

2. 問(wèn)題描述

本問(wèn)題例取自于Ref1-Chapter1.

問(wèn)題描述：假設(shè)有一個(gè)制作燈泡的機(jī)器。你想知道機(jī)器是正常工作還是有問(wèn)題。為了得到答案你可以測(cè)試每一個(gè)燈泡，但是燈泡數(shù)量很多，每一個(gè)都測(cè)試在實(shí)際生產(chǎn)過(guò)程中可能是無(wú)法承受的。使用貝葉斯推斷，你可以基于少量樣本（比如說(shuō)抽檢結(jié)果）來(lái)估計(jì)機(jī)器是否在正常地工作（in probabilistic way）。

構(gòu)建貝葉斯推斷時(shí)，首先需要兩個(gè)要素:

(1) 先驗(yàn)分布

(2) 似然率

先驗(yàn)分布是我們關(guān)于機(jī)器工作狀態(tài)的初始信念。首先我們確定第一個(gè)刻畫(huà)機(jī)器工作狀態(tài)的隨機(jī)變量，記為M。這個(gè)隨機(jī)變量有兩個(gè)工作狀態(tài)：{working, broken}，以下簡(jiǎn)寫(xiě)成{w, br}(縮寫(xiě)成br是為了與下面的Bad縮寫(xiě)成b區(qū)分開(kāi)來(lái))。作為初始信念，我們相信機(jī)器是好的，是可以正常工作的，定義先驗(yàn)分布如下：

P(M=working) = 0.99

P(M=broken ) = 0.01

這表明我們對(duì)于機(jī)器正常工作的信念度很高，有99%的概率能夠正常工作。

第二個(gè)隨機(jī)變量是L，表示機(jī)器生產(chǎn)的燈泡的工作狀態(tài)。燈泡可能是好，也可能是壞的，包含兩個(gè)狀態(tài)：{good, bad}，以下簡(jiǎn)寫(xiě)成{g,b}，注意br與b的區(qū)別。

我們需要基于機(jī)器工作狀態(tài)給出L的先驗(yàn)分布，也就是條件概率P(L|M)，在貝葉斯公式中它代表似然概率(likelihood)。

定義這個(gè)似然概率分布（由于M和L各有兩種狀態(tài)，所以一共包含4個(gè)條件概率）如下：

P(L=Good|M=w) = 0.99

P(L=Bad |M=w) = 0.01

P(L=Good|M=br ) = 0.6

P(L=Bad |M=br ) = 0.4

以上似然概率表明，在機(jī)器正常時(shí)我們相信每生成100個(gè)燈泡只會(huì)有一個(gè)壞的，而機(jī)器不正常時(shí)也不是所有燈泡都是壞的，而是有40%會(huì)是壞的。為了實(shí)現(xiàn)的方便，可以寫(xiě)成如下的矩陣形式：

現(xiàn)在，我們已經(jīng)完整地刻畫(huà)了貝葉斯模型，可以用它來(lái)做一些神奇的估計(jì)和預(yù)測(cè)的工作了。

我們的輸入是一些燈泡的抽檢結(jié)果。假設(shè)我們抽檢了十個(gè)燈泡其抽檢結(jié)果如下：

{bad, good, good, good, good, good, good, good, good, good}

讓我們來(lái)看看基于貝葉斯推斷的我們對(duì)于機(jī)器工作狀態(tài)的信念（后驗(yàn)概率）如何變化。

3. 貝葉斯規(guī)則

貝葉斯推斷規(guī)則以貝葉斯公式的形式表示為：

具體映射到本問(wèn)題中可以表達(dá)如下：

貝

葉斯推斷的優(yōu)先在于可以以在線（online）的方式進(jìn)行，即觀測(cè)數(shù)據(jù)可以一個(gè)一個(gè)地到來(lái)，每次受到一個(gè)新的觀測(cè)數(shù)據(jù)，就進(jìn)行一次基于貝葉斯公式的后驗(yàn)概率的計(jì)算更新，而更新后的后驗(yàn)概率又作為下一貝葉斯推斷的先驗(yàn)概率使用。因此在線的貝葉斯推斷的基本處理流程如下所示：

4. Bayes engine: scalar implementation

首先，我們以標(biāo)量運(yùn)算的方式寫(xiě)一個(gè)函數(shù)來(lái)進(jìn)行bayes推斷處理。

prior以向量的形式存儲(chǔ)先驗(yàn)概率分布，prior[0]表示P(M=working)，prior[1]表示P(M=broken)。

likelihood以矩陣的形式方式存儲(chǔ)似然概率分布。其中第1行表示P(L/M=working)，第2行表示P(L/M=broken).

在本例中，當(dāng)輸入時(shí)，evidence的計(jì)算式（注意evidence是依賴于輸入的觀測(cè)數(shù)據(jù)的）是:

注意，當(dāng)我寫(xiě)P(w)其實(shí)是表示P(M=w)，而其實(shí)是表示，余者類推。根據(jù)上下文，這些應(yīng)該不會(huì)導(dǎo)致混淆。

第一個(gè)函數(shù)的代碼如下：

def bayes_scalar(prior, likelihood, data):
    """
    Bayesian inference function example.
    Parameters
    ----------
    prior : float, 1-D vector
        prior information, P(X).
    likelihood : float 2-D matrix
        likelihood function, P(Y|X).
    data : List of strings. Value: 'Good','Bad'
        Observed data samples sequence
    Returns
    -------
    posterior : float
        P(X,Y), posterior sequence.
    """
    
    posterior = np.zeros((len(data)+1,2))
    posterior[0,:] = prior  # Not used in computation, just for the later plotting
 
    for k,L in enumerate(data):
        if L == 'good':
            L_value = 0
        else:
            L_value = 1
        #print(L, L_value, likelihood[:,L_value])
 
        evidence      = likelihood[0,L_value] * prior[0] + likelihood[1,L_value] * prior[1]
        LL0_prior_prod= likelihood[0,L_value] * prior[0]                
        posterior[k+1,0]  = LL0_prior_prod / evidence
 
        LL1_prior_prod= likelihood[1,L_value] * prior[1]        
        posterior[k+1,1]  = LL1_prior_prod / evidence
        
        prior = posterior[k+1,:] # Using the calculated posterior at this step as the prior for the next step
                
    return posterior

5. Bayes engine: vectorization

我們注意到，evidence的計(jì)算可以表示成兩個(gè)向量的點(diǎn)積，如下所示。

這樣就非常方便用numpy來(lái)實(shí)現(xiàn)了。本例中每個(gè)隨機(jī)變量只有兩種取值，在復(fù)雜的情況下，每個(gè)隨機(jī)變量有很多種取值時(shí)，有效利用向量或矩陣的運(yùn)算是簡(jiǎn)潔的運(yùn)算實(shí)現(xiàn)的必不可缺的要素。以上這兩個(gè)向量的點(diǎn)積可以用numpy.dot()來(lái)實(shí)現(xiàn)。

另外，likelihood和prior的乘積是分別針對(duì)M的兩種狀態(tài)進(jìn)行計(jì)算（注意，我們需要針對(duì)M的兩種不同狀態(tài)分別計(jì)算posterior），不是用向量的點(diǎn)積進(jìn)行計(jì)算，而是一種element-wise multiplication，可以用numpy.multiply()進(jìn)行計(jì)算。所以在vectorization版本中貝葉斯更新處理削減為兩條語(yǔ)句，與上面的scalar版本相比顯得非常優(yōu)雅簡(jiǎn)潔（好吧，也許這個(gè)簡(jiǎn)單例子中還顯不出那么明顯的優(yōu)勢(shì)，但是隨著問(wèn)題的復(fù)雜度的增加，這種優(yōu)勢(shì)就會(huì)越來(lái)越明顯了。）

由此我們得到向量化處理的函數(shù)如下：

def bayes_vector(prior, likelihood, data):
    """
    Bayesian inference function example.
    Parameters
    ----------
    prior : float, 1-D vector
        prior information, P(X).
    likelihood : float 2-D matrix
        likelihood function, P(Y|X).
    data : List of strings. Value: 'Good','Bad'
        Observed data samples sequence
    Returns
    -------
    posterior : float
        P(X,Y), posterior sequence.
    """
    
    posterior = np.zeros((len(data)+1,2))
    posterior[0,:] = prior  # Not used in computation, just for the later plotting
 
    for k,L in enumerate(data):
        if L == 'good':
            L_value = 0
        else:
            L_value = 1
        #print(L, L_value, likelihood[:,L_value])
 
        evidence          = np.dot(likelihood[:,L_value], prior[:])
        posterior[k+1,:]  = np.multiply(likelihood[:,L_value],prior)/evidence
 
        prior = posterior[k+1,:] # Using the calculated posterior at this step as the prior for the next step
        
    return posterior

6. 測(cè)試

讓我們來(lái)看看利用以上函數(shù)對(duì)我們的觀測(cè)數(shù)據(jù)進(jìn)行處理，后驗(yàn)概率將會(huì)如何變化。

import numpy as np
import matplotlib.pyplot as plt
 
prior      = np.array([0.99,0.01])
likelihood = np.array([[0.99,0.01],[0.6,0.4]])
data       = ['bad','good','good','good','good','good','good','good']        
 
posterior1 = bayes_scalar(prior,likelihood,data)
posterior2 = bayes_vector(prior,likelihood,data)
 
if np.allclose(posterior1,posterior2):
    print('posterior1 and posterior2 are identical!')
 
fig, ax = plt.subplots()
ax.plot(posterior1[:,0])
ax.plot(posterior1[:,1])
ax.grid()
# fig.suptitle('Poeterior curve vs observed data')
ax.set_title('Posterior curve vs observed data')
plt.show()

運(yùn)行以上代碼可以得到后驗(yàn)概率的變化如下圖所示（注意第一個(gè)點(diǎn)是prior）：