本節書摘來自異步社群《貝葉斯方法:機率程式設計與貝葉斯推斷》一書中的第1章,第1.4節,作者 【加】cameron davidson-pilon(卡梅隆 戴維森-皮隆),更多章節内容可以通路雲栖社群“異步社群”公衆号檢視
接下來模拟一個有趣的執行個體,它是一個有關使用者發送和收到短信的例子。
你得到了系統中一個使用者每天的短信條數資料,如圖1.4.1中所示。你很好奇這個使用者的短信使用行為是否随着時間有所改變,不管是循序漸進地還是突然地變化。怎麼模拟呢?(這實際上是我自己的短信資料。盡情地判斷我的受歡迎程度吧。)

在模組化之前,僅僅從圖1.4.1中你能猜到什麼嗎?你能說在這一段時間内使用者行為出現了什麼變化嗎?
我們怎麼模拟呢?像前文提到的那樣,possion随機變量能很好地模拟這種計數類型的資料。用ci表示第i天的短信條數。
我們不能确定參數λ的真實取值,然而,在圖1.4.1中,整個觀察周期的後期收到短信的幾率升高了,也可以說,λ在某些時段增加了(在前文中有提到過,當λ取大值的時候更容易得到較大的結果值。在這裡,也就是說一天收到短信條數比較多的機率增大了)。
怎麼用資料表示這種觀察呢?假設在觀察期的某些天(稱它為τ),參數λ的取值突然變得比較大。是以參數λ存在兩個取值:在時段τ之前有一個,在其他時段有另外一個。在一些資料中,像這樣的一個轉變稱之為轉換點:
如果實際上不存在這樣的情況,即λ1=λ2,那麼λ的先驗分布應該是均等的。
對于這些不知道的λ我們充滿了興趣。在貝葉斯推斷下,我們需要對不同可能值的λ配置設定相應的先驗機率。對參數λ1和λ2來說什麼才是一個好的先驗機率分布呢?前面提到過λ可以取任意正數。像我們前面見到的那樣,指數分布對任意正數都存在一個連續密度函數,這或許對模拟λi來說是一個很好的選擇。但也像前文提到的那樣,指數分布也有它自己對應的參數,是以這個參數也需要包括在我們的模型裡面。稱它為參數α。
α被稱為超參數或者是父變量。按字面意思了解,它是一個對其他參數有影響的參數。按照我們最初的設想,α應該對模型的影響不會太大,是以可以對它進行一些靈活的設定。在我們的模型中,我們不希望對這個參數賦予太多的主觀色彩。但這裡建議将其設定為樣本中計算平均值的逆。為什麼這樣做呢?既然我們用指數分布模拟參數λ,那這裡就可以使用指數函數的期望值公式得到:
使用這個值,我們是比較客觀的,這樣做的話可以減少超參數對模拟造成的影響。另外,我也非常建議大家對上面提到的不同時段的兩個λi使用不同的先驗。建構兩個不同α值的指數分布反映出我們的先驗估計,即在整個觀測過程中,收到短信的條數出現了變化。
對于參數τ,由于受到噪聲資料的影響,很難為它挑選合适的先驗。我們假定每一天的先驗估計都是一緻的。用公式表達如下:
做了這麼多了,那麼未知變量的整體先驗分布情況應該是什麼樣的呢?老實說,這并不重要。我們要知道的是,它并不好看,包括很多隻有數學家才會喜歡的符号。而且我們的模型會是以變得更加複雜。不管這些啦,畢竟我們關注的隻是先驗分布而已。
下面會介紹pymc,它是一個由數學奇才們建立起來的關于貝葉斯分析的python庫。
pymc是一個做貝葉斯分析使用的python庫。它是一個運作速度快、維護得很好的庫。它唯一的缺點是,它的說明文檔在某些領域有所缺失,尤其是在一些能區分菜鳥和高手的領域。本書的主要目标就是解決問題,并展示pymc庫有多酷。
下面用pymc模拟上面的問題。這種類型的程式設計被稱之為機率程式設計,對此的誤稱包括随機産生代碼,而且這個名字容易使得使用者誤解或者讓他們對這個領域産生恐懼。代碼當然不是随機的,之是以名字裡面包含機率是因為使用編譯變量作為模型的元件建立了機率模型。在pymc中模型元件為第一類原語。
cronin對機率程式設計有一段激動人心的描述:
“換一種方式考慮這件事情,跟傳統的程式設計僅僅向前運作不同的是,機率程式設計既可以向前也可以向後運作。它通過向前運作來計算其中包含的所有關于世界的假設結果(例如,它對模型空間的描述),但它也從資料中向後運作,以限制可能的解釋。在實踐中,許多機率程式設計系統将這些向前和向後的操作巧妙地交織在一起,以給出有效的最佳的解釋。
由于“機率程式設計”一詞會産生很多不必要的困惑,我會克制自己使用它。相反,我會簡單地使用“程式設計”,因為它實際上就是程式設計。
pymc代碼很容易閱讀。唯一的新東西應該是文法,我會截取代碼來解釋各個部分。隻要記住我們模型的元件(τ,λ1,λ2)為變量:
在這段代碼中,我們産生了對應于參數λ1和λ2的pymc變量,并令他們為pymc中的随機變量,之是以這樣稱呼它們是因為它們是由背景的随機數産生器生成的。為了表示這個過程,我們稱它們由random()方法建構。在整個訓練階段,我們會發現更好的tau值。
這段代碼産生了一個新的函數lambda_,但事實上我們可以把它想象成為一個随機變量——之前的随機參數λ。注意,因為lambda_1、lambda_2、tau是随機的,是以lambda_也會是随機的。目前我們還沒有計算出任何變量。
@pm.deterministic是一個辨別符,它可以告訴pymc這是一個确定性函數,即如果函數的輸入為确定的話(當然這裡它們不是),那函數的結果也是确定的。
model = pm.model([observation, lambda_1, lambda_2, tau])
變量observation包含我們的資料count_data,它是由變量lambda_用我們的資料産生機制得到。我們将observed設定為true來告訴pymc這在我們的分析中是一個定值。最後,pymc 希望我們收集所有變量資訊并從中産生一個model執行個體。當我們拿到結果時就會比較好處理了。
下面的代碼将在第3章中解釋,但在這裡我們展示我們的結果是從哪裡來的。可以把它想象成為一個不斷學習的過程。這裡使用的理論稱為馬爾科夫鍊蒙特卡洛(mcmc),在第3章中會給出進一步的解釋。利用它可以得到參數λ1、λ2和τ先驗分布中随機變量的門檻值。我們對這些随機變量作直方圖,觀測他們的先驗分布。接下來,将樣本(在mcmc中我們稱之為迹)放入直方圖中。結果如圖1.4.2所示。
回想一下,貝葉斯方法傳回一個分布。是以,我們現在有分布描述未知的λ和τ。我們得到了什麼?馬上,我們可以看我們估計的不确定性:分布越廣,我們的後驗機率越小。我們也可以看到參數的合理值:λ1大概為18,λ2大概是23。這兩個λ的後驗分布明顯是不同的,這表明使用者接收短信的行為确實發生了變化。(請參閱1.6補充說明中的正式參數。)
你還能做哪些其他的觀測呢?再看看原始資料,你是否覺得這些結果合理呢?
還要注意到λ的後驗分布并不像是指數分布,事實上,後驗分布并不是我們從原始模型中可以辨識的任何分布。但這挺好的!這是用計算機處理的一個好處。如果我們不這麼做而改用數學方式處理這個問題,将會非常的棘手和混亂。使用計算數學的方式可以讓我們不用在友善數學處理上考慮太多。
我們的分析頁傳回了τ的一個分布。由于它是一個離散變量,是以它的後驗看起來和其他兩個參數有點不同,它不存在機率區間。我們可以看到,在45天左右,有50%的把握可以說使用者的行為是有所改變的。沒有發生變化,或者随着時間有了慢慢的變化,τ的後驗分布會更加的分散,這也反映出在很多天中τ是不太好确定的。相比之下,從真實的結果中看到,僅僅有三到四天可以認為是潛在的轉折點。
在這本書的其餘部分,我們會面對這樣一個問題,它是一個能帶領我們得到強大結果的說明。現在,用另外一個執行個體結束這一章。
我們會用後驗樣本回答下面的問題:在第t(0≤t≤70)天中,期望收到多少條短信呢?poisson分布的期望值等于它的參數λ。是以問題相當于:在時間t中,參數λ的期望值是多少。
在下面的代碼中,令i訓示後驗分布中的變量。給定某天t,我們對所有可能的λ求均值,如果 t < τi (也就是說,如果并沒有發生什麼變化),令λi=λ1,i,否則我們令λi = λ2,i。
在圖1.4.3中展示的結果,很明顯地說明了轉折點的重要影響。但是對此我們應該保持謹慎的态度,從這個觀點看,這裡并不存在我們非常希望看到的“短信的期望數量”這樣的一條直線。我們的分析結果非常符合之前的估計——使用者行為确實發生了改變(如不然,λ1和λ2的值應該比較接近),而且這是一個突然的變化,而非一種循序漸進的變化(τ的先驗分布突然出現了峰值)。我們可以推測這種情況産生的原因是:短信費用的降低,天氣提醒短信的訂閱,或者是一段新的感情。