程序員如何用貝葉斯方式思考（貝葉斯算法的實現步驟）

1 引言

盡管你已是一個編程老手，但bug仍有可能在代碼中存在。于是，在實現了一段特別難的算法之后，你決定先來一個簡單的測試用例。這個用例通過了。接著你用了一個稍微復雜的測試用例。再次通過了。接下來更難的測試用例也通過了。這時，你開始覺得也許這段代碼已經沒有bug了。

如果你這樣想，那么恭喜你：你已經在用貝葉斯的方式思考！簡單地說，貝葉斯推斷是通過新得到的證據不斷地更新你的信念。貝葉斯推斷很少會做出絕對的判斷，但可以做出非?？尚诺呐袛?。在上面的例子中，我們永遠無法100%肯定我們的代碼是無缺陷的，除非我們測試每一種可能出現的情形，這在實踐中幾乎不可能。但是，我們可以對代碼進行大量的測試，如果每一次測試都通過了，我們更有把握覺得這段代碼是沒問題的。貝葉斯推斷的工作方式就在這里：我們會隨著新的證據不斷更新之前的信念，但很少做出絕對的判斷，除非所有其他的可能都被一一排除。

1.1　貝葉斯思維

和更傳統的統計推斷不同，貝葉斯推斷會保留不確定性。乍聽起來，這像一門糟糕的統計方法，難道不是所有的統計都是期望從隨機性里推斷出確定性嗎？要協調這些不一致，我們首先需要像貝葉斯派一樣思考。 在貝葉斯派的世界觀中，概率是被解釋為我們對一件事情發生的相信程度，換句話說，這表明了我們對此事件發生的信心。事實上，我們一會兒就會發現，這就是概率的自然的解釋。 為了更清晰地論述，讓我們看看來自頻率派關于概率的另一種解釋。頻率派是更古典的統計學派，他們認為概率是事件在長時間內發生的頻率。例如，在頻率派的哲學語境里，飛機事故的概率指的是長期來看，飛機事故的頻率值。對許多事件來說，這樣解釋概率是有邏輯的，但對某些沒有長期頻率的事件來說，這樣解釋是難以理解的。試想一下：在總統選舉時，我們經常提及某某候選人獲選的概率，但選舉本身只發生一次！頻率論為了解決這個問題，提出了“替代現實”的說法，從而用在所有這些替代的“現實”中發生的頻率定義了這個概率。

貝葉斯派，在另一方面，有更直觀的方式。它把概率解釋成是對事件發生的信心。簡單地說，概率是觀點的概述。某人把概率0賦給某個事件的發生，表明他完全確定此事不會發生；相反，如果賦的概率值是1，則表明他十分肯定此事一定會發生。0和1之間的概率值可以表明此事可能發生的權重。這個概率定義可以和飛機事故概率一致。如果除去所有外部信息，一個人對飛機事故發生的信心應該等同于他了解到的飛機事故的頻率。同樣，貝葉斯概率的定義也能適用于總統選舉，并且使得概率（信心）更加有意義：你對候選人A獲勝的信心有多少？

請注意，在之前，我們提到每個人都可以給事件賦概率值，而不是存在某個唯一的概率值。這很有趣，因為這個定義為個人之間的差異留有余地。這正和現實天然契合：不同的人即便對同一事件發生的信心也可以有不同的值，因為他們擁有不同的信息。這些不同并不說明其他人是錯誤的。考慮下面的例子。

1．在拋硬幣中我們同時猜測結果。假設硬幣沒有被做手腳，我和你應該都相信正反面的概率都是0.5。但假設我偷看了這個結果，不管是正面還是反面，我都賦予這個結果1的概率值?，F在，你認為正面的概率是多少？很顯然，我額外的信息（偷看）并不能改變硬幣的結果，但使得我們對結果賦予的概率值不同了。

2．你的代碼中也許有一個bug，但我們都無法確定它的存在，盡管對于它是否存在我們有不同的信心。

3．一個病人表現出x、y、z三種癥狀，有很多疾病都會表現出這三種癥狀，但病人只患了一種病。不同的醫生對于到底是何種疾病導致了這些癥狀可能會有稍微不同的看法。

對我們而言，將對一個事件發生的信心等同于概率是十分自然的，這正是我們長期以來同世界打交道的方式。我們只能了解到部分的真相，但可以通過不斷收集證據來完善我們對事物的觀念。與此相對的是，你需要通過訓練才能夠以頻率派的方式思考事物。

為了和傳統的概率術語對齊，我們把對一個事件A發生的信念記為P(A)，這個值我們稱為先驗概率。

偉大的經濟學家和思想者John Maynard Keynes曾經說過（也有可能是杜撰的）：“當事實改變，我的觀念也跟著改變，你呢？”這句名言反映了貝葉斯派思考事物的方式，即隨著證據而更新信念。甚至是，即便證據和初始的信念相反，也不能忽視了證據。我們用P(A|X)表示更新后的信念，其含義為在得到證據X后，A事件的概率。為了和先驗概率相對，我們稱更新后的信念為后驗概率?？紤]在觀察到證據X后，以下例子的后驗概率。

1．P(A)：硬幣有50%的幾率是正面。P(A|X)：你觀察到硬幣落地后是正面，把這個觀察到的信息記為X，那么現在你會賦100%的概率給正面，0%的概率給負面。

2．P(A)：一段很長很復雜的代碼可能含有bug。P(A|X):代碼通過了所有的X個測試；現在代碼可能仍然有bug，不過這個概率現在變得非常小了。

3．P(A)：病人可能有多種疾病。P(A|X)：做了血液測試之后，得到了證據X，排除了之前可能的一些疾病。

在上述例子中，即便獲得了新的證據，我們也并沒有完全地放棄初始的信念，但我們重新調整了信念使之更符合目前的證據（也就是說，證據讓我們對某些結果更有信心）。

通過引入先驗的不確定性，我們事實上允許了我們的初始信念可能是錯誤的。在觀察數據、證據或其他信息之后，我們不斷更新我們的信念使得它錯得不那么離譜。這和硬幣預測正相反，我們通常希望預測得更準確。

1.2 貝葉斯推斷在實踐中的運用

如果頻率推斷和貝葉斯推斷是一種編程函數，輸入是各種統計問題，那么這兩個函數返回的結果可能是不同的。頻率推斷返回一個估計值（通常是統計量，平均值是一個典型的例子），而貝葉斯推斷則會返回概率值。

例如，在代碼測試的例子中，如果你問頻率函數：“我的代碼通過了所有測試，它現在沒有bug了嗎？”頻率函數會給出“yes”的回答。但如果你問貝葉斯函數：“通常我的代碼有bug，現在我的代碼通過了所有測試，它是不是沒有bug了？”貝葉斯函數會給出非常不同的回答，它會給出“yes”和“no”的概率，例如“‘yes’的概率是80%，‘no’的概率是20%?！?/p>

這和頻率函數返回的結果是非常不同的。注意到貝葉斯函數還有一個額外的信息——“通常的我的代碼有bug”，這個參數就是先驗信念。把這個參數加進去，貝葉斯函數會將我們的先驗概率納入考慮范圍。通常這個參數是可省的，但我們將會發現缺省它會產生什么樣的結果。

加入證據　當我們添加更多的證據，初始的信念會不斷地被“洗刷”。這是符合期望的，例如如果你的先驗是非?；闹嚨男拍铑愃啤疤柦裉鞎ā保敲茨忝恳惶於紩淮蚰?，這時候你希望某種統計推斷的方法會糾正初始的信念，或者至少讓初始的信念變得不那么荒謬。貝葉斯推斷就是你想要的。

讓N表示我們擁有的證據的數量。如果N趨于無窮大，那么貝葉斯的結果通常和頻率派的結果一致。因此，對于足夠大的N，統計推斷多多少少都還是比較客觀的。另一方面，對于較小的N，統計推斷相對而言不穩定，而頻率派的結果有更大的方差和置信區間。貝葉斯在這方面要勝出了。通過引入先驗概率和返回概率結果（而不是某個固定值），我們保留了不確定性，這種不確定性正是小數據集的不穩定性的反映。

有一種觀點認為，對于大的N來說，兩種統計技術是無差別的，因為結果類似，并且頻率派的計算要更為簡單，因而傾向于頻率派的方法。在采納這種觀點之前，也許應該先參考Andrew Gelman的論述：

“樣本從來都不是足夠大的。如果N太小不足以進行足夠精確的估計，你需要獲得更多的數據。但當N“足夠大”，你可以開始通過劃分數據研究更多的問題（例如在民意調查中，當你已經對全國的民意有了較好的估計，你可以開始分性別、地域、年齡進行更多的統計）。N從來無法做到足夠大，因為當它一旦大了，你總是可以開始研究下一個問題從而需要更多的數據?！?（參見http://andrewgelman.com/2005/07/31/n_is_never_largl.）

1.3　頻率派的模型是錯誤的嗎？

不。頻率方法仍然非常有用，在很多領域可能都是最好的辦法。例如最小方差線性回歸、LASSO回歸、EM算法，這些都是非常有效和快速的方法。貝葉斯方法能夠在這些方法失效的場景發揮作用，或者是作為更有彈性的模型而補充。

1.4　關于大數據

出乎意料的是，通常解決大數據預測型問題的方法都是些相對簡單的算法。因此，我們會認為大數據預測的難點不在于算法，而在于大規模數據的存儲和計算。（讓我們再次回顧Gelman的關于樣本大小的名言，并且提問：“我真的有大數據嗎？”）

中等規?；蛘吒∫幠５臄祿沟梅治鲎兊酶鼮槔щy。用類似Gelman的話說，如果大數據問題能夠被大數據方法簡單直接地解決，那么我們應該更關注不那么大的數據集上的問題。

2 我們的貝葉斯框架

我們感興趣的估計，可以通過貝葉斯的思想被解釋為概率。我們對事件A有一個先驗估計——例如，在準備測試之前，我們對代碼中的漏洞就有了一個先驗的估計。

接下來，觀察我們的證據。繼續拿代碼漏洞為例：如果我們的代碼通過了X個測試，我們會相應地調整心里的估計。我們稱這個調整過后的新估計為后驗概率。調整這個估計值可以通過下面的公式完成，這個公式被稱為貝葉斯定理，得名于它的創立者托馬斯·貝葉斯。

∝P(X|A)P(A) (∝ 的意思是“與之成比例”)

上面的公式并不等同于貝葉斯推論，它是一個存在于貝葉斯推論之外的數學真理。在貝葉斯推論里它僅僅被用來連接先驗概率P(A)和更新的后驗概率P(A|X)。

2.1　不得不講的實例：拋硬幣

幾乎所有統計書籍都包含一個拋硬幣的實例，那我也從這個開始著手吧。假設你不確定在一次拋硬幣中得到正面的概率（劇透警告：它是50%），你認為這里肯定是存在某個比例的，稱之為p，但是你事先并不清楚p大概會是多少。

我們開始拋硬幣，并記錄下每一次拋出的結果——正面或反面，這就是我們的觀測數據。一個有趣的問題是：“隨著收集到越來越多的數據，我們對p的推測是怎么變化的呢？”

說得更具體一些，當面對著很少量的數據或擁有大量數據時，我們的后驗概率是怎么樣的呢？下面，我們按照觀測到的越來越多的數據（拋硬幣數據），逐次更新我們的后驗概率圖。

在圖中我們用曲線表示我們的后驗概率，曲線越寬，我們的不確定性越大。如圖2.1所示，當我們剛剛開始觀測的時候，我們的后驗概率的變化是不穩定的。但是最終，隨著觀測數據（拋硬幣數據）越來越多，這個概率會越來越接近它的真實值p=0.5（圖中用虛線標出）。

注意到圖中的波峰不一定都出現在0.5那里，當然它也沒有必要都這樣。應該明白的是我們事前并不知道p會是多少。事實上，如果我們的觀測十分的極端，比如說拋了8次只有1次結果是正面的，這種情況我們的分布會離0.5偏差很多（如果缺少先驗的知識，當出現8次反面1次正面時，你真的會認為拋硬幣結果是公平的嗎？）。隨著數據的累積，我們可以觀察到，雖然不是每個時候都這樣，但越來越多地，概率值會出現在p=0.5。

下面這個實例就簡單地從數據角度演示一下貝葉斯推斷。

圖2.1　后驗概率的貝葉斯變換

2.2　實例：圖書管理員還是農民

下面這個故事靈感來自于Daniel Kahneman的《思考，快與慢》一書，史蒂文被描述為一個害羞的人，他樂于助人，但是他對其他人不太關注。他非常樂見事情處于合理的順序，并對他的工作非常細心。你會認為史蒂文是一個圖書管理員還是一個農民呢？從上面的描述來看大多數人都會認為史蒂文看上去更像是圖書管理員，但是這里卻忽略了一個關于圖書管理員和農民的事實：男性圖書管理員的人數只有男性農民的1/20。所以從統計學來看史蒂文更有可能是一個農民。

怎么正確地看待這個問題呢？史蒂文實際上更有可能是一個農民還是一個圖書管理員呢？把問題簡化，假設世上只有兩種職業——圖書管理員和農民，并且農民的數量確實是圖書管理員的20倍。

設事件A為史蒂文是一個圖書管理員。如果我們沒有史蒂文的任何信息，那么P（A）=1/21=0.047。這是我們的先驗?，F在假設從史蒂文的鄰居們那里我們獲得了關于他的一些信息，我們稱它們為X。我們想知道的就是P（A|X）。由貝葉斯定理得：

我們知道P（A）是什么意思，那P（X|A）是什么呢？它可以被定義為在史蒂文真的是一個圖書管理員的情況下，史蒂文的鄰居們給出的某種描述的概率，即如果史蒂文真的是一個圖書管理員，他的鄰居們將他描述為一個圖書管理員的概率。這個值很可能接近于1。假設它為0.95。

P（X）可以解釋為：任何人對史蒂文的描述和史蒂文鄰居的描述一致的概率?，F在這種形式有點難以理解，我們將其做一些邏輯上的改造：

P(X)=P(X?and?A)+P(X?and?～A) 　　　　　　　　　　　　　　　　　=P(X|A)P(A)+P(X|～A)P(～A)

其中～A表示史蒂文不是一個圖書管理員的事件，那么他一定是一個農民?，F在我們知道P（X|A）和P(A)，另外也可知P(～A)=1-P(A)=20/21?，F在我們只需要知道P(X|～A)，即在史蒂文為一個農民的情況下，史蒂文的鄰居們給出的某種描述的概率即可。假設它為0.5，這樣，

結合以上：

這個值并不算高，但是考慮到農民的數量比圖書管理員的數量多這么多，這個結果也非常合理了。在圖2.2中，對比了在史蒂文為農民和史蒂文為圖書管理員時的先驗和后驗概率。

%matplotlib?inlinefrom?IPython.core.pylabtools?import?figsizeimport?numpy?as?np from?matplotlib?import?pyplot?as?pltfigsize(12.5,?4)plt.rcParams['savefig.dpi']?=?300plt.rcParams['figure.dpi']?=?300　 colors?=?["#348ABD",?"#A60628"] prior?=?[1/21.,?20/21.] posterior?=?[0.087,1-0.087] plt.bar([0,?.7],?prior,?alpha=0.70,?width=0.25, ???????color=colors[0],?label="prior?distribution", ???????lw="3",?edgecolor="#348ABD") 　 plt.bar([0+0.25,?.7+0.25],?posterior,?alpha=0.7, ???????width=0.25,?color=colors[1], ???????label="posterior?distribution", ???????lw="3",?edgecolor="#A60628") 　 plt.xticks([0.20,?0.95],?["Librarian",?"Farmer"]) plt.title("Prior?and?posterior?probabilities?of?Steve's\ ?????????occupation") plt.ylabel("Probability") plt.legend(loc="upper?left");

在我們得到X的觀測值之后，史蒂文為圖書管理員的概率增加了，雖然增加的不是很多，史蒂文為農民的可能性依舊是相當大的。

這是一個關于貝葉斯推斷和貝葉斯法則的一個簡單的實例。不幸的是，除了在人工結構的情況下，要執行更加復雜的貝葉斯推斷所使用到的數學只會變得更加的復雜。在后面我們將看到執行這種復雜的屬性分析并沒有必要。首先，我們必須擴充我們的建模工具。

圖2.2　史蒂文職業的先驗及后驗概率

3　概率分布

首先定義以下希臘文字的發音：

α = alpha

β = beta

λ = lambda

μ = mu

σ = sigma

τ = tau

很好。接下來正式開始講概率分布。首先快速地回憶一下什么是概率分布。設Z為一個隨機變量，那么就存在一個跟Z相關的概率分布函數，給定Z任何取值，它都得到一個相應的概率值。

我們把隨機變量分為3種類別。

Z為離散的。離散隨機變量的取值只能是在特定的列表中。像貨幣、電影收視率、投票數都是離散隨機變量。當我們將它與連續型隨機變量對比時，這個概念就比較清楚了。

Z為連續的。連續型隨機變量的值可以是任意精確數值。像溫度、速度和時間都是連續型變量，因為對于這些數值你可以精確到任何程度。

Z為混合的?；旌闲碗S機變量的值可以為以上兩種形式，即結合了以上兩種隨機變量的形式。

離散情況

如果Z是離散的，那么它的分布為概率質量函數，它度量的是當Z取值為k時的概率，用P(Z=k)表示?？梢宰⒁獾?，概率質量函數完完全全地描述了隨機變量Z，即如果知道Z的概率質量方程，那么Z會怎么表現都是可知的。下面介紹一些經常會碰到的概率質量方程，學習它們是十分有必要的。第一個要介紹的概率質量方程十分重要，設Z服從于Poisson分布：

λ被稱為此分布的一個參數，它決定了這個分布的形式。對于Poisson分布來說，λ可以為任意正數。隨著λ的增大，得到大值的概率會增大；相反地，當λ減小時，得到小值的概率會增大。λ可以被稱為Poisson分布的強度。

跟λ可以為任意正數不同，值k可以為任意非負整數，即k必須為0、1、2之類的值。這個是十分重要的，因為如果你要模擬人口分布，你是不可以假設有4.25個或是5.612個人的。

如果一個變量Z存在一個Poisson質量分布，我們可以表示為：

Z～Poi(λ)

Poisson分布的一個重要性質是：它的期望值等于它的參數。即：

E[Z|λ]=λ

這條性質以后經常會被用到，所以記住它是很有用的。在圖3.1中，展示了不同λ取值下的概率質量分布。首先值得注意的是，增大λ的取值，k取大值的概率會增加。其次值得注意的是，雖然x軸在15的時候截止，但是分布并沒有截止，它可以延伸到任意非負的整數。

figsize(12.5,?4) 　import?scipy.stats?as?stats a?=?np.arange(16) poi?=?stats.poisson lambda_?=?[1.5,?4.25] colors?=?["#348ABD",?"#A60628"] 　 plt.bar(a,?poi.pmf(a,?lambda_[0]),?color=colors[0], ???????label="$\lambda?=?%.1f$"?%?lambda_[0],?alpha=0.60, ???????edgecolor=colors[0],?lw="3") 　 plt.bar(a,?poi.pmf(a,?lambda_[1]),?color=colors[1], ???????label="$\lambda?=?%.1f$"?%?lambda_[1],?alpha=0.60, ???????edgecolor=colors[1],?lw="3") plt.xticks(a?+?0.4,?a) plt.legend() plt.ylabel("Probability?of?$k$") plt.xlabel("$k$") plt.title("Probability?mass?function?of?a?Poisson?random?variable,\ ?????????differing?$\lambda$?values");

圖3.1　不同λ取值情況下，Poisson隨機變量的概率質量函數

連續情況

對應于離散情況下的概率質量函數，連續情況下概率分布函數被稱為概率密度函數。雖然在命名上作這樣的區分看起來是沒有必要的，但是概率質量函數和概率密度函數有著本質的不同。舉一個連續型隨機變量的例子：指數密度。指數隨機變量的密度函數如下式：

fZ(z|λ)=λe-λz,z≥0

類似于Poisson隨機變量，指數隨機變量只可以取非負值。但是和Poisson分布不同的是，這里的指數可以取任意非負值，包括非整數，例如4.25或是5.612401。這個性質使得計數數據（必須為整數）在這里不太適用；而對于類似時間數據、溫度數據（當然是以開氏溫標計量）或其他任意對精度有要求的正數數據，它是一種不錯的選擇。圖3.2展示了λ取不同值時的概率密度函數。

當隨機變量Z擁有參數為λ的指數分布時，我們稱Z服從于指數分布，并記為：

Z～Exp(λ)

對指定的參數λ，指數型隨機變量的期望值為λ的逆，即

E[Z|λ]=1/λ

a?=?np.linspace(0,?4,?100) expo?=?stats.expon lambda_?=?[0.5,?1] 　for?l,?c?in?zip(lambda_,?colors): ???plt.plot(a,?expo.pdf(a,?scale=1./l),?lw=3, ??????????????color=c,?label="$\lambda?=?%.1f$"?%?l) ???plt.fill_between(a,?expo.pdf(a,?scale=1./l),?color=c,?alpha=.33) 　 plt.legend() plt.ylabel("Probability?density?function?at?$z$") plt.xlabel("$z$") plt.ylim(0,1.2) plt.title("Probability?density?function?of?an?exponential?random\ ?????????variable,?differing?$\lambda$?values");

圖3.2　不同λ取值情況下，指數分布的概率密度函數

這里值得注意的是，概率密度方程在某一點的值并不等于它在這一點的概率。這個將會在后面講到，當然如果你對它感興趣，可以加入下面的討論：http://stats.stackexchange.com/questions/4220/a-probability-distribution-value-exceeding-1-is-ok。

3.3　什么是λ

這個問題可以理解為統計的動機是什么。在現實世界，我們并不知道λ的存在，我們只能直觀地感受到變量Z，如果確定參數λ的話，那就一定要深入到整個事件的背景中去。這個問題其實很難，因為并不存在Z到λ的一一對應關系。對于λ的估計有很多的設計好的方法，但因為λ不是一個可以真實觀察到的東西，誰都不能說哪種方式一定是最好的。

貝葉斯推斷圍繞著對λ取值的估計。與其不斷猜測λ的精確取值，不如用一個概率分布來描述λ的可能取值。

起初這看起來或許有些奇怪。畢竟，λ是一個定值，它不一定是隨機的！我們怎么可能對一個非隨機變量值賦予一個概率呢？不，這樣的考慮是老舊的頻率派思考方式。如果我們把它們理解為估計值的話，在貝葉斯的哲學體系下，它們是可以被賦予概率的。因此對參數λ估計是完全可以接受的。

4　使用計算機執行貝葉斯推斷

接下來模擬一個有趣的實例，它是一個有關用戶發送和收到短信的例子。

4.1　實例：從短信數據推斷行為

你得到了系統中一個用戶每天的短信條數數據，如圖1.4.1中所示。你很好奇這個用戶的短信使用行為是否隨著時間有所改變，不管是循序漸進地還是突然地變化。怎么模擬呢？（這實際上是我自己的短信數據。盡情地判斷我的受歡迎程度吧。）

figsize(12.5,?3.5) count_data?=?np.loadtxt("data/txtdata.csv") n_count_data?=?len(count_data) plt.bar(np.arange(n_count_data),?count_data,?color="#348ABD") plt.xlabel("Time?(days)") plt.ylabel("Text?messages?received") plt.title("Did?the?user's?texting?habits?change?over?time?") plt.xlim(0,?n_count_data);

圖4.1　用戶的短信使用行為是否隨著時間發生改變?

在建模之前，僅僅從圖4.1中你能猜到什么嗎？你能說在這一段時間內用戶行為出現了什么變化嗎？

我們怎么模擬呢？像前文提到的那樣，Possion隨機變量能很好地模擬這種計數類型的數據。用Ci表示第i天的短信條數。

Ci～Poi(λ)

我們不能確定參數λ的真實取值，然而，在圖4.1中，整個觀察周期的后期收到短信的幾率升高了，也可以說，λ在某些時段增加了（當λ取大值的時候更容易得到較大的結果值。在這里，也就是說一天收到短信條數比較多的概率增大了）。

怎么用數據表示這種觀察呢？假設在觀察期的某些天（稱它為τ），參數λ的取值突然變得比較大。所以參數λ存在兩個取值：在時段τ之前有一個，在其他時段有另外一個。在一些資料中，像這樣的一個轉變稱之為轉換點：

如果實際上不存在這樣的情況，即λ1=λ2，那么λ的先驗分布應該是均等的。

對于這些不知道的λ我們充滿了興趣。在貝葉斯推斷下，我們需要對不同可能值的λ分配相應的先驗概率。對參數λ1和λ2來說什么才是一個好的先驗概率分布呢？前面提到過λ可以取任意正數。像我們前面見到的那樣，指數分布對任意正數都存在一個連續密度函數，這或許對模擬λi來說是一個很好的選擇。但也像前文提到的那樣，指數分布也有它自己對應的參數，所以這個參數也需要包括在我們的模型里面。稱它為參數α。

λ1～Exp(α)

λ2～Exp(α)

α被稱為超參數或者是父變量。按字面意思理解，它是一個對其他參數有影響的參數。按照我們最初的設想，α應該對模型的影響不會太大，所以可以對它進行一些靈活的設定。在我們的模型中，我們不希望對這個參數賦予太多的主觀色彩。但這里建議將其設定為樣本中計算平均值的逆。為什么這樣做呢？既然我們用指數分布模擬參數λ，那這里就可以使用指數函數的期望值公式得到：

使用這個值，我們是比較客觀的，這樣做的話可以減少超參數對模擬造成的影響。另外，我也非常建議大家對上面提到的不同時段的兩個λi使用不同的先驗。構建兩個不同α值的指數分布反映出我們的先驗估計，即在整個觀測過程中，收到短信的條數出現了變化。

對于參數τ，由于受到噪聲數據的影響，很難為它挑選合適的先驗。我們假定每一天的先驗估計都是一致的。用公式表達如下：

τ～DiscreteUniform(1,70)

＝﹥P(τ=k)=1/70

做了這么多了，那么未知變量的整體先驗分布情況應該是什么樣的呢？老實說，這并不重要。我們要知道的是，它并不好看，包括很多只有數學家才會喜歡的符號。而且我們的模型會因此變得更加復雜。不管這些啦，畢竟我們關注的只是先驗分布而已。

下面會介紹PyMC，它是一個由數學奇才們建立起來的關于貝葉斯分析的Python庫。

PyMC是一個做貝葉斯分析使用的Python庫。它是一個運行速度快、維護得很好的庫。它唯一的缺點是，它的說明文檔在某些領域有所缺失，尤其是在一些能區分菜鳥和高手的領域。本書的主要目標就是解決問題，并展示PyMC庫有多酷。

下面用PyMC模擬上面的問題。這種類型的編程被稱之為概率編程，對此的誤稱包括隨機產生代碼，而且這個名字容易使得使用者誤解或者讓他們對這個領域產生恐懼。代碼當然不是隨機的，之所以名字里面包含概率是因為使用編譯變量作為模型的組件創建了概率模型。在PyMC中模型組件為第一類原語。

Cronin對概率編程有一段激動人心的描述：

“換一種方式考慮這件事情，跟傳統的編程僅僅向前運行不同的是，概率編程既可以向前也可以向后運行。它通過向前運行來計算其中包含的所有關于世界的假設結果（例如，它對模型空間的描述），但它也從數據中向后運行，以約束可能的解釋。在實踐中，許多概率編程系統將這些向前和向后的操作巧妙地交織在一起，以給出有效的最佳的解釋。

由于“概率編程”一詞會產生很多不必要的困惑，我會克制自己使用它。相反，我會簡單地使用“編程”，因為它實際上就是編程。

PyMC代碼很容易閱讀。唯一的新東西應該是語法，我會截取代碼來解釋各個部分。只要記住我們模型的組件（τ，λ1，λ2）為變量：

import?pymc?as?pm 　 alpha?=?1.0/count_data.mean()#?Recall?that?count_data?is?the ???????????????????????????????????#?variable?that?holds?our?text?counts. lambda_1?=?pm.Exponential("lambda_1",?alpha) lambda_2?=?pm.Exponential("lambda_2",?alpha) 　 tau?=?pm.DiscreteUniform("tau",?lower=0,?upper=n_count_data)

在這段代碼中，我們產生了對應于參數λ1和λ2的PyMC變量，并令他們為PyMC中的隨機變量，之所以這樣稱呼它們是因為它們是由后臺的隨機數產生器生成的。為了表示這個過程，我們稱它們由random（）方法構建。在整個訓練階段，我們會發現更好的tau值。

print?"Random?output:",?tau.random(),?tau.random(),?tau.random() [Output]: 　 Random?output:?53?21?42@pm.deterministic def?lambda_(tau=tau,?lambda_1=lambda_1,?lambda_2=lambda_2): ?????out?=?np.zeros(n_count_data)?#?number?of?data?points ?????out[:tau]?=?lambda_1?#?lambda?before?tau?is?lambda_1 ?????out[tau:]?=?lambda_2?#?lambda?after?(and?including)?tau?is ??????????????????????????#?lambda_2?????return?out

這段代碼產生了一個新的函數lambda，但事實上我們可以把它想象成為一個隨機變量——之前的隨機參數λ。注意，因為lambda_1、lambda_2、tau是隨機的，所以lambda也會是隨機的。目前我們還沒有計算出任何變量。

@pm.deterministic是一個標識符，它可以告訴PyMC這是一個確定性函數，即如果函數的輸入為確定的話（當然這里它們不是），那函數的結果也是確定的。

observation?=?pm.Poisson("obs",?lambda_,?value=count_data,

其中

λ1～Exp(α)

λ2～Exp(α)

λ3～Exp(α)

并且

τ1～DiscreteUniform(1,69)

τ2～DiscreteUniform(τ1,70)

我們把這個模型也編譯成代碼，跟前面的代碼看上去差不多。

lambda_1?=?pm.Exponential("lambda_1",?alpha) lambda_2?=?pm.Exponential("lambda_2",?alpha) lambda_3?=?pm.Exponential("lambda_3",?alpha) 　 tau_1?=?pm.DiscreteUniform("tau_1",?lower=0,?upper=n_count_data-1) tau_2?=?pm.DiscreteUniform("tau_2",?lower=tau_1,?upper=n_count_data) 　 @pm.deterministic def?lambda_(tau_1=tau_1,?tau_2=tau_2, ??????????????lamb-da_1=lambda_1,?lambda_2=lambda_2,?lambda_3?=?lambda_3): ????out?=?np.zeros(n_count_data)?#?number?of?data?points ????out[:tau_1]?=?lambda_1?#?lambda?before?tau?is?lambda_1 ????out[tau_1:tau_2]?=?lambda_2 ????out[tau_2:]?=?lambda_3????#?lambda?after?(and?including)?tau ????????????????????????????????????#?is?lambda_2????return?out observa-tion?=?pm.Poisson("obs",?lambda_,?value=count_data,observed=True) mod-el?=?pm.Model([observation,?lambda_1,?lambda_2,?lambda_3,?tau_1, ??????????????????????tau_2]) mcmc?=?pm.MCMC(model) mcmc.sample(40000,?10000)

圖6.1展示了5個未知數的后驗。我們可以看到模型的轉折點大致在第45天和第47天的時候取得。對此你怎么認為呢？我們的模型是否對數據過擬合呢？

確實，我們都可能對數據中有多少個轉折點抱有疑惑的態度。例如，我就認為一個轉折點好過兩個轉折點，兩個轉折點好過三個轉折點，以此類推。這意味著對于應該有多少個轉折點可以設置一個先驗分布并讓模型自己做決定！在對模型進行調整之后，答案是肯定的，一個轉折點確實比較適合。代碼在本章就不展示了，這里我只是想介紹一種思想：用懷疑數據那樣的眼光審視我們的模型。

圖6.1　擴充后的短信模型中5個未知參數的后驗分布

7習題

1．利用 lambda_1_samples 和 lambda_2_samples，怎么獲得參數λ1和λ2后驗分布的平均值 ?

2．計算短信頻率提升比例的期望值是多少。提示：需要計算?(lambda_2_samples-lambda_1_samples)/lambda_1_samples的均值。注意這個結果和(lambda_2_ samples.mean()-lambda_1_samples.mean())/ lambda_1_samples.mean()計算出來的結果是有很大區別的。

3．在τ小于45的前提下，計算λ1的均值。也就是說，在我們被告知行為的變化發生在第45天之前時，對λ1的期望會是多少? （不需要重復PyMC 那部分，只需要考慮當 tau_samples < 45時所有可能的情況。）

8答案

1．計算后驗的均值（即后驗的期望值），我們只需要用到樣本和a.mean函數。

print?lambda_1_samples.mean() print?lambda_2_samples.mean()

2．給定兩個數a 和 b，相對增長可以由 (a ? b)/b給出。在我們的實例中，我們并不能確定λ1和λ2的值是多少。通過計算

(lambda_2_samples-lambda_1_samples)/lambda_1_samples

我們得到另外一個向量，它表示相對增長的后驗，如圖7.1所示。

relative_increase_samples?=?(lambda_2_samples-lambda_1_samples) ???????????????????????????????????/lambda_1_samples print?relative_increase_samples [Output]: 　 [?0.263?0.263?0.263?0.263?...,?0.1622?0.1898?0.1883?0.1883] figsize(12.5,4) plt.hist(relative_increase_samples,?histtype='stepfilled', ???????????bins=30,?alpha=0.85,?color="#7A68A6",?normed=True, ???????????label='posterior?of?relative?increase') plt.xlabel("Relative?increase") plt.ylabel("Density?of?relative?increase") plt.title("Posterior?of?relative?increase") plt.legend();

為了計算這個均值，需要用到新向量的均值：

print?relative_increase_samples.mean() [Output]: 　0.280845247899

圖7.1　相對增長的后驗

3．如果已知 τ < 45，那么所有樣本都需要考慮到這點：

ix?=?tau_samples?

本文節選自《貝葉斯方法：概率編程與貝葉斯推斷》

軟件開發 人工智能 云計算

版權聲明：本文內容由網絡用戶投稿，版權歸原作者所有，本站不擁有其著作權，亦不承擔相應法律責任。如果您發現本站中有涉嫌抄襲或描述失實的內容，請聯系我們jiasou666@gmail.com 處理，核實后本網站將在24小時內刪除侵權內容。