(華南農(nóng)業(yè)大學(xué) 廣東 廣州 510000)

基于概率編程的貝葉斯推斷

林鎮(zhèn)溪

(華南農(nóng)業(yè)大學(xué)廣東廣州510000)

本文介紹了基于概率編程的貝葉斯推斷的建模過程和利用貝葉斯推斷對挑戰(zhàn)者號事故進(jìn)行分析和推斷，建立事故發(fā)生概率關(guān)于溫度的邏輯分布模型p(t)=1/(1+eβt+α)，其中α，β采用正態(tài)分布，最后采用伯努利分布將事故發(fā)生的概率和事故發(fā)生的結(jié)果聯(lián)系起來，使用PyMC計(jì)算后驗(yàn)分布，分析得出當(dāng)溫度為31度時事故發(fā)生的概率幾乎是必然的。

概率編程；貝葉斯推斷；事故發(fā)生概率；PyMC

一、實(shí)驗(yàn)背景

貝葉斯推斷是一種統(tǒng)計(jì)學(xué)方法，用來估計(jì)統(tǒng)計(jì)量的性質(zhì)。貝葉斯推斷與其他統(tǒng)計(jì)學(xué)推斷方法截然不同，他建立在主觀判斷的基礎(chǔ)上，也就是說，你不需要客觀證據(jù)，先估計(jì)一個值，然后根據(jù)實(shí)際結(jié)果不斷的修正。概率編程是對問題建模，然后利用計(jì)算機(jī)采樣的方法進(jìn)行自動的貝葉斯推理，得出未知參數(shù)的概率分布。PyMC是一個馬爾科夫鏈蒙特卡洛采樣(Markov chain Monte Carlo Sampling)工具包，是一個做貝葉斯分析使用的Python庫,包含了畫圖，擬合優(yōu)化度和聚類診斷的方法。本文將利用PyMC來進(jìn)行貝葉斯的推斷[1]。

二、實(shí)驗(yàn)方法

貝葉斯是通過引入先驗(yàn)的不確定性，我們事實(shí)上允許了我們的主觀判斷可能是錯誤的。在觀察數(shù)據(jù)，證據(jù)或其他信息后，我們不斷的更新我們的判斷使得它錯的不那么離譜。所以在貝葉斯推斷的過程中，最重要的就是先驗(yàn)分布的構(gòu)建和確立。第一步，先確定需要推斷的對象以及和對象相關(guān)的參數(shù)。例如我們要推斷一個西瓜是否是好瓜，那么我們需要推斷的對象就是西瓜的好壞，而和西瓜好壞相關(guān)的參數(shù)可以是西瓜的色澤，根蒂，敲聲等等[2]。第二步，確定參數(shù)和對象的模型。根據(jù)對象以及參數(shù)的類型以及取值范圍構(gòu)建合理的模型。例如對象的取值為0或者1，那么可以猜測參數(shù)和對象之間的模型為邏輯模型。第三步，確定先驗(yàn)分布的模型。根據(jù)對象的數(shù)據(jù)類型我們需要確定一個先驗(yàn)分布，例如如果數(shù)據(jù)是離散型的，那么對象的分布可能符合Poisson分布；如果數(shù)據(jù)是連續(xù)的，那么對象的分布可能符合指數(shù)分布。

三、實(shí)驗(yàn)過程

(一)實(shí)驗(yàn)背景

1986年1月28號，挑戰(zhàn)者號起飛不久后一個火箭推動器發(fā)生了爆炸。事故的起因是因?yàn)檫B接在火箭推進(jìn)器上的O型圈有缺陷，這種缺陷來自于設(shè)計(jì)的不合理，這種設(shè)計(jì)使得O型圈對很多因素包括外界溫度都是非常敏感的。我們將利用挑戰(zhàn)者號的數(shù)據(jù)來分析最后一次飛行發(fā)生事故的概率,即當(dāng)溫度為31度時，事故發(fā)生的概率。

(二)數(shù)據(jù)整理

數(shù)據(jù)中包括了日期，溫度以及事故是否發(fā)生。由于數(shù)據(jù)中包含了缺失值，可能會影響推斷的結(jié)果，所以我們需要剔除這些包含缺失值的樣本，而且根據(jù)背景我們知道事故的發(fā)生和溫度有關(guān)，而和日期的關(guān)系不太大，所以我們也需要剔除那些對結(jié)果沒什么影響的因素。

(三)構(gòu)建模型

1.數(shù)據(jù)觀察

由于直觀數(shù)據(jù)很難展現(xiàn)溫度和事故發(fā)生的關(guān)系，所以我們首先作圖來觀察它們之間的關(guān)系。因?yàn)槭鹿拾l(fā)生的情況只有兩種，0(不發(fā)生)和1(發(fā)生)。我們可以清晰地看出，隨著外界溫度的下降，發(fā)生事故的概率變得更高。因?yàn)闇囟群褪鹿拾l(fā)生之間沒有一個嚴(yán)格的轉(zhuǎn)折點(diǎn)，所以我們需要對事故發(fā)生的概率p建模。

2.構(gòu)建先驗(yàn)分布

在這個模型模型中，t為溫度，p(t)為事故發(fā)生的概率，β是個不確定變量。

(2)確定β和α的模型。雖然我們確定了溫度和事故發(fā)生概率的模型，但是我們引進(jìn)了兩個輔助參數(shù)，然而我們對這兩個參數(shù)沒有任何的信息，而且可以看出β和α取值可為正可為負(fù)，一般這種情況我們可以想到的正態(tài)分布來模擬，因?yàn)槭篱g萬物很多都滿足正態(tài)分布。

(3)確定先驗(yàn)分布模型。我們前面分析了兩個模型，事故發(fā)生的概率和溫度的模型，事故發(fā)生的概率中參數(shù)β和α的模型，相當(dāng)于如果我們確定了β和α的值，我們只能算出一個概率，但是最后的推斷結(jié)果為0或者1。所以我們需要構(gòu)建概率和推斷結(jié)果的模型。這里可以很容易想到伯努利分布。

所以我們最終確定的先驗(yàn)分布模型如下：

Di～Ber(p(ti)),i=1,2,…,N

其中，

α～N(μ,τ)

β～N(μ,τ)

3.利用PyMC3計(jì)算后驗(yàn)分布

PyMC3可以很方便地幫我們算出模型的后驗(yàn)分布，只需要我們確定參數(shù)以及模型即可。

(四)結(jié)果分析

根據(jù)MCMC采樣結(jié)果，下面我們作出了事故發(fā)生的概率關(guān)于溫度的曲線圖。

上面的圖描繪了事故發(fā)生概率的后驗(yàn)分布圖，由于受到溫度，α，β的影響，我們不能從圖中判斷在哪個溫度時概率發(fā)生明顯的變化。但是從圖中的曲線斜率的變化趨勢也可以看出在溫度在[60,70]區(qū)間變化最快，所以我們猜測溫度的影響值在[60,70]之間。我們作出當(dāng)溫度t=31時事故發(fā)生的后驗(yàn)分布如下。

從上面的直方圖可以看出當(dāng)溫度t=31°時事故發(fā)生的概率基本為99.9%以上，說明1986年1月28號O型圈發(fā)生事故幾乎是必然的了。除此之外，我們還可以預(yù)測其他溫度發(fā)生事故的概率。

[1]貝葉斯方法.概率編程與貝葉斯推斷，Cameron Davidson-Pilon

[2]周志華.機(jī)器學(xué)習(xí)

林鎮(zhèn)溪，男，漢族，廣東廉江，本科，華南農(nóng)業(yè)大學(xué)，研究方向機(jī)器學(xué)習(xí)。