1. 寫在之前的話
0x1：貝葉斯推斷的思想
我們從一個例子開始我們本文的討論。小明是一個程式設計老手，但是依然堅信bug仍有可能在程式碼中存在。於是，在實現了一段特別難的演算法之後，他開始決定先來一個簡單的測試用例，這個用例通過了。接著，他用了一個稍微複雜的測試用例，再次通過了。接下來更難的測試用例也通過了，這時，小明開始覺得這段程式碼出現bug的可能性大大大大降低了....
上面這段白話文中，已經包含了最質樸的貝葉斯思想了！簡單來說，貝葉斯推斷是通過新得到的證據不斷地更新我們的信念。
貝葉斯推斷很少會做出絕對的判斷，但可以做出非常可信的判斷。在上面的例子中，小明永遠無法100%肯定程式碼是無缺陷的，除非他窮舉測試了每一種可能出現的情形，這在實踐中幾乎是不可能的。但是，我們可以對程式碼進行大量的測試，如果每一次測試都通過了，我們“更”有把握覺得這段程式碼是沒問題的。
貝葉斯推斷的工作方式就在這裡：我們會隨著新的證據不斷更新之前的信念，但很少做出絕對的判斷，除非所有其他的可能都被一一排除。
0x2：貝葉斯思維和頻率派思維的區別
1. 頻率派關於概率的解釋 - 統計
頻率派是更古典的統計學派，它們認為概率是事件在長時間內發生的概率。例如，在頻率派的哲學語境裡，飛機事故的概率指的是長期來看，飛機事故的概率值。
對許多事件來說，這樣解釋是有邏輯的。但是對某些沒有長期頻率的事件來說，這樣的解釋是難以理解的。例如，在總統選舉時，某個候選人的獲選的概率是難以用統計的方式來回答的，因為這個候選人的選舉可能只發生一次。頻率論為了解決這個問題，提出了“替代現實”的說法，從而用在所有這些替代的“現實”中發生的概率定義了這個概率。
2. 貝葉斯派關於概率的解釋 - 資訊
貝葉斯派對概率有更直觀的解釋。它把概率解釋成是對事件發生的信心。簡單地說，概率是觀點的概述。
如果某人把概率 0 賦給某個事件的發生，表明他完全確定此事不會發生；相反，如果某人賦的概率值是 1，則表明他十分肯定此事一定會發生。
0 和 1 之間的概率值可以表明此事可能發生的權重。
這個概率定義可以和飛機事故概率一致，如果除去所有外部資訊，一個人對飛機事故發生的信心應該等同於他了解到的飛機事故的頻率。同樣，貝葉斯概率定義也能適用於總統選舉，並且使得概率（信心）更加有意義：即你對候選人 A 獲勝的資訊有多少？
對我們而言，將對一個事件發生的資訊等同於概率是十分自然的，這正是我們長期以來同世界打交道的方式。我們只能瞭解到部分的真相，但可以通過不斷收集證據來完善我們對事物的觀念。
3. 貝葉斯推斷的先驗概率和後驗概率
1）先驗概率
對於統計學派來說，在進行統計之前，我們對概率是一無所知的。對於統計學派和統計方法來說，所有的概率都必須在完成統計之後（包括最大似然估計）才能得到一個統計概率。
但是對於貝葉斯推斷來說，在進行統計之前，我們對事物已經擁有了一個先驗的判斷，即先驗概率，記為 P(A)。
例如：P(A)：硬幣有50%的概率是正面；P(A)：一段很長的很複雜的程式碼可能含有bug
2）後驗概率
我們用 P(A|X) 表示更新後的信念，其含義為在得到證據 X 後，A 事件的概率。考慮下面的例子：
P(A|X)：我們觀察到硬幣落地後是正面，那麼硬幣是正面的後驗概率會比50%大；
P(A|X)：程式碼通過了所有的 X 個測試，現在程式碼可能仍然有bug，但是現在這個概率變得非常小了。
3）貝葉斯推斷柔和地融合了我們對事物的先驗領域知識以及實際的觀察樣本的統計結果
在上述的例子中，即便獲得了新的證據，我們也並沒有完全地放棄初始的信念，但是我們重新調整了信念使之符合目前的證據（訓練樣本），也就是說，證據讓我們對某些結果更有資訊。
通過引入先驗的不確定性，我們事實上允許了我們的初始信念可能是錯誤的。在觀察資料、證據或其他資訊之後，我們不斷更新我們的信念使得它錯得不那麼離譜
0x3：在大資料情況下統計思想和貝葉斯推斷思想是相等的？
頻率推斷和貝葉斯推斷都是一種程式設計函式，輸入的是各種統計問題。頻率推斷返回的是一個統計值（通常是統計量，平均值是一個典型的例子），而貝葉斯推斷則會返回一個概率值。
這個小節我們來討論下，當樣本數量 N 不斷增大時，頻率推斷和貝葉斯推斷會發生什麼。
1. 頻率推斷：N越大，頻率統計越接近真實概率
頻率統計學派假設世界上存在一個真理，即概率。統計做的事只是通過大量的實驗來揭示這個真理本身而已。當 N 趨於無窮時，頻率統計得到的統計值會無限接近於概率本身，即規律本身。
2. 貝葉斯推斷：N越大，先驗的離譜程度會不斷降低，即越靠近真實的規律本身
當我們增加 N 時，即增加更多的證據時，初始的信念會不斷地被“洗刷”。這是符合期望的。例如如果你的先驗是非常荒謬的信念類似“太陽今天會爆炸”，那麼每一天的結果都會不斷“洗刷”這個荒謬的先驗，使得後驗概率不斷下降。
3. N趨於無窮大時，頻率推斷和貝葉斯推斷趨近相等
讓 N 表示我們擁有的證據的數量。如果 N 趨於無窮大，那麼貝葉斯的結果通常和頻率派的結果一致。
因此，對於足夠大的 N，統計推斷多多少少還是比較客觀的。
另一方面，對於較小的 N，統計推斷相對而言不太穩定，頻率統計的結果有更大的方差和置信區間。所以才會有相關性檢驗、卡方檢驗、皮森相關係數這些評價演算法。
貝葉斯推斷在小資料這方面效果要更好，通過引入先驗概率和返回概率結果（而不是某個固定的值），貝葉斯推斷保留了不確定性，這種不確定性正式小資料集的不穩定性的反映。
4. 頻率統計仍然是十分有用的！
雖然頻率統計的有更大的方差和置信區間，在小資料集下表現也不一定穩定。但是頻率統計仍然非常有用，在很多領域可能都是最好的辦法。例如最小方差線性迴歸、LASSO迴歸、EM演算法，這些都是非常有效和快速的方法。
貝葉斯方法能夠在這些方法失效的場景發揮作用，或者是作為更有彈性的模型而補充。
0x4：貝葉斯推斷和極大似然估計、SGD隨機梯度下降是什麼關係？
貝葉斯推斷和極大似然估計、SGD隨機梯度下降一樣，都是一種求解演算法模型超引數（模型分佈）的方法。
貝葉斯推斷更側重於先驗和證據的動態融合，得到極大後驗概率估計；
極大似然估計更側重於對訓練樣本的頻率統計；
而SGD在各種場景下都能表現出較好的效能，通過反向梯度的思想，不斷進行區域性最優的求解，從而獲取近似的全域性最優解。
2. 貝葉斯框架
我們知道，對於我們感興趣的估計，可以通過貝葉斯思想被解釋為概率。在頻率統計中，這個估計時通過統計得到的，而在貝葉斯估計中，估計時通過後驗概率得到的。
我們已經從抽象概念上了解了何謂後驗概率，後驗概率即為：新證據對初始信念洗刷後的結果，即修正後的信念。這個概念怎麼落實到數學計算上呢？
0x1：貝葉斯定理
我們需要一個方法能夠將幾個因素連線起來，即先驗概率、我們的證據、後驗概率。這個方法就是貝葉斯定理，得名於它的創立者托馬斯.貝葉斯。
上面的公式並不等同於貝葉斯推斷，它是一個存在於貝葉斯推斷之外的數學真理。或者應該說貝葉斯推斷是在貝葉斯定理的基石上建立起來的。
在貝葉斯推斷裡，貝葉斯定理公式僅僅被用來連線先驗概率 P(A) 和更新後的後驗概率 P(A|B)
0x2：一個貝葉斯推斷的例子 - 圖書管理員還是農民
斯蒂文是一個害羞的人，他樂於助人，但是他對其他人不太關注。他非常樂見事情處於合理的順序，並對他的工作非常細心。從特徵工程的角度來看，似乎我們更應該判定史蒂文為已個圖書管理員。
但是這裡忽略了一個非常重要的先驗知識，即男性圖書管理員的人數只有男性農民的 1/20，這個先驗知識將極大改變史蒂文是否是圖書管理員這個事件的後驗概率。
1. 形式化定義這個問題
把問題簡化，假設世上只有兩種職業，圖書管理員和農民，並且農民的數量確實是圖書管理員的20倍。
設事件 A 為史蒂文是一個圖書管理員。如果我們沒有史蒂文的任何資訊，那麼 P(A) = 1/21。這是我們的先驗。
我們假設從史蒂文的鄰居們那裡我們獲得了關於他的一些資訊，我們稱它們為 X，我們想知道的就是 P(A|X)。由貝葉斯定理得：P(A|X) = P(X|A) * P(A) / P(X)。
P(X|A)被定義為在史蒂文真的是一個圖書管理員的情況下，史蒂文的鄰居們給出的某種描述的概率，即如果史蒂文真的是一個圖書管理員，他的鄰居們將他描述為一個圖書管理員的概率。這個值很可能接近於1，假設它為 0.95。
P(X) 可以解釋為：任何人史蒂文的描述和史蒂文鄰居的描述一致的概率，我們將其做邏輯上的改造：
P(X) = P(X and A) + P(X and ~A) = P(X|A)P(A) + P(X|~A)P(~A)
P(X|A)P(A) 這兩個值我們都已經知道了。另外，~A 表示史蒂文不是一個圖書管理員的概率，那麼他一定是一個農民。所以 P(~A) = 1 - P(A) = 20 / 21。
P(X|~A) 即在史蒂文是一個農名的情況下，史蒂文的鄰居給出的某個描述的概率，我們假設其為 0.5。
綜上，P(X) = 0.95 * 1/24 + 0.5 * 20/21 = 0.515773809524
帶入貝葉斯定理公式，P(A|X) = 0.95 * 1/20 / 0.515773809524 = 0.0877091748413
可以看到，因為農民數量遠大於圖書管理員，所以後驗概率並不算高，儘管我們看到 P(X|A)的概率要大於 P(X|~A)，即我們得到的特徵描述更傾向於表明史蒂文是一個圖書管理員。
這個例子非常好的闡述了，先驗概率以及證據時如何動態影響最終的後驗概率結果的。
2. 通過視覺化的程式碼呈現該問題
# -*- coding: utf-8 -*- import numpy as np from matplotlib import pyplot as plt plt.rcParams['savefig.dpi'] = 300 plt.rcParams['figure.dpi'] = 300 colours = ["#348ABD", "#A60628"] prior = [1/20., 20/21.] # 圖書管理員和農民的先驗概率分別是 1/20 和 20/21 posterior = [0.087, 1-0.087] # 後驗概率 plt.bar([0, .7], prior, alpha=0.70, width=0.25, color=colours[0], label="prior distribution", lw="3", edgecolor=colours[0]) plt.bar([0 + 0.25, .7 + 0.25], posterior, alpha=0.7, width=0.25, color=colours[1], label="posterior distribution", lw="3", edgecolor=colours[1]) plt.ylim(0, 1) plt.xticks([0.20, 0.95], ['Librarian', 'Farmer']) plt.title("Prior and Posterior probability of Steve's occupation") plt.ylabel("Probability") plt.legend(loc="upper left") plt.show()
可以看到，在得到 X 的觀測值之後，史蒂文為圖書管理員的後臺概率增加了。同時，農民的後驗概率降低了。
但同時也看到，史蒂文為圖書管理員的後驗概率增加的並不是很多，史蒂文是農民的可能性依舊很大。
這表明，先驗概率對貝葉斯推斷影響是很大的，例如我們假設一臺伺服器平時不會被黑客攻陷（IOC），但是如果觀測到一系列的異常事件，被攻陷的後驗概率就會不斷增大，但這取決於觀測到的現象的數量和強度。
3. 一個從簡訊資料推斷行為模式的例子 - 使用計算機執行貝葉斯推斷
0x1：場景描述 - 我們的資料是什麼樣的？
下面的直方圖顯示了一段時間內的簡訊數量的統計。
# -*- coding: utf-8 -*- import numpy as np from matplotlib import pyplot as plt 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("count of text-msgs received") plt.title("Did the user's texting habits change over time?") plt.xlim(0, n_count_data) plt.show()
資料集下載連結。
0x2：問題定義 - 我們要從資料中挖掘出什麼東西？
我們的問題是：這段時間內使用者的行為出現了什麼變化嗎？或者是說，這段時間內是否出現了異常的行為模式，又具體是幾號發生的？
這特別像安全攻防領域經常面對的一個問題，即在一段日誌中進行資料探勘，從中找出是否存在異常行為模式（即入侵IOC）。
通常來說，一個直觀的想法可能會是這樣，我們通過ngram或者定長的時間視窗進行日誌聚合，然後通過人工經驗聚合中抽取一系列的特徵，並希望藉助有監督或者無監督的方法直接得到一個 0/1 判斷，即這個視窗期間內是否存在異常點。這是一種典型的端到端判斷場景。
但是其實我覺得應該更靈活地理解演算法，不一定每一個場景都適合使用判別式模型進行二分或者多分判斷，有時候，我們其實可以嘗試使用生成式模型，基於訓練資料，得到一個或多個概率分佈生成式。而多個生成式模型之間的分界點，很可能就是代表了一種模式躍遷，即行為異常。
0x3：使用什麼生成式模型進行資料擬合？
1. 二項分佈？正態分佈？泊松分佈？
二項分佈，正態分佈，泊松分佈這3個都是對自然界的事物在一定時間內的發生次數進行泛化模擬的分佈模型。實際上，它們3個之間相差並不是很大。
如果忽略分佈是離散還是連續的前提（二項分佈和泊松分佈一樣都是離散型概率分佈，正態分佈是連續型概率分佈），二項分佈與泊松分佈以及正態分佈至少在形狀上是十分接近的，也即兩邊低中部高。
當 n 足夠大，p 足夠小（事件間彼此獨立，事件發生的概率不算太大，事件發生的概率是穩定的），二項分佈逼近泊松分佈，λ=np；
一個被廣泛接受的經驗法則是如果 n≥20，p≤0.05，可以用泊松分佈來估計二項分佈值。
至於正態分佈是一個連續分佈 當實驗次數 n 再變大，幾乎可以看成連續時二項分佈和泊松分佈都可以用正態分佈來代替。
2. 單位時間內隨機事件發生的次數趨近於一個相對固定值 - 選泊松分佈
在日常生活中，大量事件我們可以預見其在時間視窗內的總數，但是無法準確知道具體某個時刻發生的次數。例如
1. 某醫院平均每小時出生3個嬰兒； 2. 某公司平均每10分鐘接到1個電話； 3. 某超市平均每天銷售4包 xx 牌奶粉； 4. 某網站平均每分鐘有2次訪問；
這些事件的特點就是，我們可以預估這些事件的總數，但是沒法具體知道發生時間。例如醫院具體下一小時的前10分鐘出現幾個嬰兒，有可能前10分鐘出現1個，也可能一下子出生3個。
泊松分佈就是描述某段時間內，事件具體的發生概率。
P 表示概率，N 表示某種函式關係，t 表示時間，n 表示數量，λ 表示事件的頻率。
從公式上看，我們得到幾個泊松分佈的基本判斷：
1. 單位時間內發生 n 次數的概率隨著時間 t 變化；
2. 概率受發生次數 n 的影響；
3. λ 相當於這件事的先驗屬性，在整個時間週期內都影響了最終的概率，λ 被稱為此分佈的一個超引數，它決定了這個分佈的形式。
λ 可以位任意正數，隨著 λ 的增加，得到大值的概率會增大。相反地，當 λ 減小時，得到小值的概率會增大。λ 可以被稱為 Poisson分佈的強度。
接下來兩個小時，一個嬰兒都不出現的概率是（已知 λ = 3）：
基本不可能發生。
接下來一個小時，至少出現兩個嬰兒的概率是：
即 80%。
泊松分佈的影象大概是這樣的：
可以看到，在頻率附近，事件發生的概率最高，然後向兩邊對稱下降。每小時出生 3 個嬰兒，這是最可能的結果，出生得越多或越少，就越不可能。
Poisson分佈的一個重要性質是：它的期望值等於它的引數。即：E[Z|λ] = λ。
下圖展示了不同 λ 取值下的概率質量分佈，可以看到，隨著 λ 的增加，取到大值的概率會增加。
# -*- coding: utf-8 -*- import numpy as np from matplotlib import pyplot as plt import scipy.stats as stats a = np.arange(16) poi = stats.poisson lambda_ = [1.5, 4.25] colours = ["#348ABD", "#A60628"] plt.bar(a, poi.pmf(a, lambda_[0]), color=colours[0], label="$\lambda = %.1f$" % lambda_[0], alpha=0.60, edgecolor=colours[0], lw="3") plt.bar(a, poi.pmf(a, lambda_[1]), color=colours[1], label="$\lambda = %.1f$" % lambda_[1], alpha=0.60, edgecolor=colours[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") plt.show()
3. 用單個泊松分佈就足夠擬合資料集了嗎？
好！我們現在已經確定要使用泊松分佈來對資料集進行擬合。接下來的問題是，我們是不是隻要用單個泊松分佈就足夠對這份資料集進行完美擬合了呢？
我們再來看看我們的問題場景。
我們這裡通過觀察資料，我們認為在某個時間段，發生了一個行為模式的突變，即在後期，收到簡訊的頻率好像提高了。
實際上，這已經進入到異常檢測的範疇內了，我們可以基於泊松分佈混合模型對這個資料集進行擬合。
4. 怎麼形式化描述泊松分佈引數 λ 的階躍情況？
我們不能確定引數 λ 的真實取值，但是可以猜測在某個時段 λ 增加了（λ取大值時，更容易得到較大的結果值），也即一天收到的簡訊條數比較多的概率增大了。
我們怎麼形式化地描述這個現象呢？假設在某一個時間段內(用τ表示)(τ可能也無法準確定義到具體的某一天)，引數?突然變大。
所以?可以取兩個值，在時段τ之前有一個，在τ時段之後有另一個。這種突然發生變化的情況稱之為轉換點。
我們以貝葉斯推斷的思想來看待這個問題：
1. 首先，存在一個 λ 階躍點只是一個先驗假設。事實上，可能並不存在這樣一個 λ 階躍點。 2. 第二，具體的階躍點 τ 是什麼時間呢？可能我們無法準確知道時間點，只能得到一個時間區間的概率分佈。
1）確定泊松分佈引數 λ 的先驗分佈
如果實際中 λ 沒有變化，則 
，那麼 λ 的先驗概率是相等的。
在貝葉斯推斷下，我們需要對不同可能取值的 λ 分配相應的先驗概率。那具體什麼樣的先驗概率對
來說才是最有可能的呢？
嚴格來說，我們對這件事幾乎沒有任何先驗知識，因為我們可能對這個資料集的屬主使用者行為模式一無所知。
注意到 λ 可以取任意正數，所以我們可以用指數分佈來表示 λ 的先驗分佈，因為指數分佈對任意正數都存在一個連續密度函式，這對模擬 λ 來說是一個很好的選擇。但是，指數分佈本身也有它自己的引數，所以這個引數也需要包括在我們的模型公式裡面。我們稱它為 α。不同的 α 決定了該指數分佈的形狀。
α 稱之為超引數。該引數能夠影響其他引數。按照最大熵演算法的原理，我們不希望對這個引數賦予太多的主觀色彩。
注意到我們計算樣本中計算平均值的公式：
，它正好就是我們要估計的超引數 α 的逆，所以，超引數 α 可以直接根據觀測樣本的計算得到。
2）確定階躍點 τ 的先驗分佈
我們對於引數τ，由於受到噪聲資料的影響，很難為它選擇合適的先驗。所以還是依照最大熵原理，我們假定每一天的先驗估計都是一致的。即：
在不知道任何先驗知識的情況下，等概論被證明是一種最好的先驗假設。
0x4：使用PyMC進行貝葉斯推斷
Cronin對概率程式設計有一段激動人心的描述：
換一種方式考慮這件事件，跟傳統的程式設計僅僅向前執行不同的是，概率程式設計既可以向前也可以向後執行。它通過向前執行來計算其中包含的所有關於世界的假設結果（例如對模型空間的描述），但它也從資料中向後執行，以約束可能的解釋。
在實踐中，許多概率程式設計系統將這些向前和向後的操作巧妙地交織在一起，以給出有效的最佳解釋。
1. 宣告模型中的隨機變數
按照前面我們的討論，α 取樣本中計算平均值的逆，而泊松分佈的的引數 λ 由 α決定。對階躍點先驗估計的 τ 引數，使用等概率估計（最大熵原理）。
# -*- coding: utf-8 -*- import numpy as np from matplotlib import pyplot as plt import pymc as pm 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("count of text-msgs received") plt.title("Did the user's texting habits change over time?") plt.xlim(0, n_count_data) alpha = 1.0 / count_data.mean() # Recall count_data is the # variable that holds our txt 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) print("Random output:", tau.random(), tau.random(), tau.random())
程式碼輸出的結果為：('tau Random output:', array(36), array(57), array(38))
表明我們的階躍點 τ 引數是等概率隨機生成的。
引數 λ 受引數 τ 影響，所以在階躍點兩邊的不同的泊松分佈引數 λ 也是隨機的，它根據 τ 隨機生成。
我們希望PyMC收集所有變數資訊，以及觀測資料，並從中產生一個model。PyMC可以藉助一個叫馬爾科夫鏈蒙特卡洛（MCMC）的演算法。得到引數 λ1，λ2，τ 的後驗估計
2. 利用MCMC根據觀測資料得到對模型引數的極大後驗估計結果
# -*- coding: utf-8 -*- import numpy as np from matplotlib import pyplot as plt import pymc as pm 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("count of text-msgs received") plt.title("Did the user's texting habits change over time?") plt.xlim(0, n_count_data) alpha = 1.0 / count_data.mean() # Recall count_data is the # variable that holds our txt 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) @pm.deterministic def lambda_(tau=tau, lambda_1=lambda_1, lambda_2=lambda_2): out = np.zeros(n_count_data) out[:tau] = lambda_1 # lambda before tau is lambda1 out[tau:] = lambda_2 # lambda after (and including) tau is lambda2 return out print("tau Random output:", tau.random(), tau.random(), tau.random()) #print("lambda distribution:", lambda_(tau, lambda_1, lambda_2)) observation = pm.Poisson("obs", lambda_, value=count_data, observed=True) model = pm.Model([observation, lambda_1, lambda_2, tau]) mcmc = pm.MCMC(model) mcmc.sample(40000, 10000, 1) lambda_1_samples = mcmc.trace('lambda_1')[:] lambda_2_samples = mcmc.trace('lambda_2')[:] tau_samples = mcmc.trace('tau')[:] # histogram of the samples: ax = plt.subplot(311) ax.set_autoscaley_on(False) plt.hist(lambda_1_samples, histtype='stepfilled', bins=30, alpha=0.85, label="posterior of $\lambda_1$", color="#A60628", normed=True) plt.legend(loc="upper left") plt.title(r"""Posterior distributions of the variables $\lambda_1,\;\lambda_2,\;\tau$""") plt.xlim([15, 30]) plt.xlabel("$\lambda_1$ value") ax = plt.subplot(312) ax.set_autoscaley_on(False) plt.hist(lambda_2_samples, histtype='stepfilled', bins=30, alpha=0.85, label="posterior of $\lambda_2$", color="#7A68A6", normed=True) plt.legend(loc="upper left") plt.xlim([15, 30]) plt.xlabel("$\lambda_2$ value") plt.subplot(313) w = 1.0 / tau_samples.shape[0] * np.ones_like(tau_samples) plt.hist(tau_samples, bins=n_count_data, alpha=1, label=r"posterior of $\tau$", color="#467821", weights=w, rwidth=2.) plt.xticks(np.arange(n_count_data)) plt.legend(loc="upper left") plt.ylim([0, .75]) plt.xlim([35, len(count_data) - 20]) plt.xlabel(r"$\tau$ (in days)") plt.ylabel("probability") plt.show()
0x5：得到的模型引數後驗分佈結果說明了什麼？
首先，我們可以看到，貝葉斯推斷得到的是一個概率分佈，不是一個準確的結果。我們使用分佈描述位置的 λ 和 τ。
分佈越廣，我們的後驗概率越小
1. 有很大概率存在異常行為模式的突變
我們也可以看到引數的合理值：λ1 大概是18，λ2 大概是23。這兩個 λ 的後驗分佈明顯是不同的，這表明該使用者接受簡訊的行為確實發生了變化。
2. 後驗分佈是我們假設的指數分佈嗎？
我們還注意到，λ 的後驗分佈並不是我們在模型建立時假設的指數分佈，事實上，λ 的後驗分佈並不是我們能從原始模型中可以辨別的任何分佈。
但是這恰恰就是概率程式設計的好處所在，我們藉助計算機繞開了複雜的數學建模過程。
3. 階躍點時間點
我們的分析結果返回了一個 τ 的一個概率分佈。在很多天中，τ 不太好確定，總體來說，我們可以認為在3~4天內可以認為是潛在的階躍點。
我們可以看到，在45天左右，有50%的把握可以說使用者的行為是有所改變的。這個點的概率分佈突然出現了一個峰值，我們可以推測產生這種情況的原因是：簡訊費用的降低，天氣提醒簡訊的訂閱，或者是一段感情的開始。
0x6：根據模型引數得到後驗樣本
在開始討論後驗樣本之前。首先需要明白的是，模型的後驗樣本並不是真實存在的事實，它更多地像是我們基於已經得到的概率分佈得到的一個期望值。因為很顯然，我們的模型不能很好地表達出真實物理中的隨機分佈以及噪音的影響。我們能得到的只能是一個對整體情況的概括。
我們用後驗樣本來回答一下下面這個問題：
在第 t （0 <= t <= 70）天中，期望收到多少條簡訊？
顯然，這個問題要求我們給出一個概括統計，即求每天簡訊數量的期望。
Poisson分佈的期望值等於它的引數 λ。因此，問題相當於：在時間 t 中，引數 λ 的期望值是多少？
# -*- coding: utf-8 -*-
import numpy as np
from matplotlib import pyplot as plt
import pymc as pm
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("count of text-msgs received")
plt.title("Did the user's texting habits change over time?")
plt.xlim(0, n_count_data)
alpha = 1.0 / count_data.mean() # Recall count_data is the
# variable that holds our txt 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)
@pm.deterministic
def lambda_(tau=tau, lambda_1=lambda_1, lambda_2=lambda_2):
out = np.zeros(n_count_data)
out[:tau] = lambda_1 # lambda before tau is lambda1
out[tau:] = lambda_2 # lambda after (and including) tau is lambda2
return out
print("tau Random output:", tau.random(), tau.random(), tau.random())
#print("lambda distribution:", lambda_(tau, lambda_1, lambda_2))
observation = pm.Poisson("obs", lambda_, value=count_data, observed=True)
model = pm.Model([observation, lambda_1, lambda_2, tau])
mcmc = pm.MCMC(model)
mcmc.sample(40000, 10000, 1)
lambda_1_samples = mcmc.trace('lambda_1')[:]
print "lambda_1_samples"
print lambda_1_samples
lambda_2_samples = mcmc.trace('lambda_2')[:]
print "lambda_2_samples"
print lambda_2_samples
tau_samples = mcmc.trace('tau')[:]
print "tau_samples"
print tau_samples
# tau_samples, lambda_1_samples, lambda_2_samples contain
# N samples from the corresponding posterior distribution
N = tau_samples.shape[0]
expected_texts_per_day = np.zeros(n_count_data)
for day in range(0, n_count_data):
# ix is a bool index of all tau samples corresponding to
# the switchpoint occurring prior to value of 'day'
ix = day < tau_samples # 如果天數在轉折點tau之前，取值lambda_1,否則取lambda_2,然後在取平均值
# Each posterior sample corresponds to a value for tau.
# for each day, that value of tau indicates whether we're "before"
# (in the lambda1 "regime") or
# "after" (in the lambda2 "regime") the switchpoint.
# by taking the posterior sample of lambda1/2 accordingly, we can average
# over all samples to get an expected value for lambda on that day.
# As explained, the "message count" random variable is Poisson distributed,
# and therefore lambda (the poisson parameter) is the expected value of
# "message count".
print "lambda_1_samples[ix]"
print lambda_1_samples[ix]
print "lambda_2_samples[~ix]"
print lambda_2_samples[~ix]
expected_texts_per_day[day] = (lambda_1_samples[ix].sum() # 這個期望值我們假設是服從泊松分佈的，分佈的期望值＝引數lambda。
+ lambda_2_samples[~ix].sum()) / N
plt.plot(range(n_count_data), expected_texts_per_day, lw=4, color="#E24A33",
label="expected number of text-messages received")
plt.xlim(0, n_count_data)
plt.xlabel("Day")
plt.ylabel("Expected # text-messages")
plt.title("Expected number of text-messages received")
plt.ylim(0, 60)
plt.bar(np.arange(len(count_data)), count_data, color="#348ABD", alpha=0.65,
label="observed texts per day")
plt.legend(loc="upper left");
plt.show()
在程式碼中可以看到，給定某天 t，我們對所有可能的 λ 求均值。
1. 基於後驗樣本，從統計學上確定兩個 λ 是否真的不一樣
我們基於引數 λ 的後驗分佈得到一個結論，λ 存在階躍點。但是如果這不是真相呢？有一部分分佈其實是重合的呢？這個可能性有多大呢？怎麼量化地評估這個可能性呢？
一個方法就是計算出 P（λ1 < λ2 | data），即在獲得觀察資料的前提下，λ1 的真實值比 λ2 小的概率。
print (lambda_1_samples < lambda_2_samples).mean() 0.999933333333
很顯然，有99.93%的概率說這兩個值是不相等的。
可以看到，這本質上是計算在 λ1 和 λ2 後驗概率分佈的重合程度。
0x7：擴充至兩個階躍點 - 能否藉助貝葉斯推理自動挖掘潛在的規律？
我們在章節的一開始就作出一個假設，即認為我們的資料集中蘊含了一個現象，即使用者的行為模式在某個時間點發生了突變。
但是是否存在可能說這個轉折點不只一個呢？
假設現在我們對一個轉折點表示很懷疑，我們現在認為使用者的行為發生了兩次改變．擴充之後，使用者的行為分為三個階段，三個泊松分佈對應三個λ1,λ2,λ3，兩個轉折點τ1,τ2。
# -*- coding: utf-8 -*- import numpy as np from matplotlib import pyplot as plt import pymc as pm 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("count of text-msgs received") plt.title("Did the user's texting habits change over time?") plt.xlim(0, n_count_data) alpha = 1.0 / count_data.mean() # Recall count_data is the # variable that holds our txt counts 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, lambda_1=lambda_1, lambda_2=lambda_2, lambda_3=lambda_3): out = np.zeros(n_count_data) out[:tau_1] = lambda_1 # lambda before tau is lambda1 out[tau_1:tau_2] = lambda_2 # lambda after (and including) tau is lambda2 out[tau_2:] = lambda_3 # lambda after (and including) tau is lambda2 return out observation = pm.Poisson("obs", lambda_, value=count_data, observed=True) model = pm.Model([observation, lambda_1, lambda_2, lambda_3, tau_1, tau_2]) mcmc = pm.MCMC(model) mcmc.sample(40000, 10000) lambda_1_samples = mcmc.trace('lambda_1')[:] lambda_2_samples = mcmc.trace('lambda_2')[:] lambda_3_samples = mcmc.trace('lambda_3')[:] tau_1_samples = mcmc.trace('tau_1')[:] tau_2_samples = mcmc.trace('tau_2')[:] # histogram of the samples: # lambda_1 ax = plt.subplot(311) ax.set_autoscaley_on(False) plt.hist(lambda_1_samples, histtype='stepfilled', bins=30, label="$\lambda_1$", normed=True) plt.legend(loc="upper left") plt.grid(True) plt.title(r"""Posterior distributions of the variables $\lambda_1,\;\lambda_2,\;\tau$""") # x軸座標範圍 plt.xlim([15, 30]) plt.xlabel("$\lambda_1$ value") # lambda_2 ax = plt.subplot(312) ax.set_autoscaley_on(False) plt.hist(lambda_2_samples, histtype='stepfilled', bins=30, label=" $\lambda_2$",color='#3009A6',normed=True) plt.legend(loc="upper left") plt.grid(True) plt.xlim([30, 90]) plt.xlabel("$\lambda_2$ value") # lambda_3 ax = plt.subplot(313) ax.set_autoscaley_on(False) plt.hist(lambda_3_samples, histtype='stepfilled', bins=30, label=" $\lambda_2$",color='#6A63A6',normed=True) plt.legend(loc="upper left") plt.grid(True) plt.xlim([15, 30]) plt.xlabel("$\lambda_3$ value") plt.show()
# -*- coding: utf-8 -*- import numpy as np from matplotlib import pyplot as plt import pymc as pm 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("count of text-msgs received") plt.title("Did the user's texting habits change over time?") plt.xlim(0, n_count_data) alpha = 1.0 / count_data.mean() # Recall count_data is the # variable that holds our txt counts 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, lambda_1=lambda_1, lambda_2=lambda_2, lambda_3=lambda_3): out = np.zeros(n_count_data) out[:tau_1] = lambda_1 # lambda before tau is lambda1 out[tau_1:tau_2] = lambda_2 # lambda after (and including) tau is lambda2 out[tau_2:] = lambda_3 # lambda after (and including) tau is lambda2 return out observation = pm.Poisson("obs", lambda_, value=count_data, observed=True) model = pm.Model([observation, lambda_1, lambda_2, lambda_3, tau_1, tau_2]) mcmc = pm.MCMC(model) mcmc.sample(40000, 10000) lambda_1_samples = mcmc.trace('lambda_1')[:] lambda_2_samples = mcmc.trace('lambda_2')[:] lambda_3_samples = mcmc.trace('lambda_3')[:] tau_1_samples = mcmc.trace('tau_1')[:] tau_2_samples = mcmc.trace('tau_2')[:] # histogram of the samples: # lambda_1 # tau_1 w = 1.0 / tau_1_samples.shape[0] * np.ones_like(tau_1_samples) plt.hist(tau_1_samples, bins=count_data, alpha=1, label=r"$\tau_1$",color="blue",weights=w ) plt.grid(True) plt.legend(loc="upper left") plt.xlabel(r"$\tau_1$ (in days)") plt.ylabel("probability") # tau_2 w = 1.0 / tau_2_samples.shape[0] * np.ones_like(tau_1_samples) plt.hist(tau_2_samples, bins=count_data, alpha=1, label=r"$\tau_2$",weights=w,color="red",) plt.xticks(np.arange(count_data)) plt.grid(True) plt.legend(loc="upper left") plt.ylim([0, 1.0]) plt.xlim([35, len(count_data) - 20]) plt.xlabel(r"$\tau_2$ (in days)") plt.ylabel("probability") plt.show()
貝葉斯推理展示了5個未知數的後驗推理結果。我們可以看到模型的轉折點大致在第45天和第47天的時候取得。這個是真實的情況嗎？我們的模型有可能過擬合了嗎？
確實，我們都可能對資料中有多少個轉折點抱有懷疑的態度。這意味著應該有多少個轉折點可以通過設定一個先驗分佈，並讓模型自己做決定。
這是一種非常重要的思想：我們對所有未知的變數都可以用一個先驗分佈來假設它，通過觀測樣本的訓練，得到最佳的後驗分佈。
Relevant Link:
https://blog.csdn.net/u014281392/article/details/79330286
4. 挑戰者號事故案例
# -*- coding: utf-8 -*- import numpy as np from matplotlib import pyplot as plt np.set_printoptions(precision=3, suppress=True) challenger_data = np.genfromtxt("data/challenger_data.csv", skip_header=1, usecols=[1, 2], missing_values="NA", delimiter=",") # drop the NA values challenger_data = challenger_data[~np.isnan(challenger_data[:, 1])] # plot it, as a function of temperature (the first column) print("Temp (F), O-Ring failure?") print(challenger_data) plt.scatter(challenger_data[:, 0], challenger_data[:, 1], s=75, color="k", alpha=0.5) plt.yticks([0, 1]) plt.ylabel("Damage Incident?") plt.xlabel("Outside temperature (Fahrenheit)") plt.title("Defects of the Space Shuttle O-Rings vs temperature") plt.show()
實驗資料如下：
Date,Temperature,Damage Incident 04/12/1981,66,0 11/12/1981,70,1 3/22/82,69,0 6/27/82,80,NA 01/11/1982,68,0 04/04/1983,67,0 6/18/83,72,0 8/30/83,73,0 11/28/83,70,0 02/03/1984,57,1 04/06/1984,63,1 8/30/84,70,1 10/05/1984,78,0 11/08/1984,67,0 1/24/85,53,1 04/12/1985,67,0 4/29/85,75,0 6/17/85,70,0 7/29/85,81,0 8/27/85,76,0 10/03/1985,79,0 10/30/85,75,1 11/26/85,76,0 01/12/1986,58,1 1/28/86,31,Challenger Accident
0x1：我們如何對資料分佈進行建模？
1. 通過非線性迴歸模擬樣本資料的函式表示
首先，通過簡單的觀察，這個樣本資料的自變數和因變數明顯不符合一個線性函式的關係。而且，我們甚至有理由懷疑它們之間存在一個階躍函式的關係，即在65°時發生階躍。但階躍函式太過武斷，不夠柔和，而且也沒有證據。
我們使用邏輯斯蒂線性迴歸，嘗試對觀測樣本進行擬合。
# -*- coding: utf-8 -*- # Code source: Gael Varoquaux # License: BSD 3 clause import numpy as np import matplotlib.pyplot as plt from sklearn import linear_model # data set X = np.array([ [66.], [70.], [69.], [80.], [68.], [67.], [72.], [73.], [70.], [57.], [63.], [70.], [78.], [67.], [53.], [67.], [75.], [70.], [81.], [76.], [79.], [75.], [76.], [58.] ]) y = np.array([ 0., 1., 0., 0., 0., 0., 0., 0., 0., 1., 1., 1., 0., 0., 1., 0., 0., 0., 0., 0., 0., 1., 0., 1. ]) print "X" print X.ravel() print "y" print y # and plot the data plt.figure(1, figsize=(8, 6)) plt.clf() plt.scatter(X.ravel(), y, color='black', zorder=20) X_test = np.linspace(50, 85, 300) # run the classifier clf = linear_model.LogisticRegression(C=1e5) clf.fit(X, y) def model(x): y_value = 1 / (1 + np.exp(-x)) return y_value # X_test * clf.coef_ + clf.intercept = 相當於: wx + b # model()負責轉換為邏輯斯蒂迴歸函式的 y 值 print "clf.coef: ", clf.coef_ print "clf.intercept_: ", clf.intercept_ y_value = model(X_test * clf.coef_ + clf.intercept_).ravel() plt.plot(X_test, y_value, color='red', linewidth=3) # 畫出模型擬合後的曲線 plt.ylabel('y') plt.xlabel('X') plt.xticks(range(50, 85)) plt.yticks([0, 0.5, 1]) plt.ylim(-.25, 1.25) plt.xlim(50, 85) plt.legend(('Logistic Regression Model'), loc="lower right", fontsize='small') plt.show()
我們得到了sklearn對樣本資料的邏輯斯蒂迴歸擬合結果，我們可以看到，一個很明顯的結論是：隨著溫度下降，出現故障的風險會不斷提高，在50°時甚至接近於1，即必然發生。
同時我們也列印出了擬合的模型引數：
clf.coef: [[-0.22842526]] clf.intercept_: [ 14.77684766]
這兩個引數決定這個邏輯函式的分佈形態，我們後面還會看到，我們用貝葉斯推斷的方式，也能獲取到這2個引數的估計結果。
注意，這裡coef是包含負號的，所以我們可以認為sklearn推匯出來的結果是：1 / 1 + e ^ 0.2284 + 14.77
2. 通過數學建模和貝葉斯推理的方式進行資料擬合
1）溫度和事故發生的概率如何建模？
? 是這個模型的引數，是我們需要從資料中進行推理得到的。下面畫出? = 1,3, −5時的?(?)的曲線
# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt def logistic(x, beta): return 1.0 / (1.0 + np.exp(beta * x)) x = np.linspace(-4, 4, 100) plt.plot(x, logistic(x, 1), label=r"$\beta = 1$") plt.plot(x, logistic(x, 3), label=r"$\beta = 3$") plt.plot(x, logistic(x, -5), label=r"$\beta = -5$") plt.legend() plt.show()
畫這個圖的意義是什麼呢？其實沒啥意義，這只是隨便估計了一下引數，但起碼我們看到，引數 ? 起碼必須是正數，而且 ? 應該是一個傾向於小的值。當然，這只是我們的猜測，具體的結果要進行貝葉斯程式設計後推導得到。
此外，在上圖中，在 ? = 0 附件概率會發生改變，但是在實際的資料中在 65 到 70 這個範圍內概率發生改變。所以還需要在 logistic 函式中新增一個偏置項，函式變為：
不同的 ? 對應的 logistic 函式如下
# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt def logistic(x, beta, alpha=0): return 1.0 / (1.0 + np.exp(np.dot(beta, x) + alpha)) x = np.linspace(-4, 4, 100) plt.plot(x, logistic(x, 1), label=r"$\beta = 1$", ls="--", lw=1) plt.plot(x, logistic(x, 3), label=r"$\beta = 3$", ls="--", lw=1) plt.plot(x, logistic(x, -5), label=r"$\beta = -5$", ls="--", lw=1) plt.plot(x, logistic(x, 1, 1), label=r"$\beta = 1, \alpha = 1$",color="#348ABD") plt.plot(x, logistic(x, 3, -2), label=r"$\beta = 3, \alpha = -2$",color="#A60628") plt.plot(x, logistic(x, -5, 7), label=r"$\beta = -5, \alpha = 7$",color="#7A68A6") plt.legend(loc="lower left"); plt.show()
? 引數控制著曲線在橫軸上左右移動，這也就是把 ? 稱為偏置的原因。
2）邏輯函式的引數如何建模？
讀者小夥伴們注意，在貝葉斯概率程式設計領域，我們要習慣用概率分佈的思維方式來思考問題。對於邏輯函式的引數 ? 和引數 ? 來說，並沒有限定要為正數，或者在某個範圍內要相對的大。基本上來說，我們不知道這兩個引數的先驗，因此我們用正態分佈來模擬它們是一個相對比較合適的選擇。
，正態分佈有兩個引數均值?，方差?。也可以用?2代替1⁄?，它們二者互為倒數。?越小，正態分佈的曲線越平坦(不確定性越大)，反之亦然(不確定性越小)。?的取值範圍為正數。對模型引數建模之後，就可以代入上一小節溫度和事故發生概率的模型中。
import pymc as pm temperature = challenger_data[:, 0] D = challenger_data[:, 1] # defect or not? # notice the`value` here. We explain why below. beta = pm.Normal("beta", 0, 0.001, value=0) alpha = pm.Normal("alpha", 0, 0.001, value=0) @pm.deterministic def p(t=temperature, alpha=alpha, beta=beta): return 1.0 / (1. + np.exp(beta * t + alpha))
3）我們的觀測資料如何和模型引數融合建模？
前兩個小節，邏輯函式和正態分佈都是針對單個樣本的輸入和輸出進行建模，我們的實驗觀察是多個樣本，那麼我們怎麼將我們的觀測資料和我們的概率值（通過邏輯函式模擬的）聯絡起來呢？
如果按照我們在機器學習程式設計中的做法，這個時候，我們已經對輸入和輸出進行了建模（邏輯函式），那接下來的做法就是選擇一個損失函式（目標函式）（這裡很可能就是0-1損失），通過計算觀察樣本的整體損失並通過極值求導的方式來得到模型的最優擬合引數。
但是記住，我們現在是在貝葉斯概率程式設計領域中，貝葉斯推斷理論認為，輸入變數和輸出變數之間是不存在絕對的因果推導關係的，我們所謂的觀測結果只是在在一個概率分佈下的一次（一批）具體實驗具現化結果，我們對任何的因果關係都要持懷疑態度，即採集概率分佈的方式來進行建模。
回到我們這個具體的問題場景，在一定的溫度下，我們可以得到事件發生的概率，但是事件是否真的發生，有兩種可能，發生 or 不發生。因此我們可以用伯努利分佈來對“事件發生概率 -> 是否真的發生”進行概率建模。
，其中 p(t) 是我們的邏輯函式，ti 是觀測到的溫度值。4）通過pymc對觀測樣本進行模型引數的後驗估計
# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt import pymc as pm np.set_printoptions(precision=3, suppress=True) challenger_data = np.genfromtxt("data/challenger_data.csv", skip_header=1, usecols=[1, 2], missing_values="NA", delimiter=",") # drop the NA values challenger_data = challenger_data[~np.isnan(challenger_data[:, 1])] # plot it, as a function of temperature (the first column) print("Temp (F), O-Ring failure?") print(challenger_data) plt.scatter(challenger_data[:, 0], challenger_data[:, 1], s=75, color="k", alpha=0.5) plt.yticks([0, 1]) plt.ylabel("Damage Incident?") plt.xlabel("Outside temperature (Fahrenheit)") plt.title("Defects of the Space Shuttle O-Rings vs temperature") temperature = challenger_data[:, 0] D = challenger_data[:, 1] # defect or not? # notice the`value` here. We explain why below. beta = pm.Normal("beta", 0, 0.001, value=0) alpha = pm.Normal("alpha", 0, 0.001, value=0) @pm.deterministic def p(t=temperature, alpha=alpha, beta=beta): return 1.0 / (1. + np.exp(beta * t + alpha)) # connect the probabilities in `p` with our observations through a # Bernoulli random variable. observed = pm.Bernoulli("bernoulli_obs", p, value=D, observed=True) model = pm.Model([observed, beta, alpha]) # Mysterious code to be explained in Chapter 3 map_ = pm.MAP(model) map_.fit() mcmc = pm.MCMC(model) mcmc.sample(120000, 100000, 2) alpha_samples = mcmc.trace('alpha')[:, None] # best to make them 1d beta_samples = mcmc.trace('beta')[:, None] # histogram of the samples: plt.subplot(211) plt.title(r"Posterior distributions of the variables $\alpha, \beta$") plt.hist(beta_samples, histtype='stepfilled', bins=35, alpha=0.85, label=r"posterior of $\beta$", color="#7A68A6", normed=True) plt.legend() plt.subplot(212) plt.hist(alpha_samples, histtype='stepfilled', bins=35, alpha=0.85, label=r"posterior of $\alpha$", color="#A60628", normed=True) plt.legend() plt.show()
0x2：我們從模型引數的後驗概率估計中可以得到什麼資訊？
1. 溫度對事故的發生肯定是有影響的
所有的 ? 樣本值都大於 0，說明溫度這一自變數是因變數，一定有所影響。
如果後驗概率集中在 0 附件，假設? = 0，說明溫度對失敗的概率沒有影響。
2. ? 一定是小於 0 的
所有的 ? 後驗值都是負的且遠離 0，可以確信 ? 一定是小於 0 的。
3. 我們對模型引數的真實取值依然無法 100% 確定
鑑於資料的分佈情況，我們不確定真實的引數?, ?到底是什麼。當然，考慮到樣本集比較小，以及發生缺陷和沒有缺陷事故的大量重疊，得到這樣的結果也是正常的。
0x3：基於模型引數的後驗估計反推事件發生的概率的後驗估計
像之前的簡訊行為後驗估計一樣，我們現在得到了邏輯函式的模型引數的後驗估計，我們現在可以反過來做，把估計得到的結果當做輸入，得到模型的輸出結果，即事件發生的後驗概率。
# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt import pymc as pm np.set_printoptions(precision=3, suppress=True) challenger_data = np.genfromtxt("data/challenger_data.csv", skip_header=1, usecols=[1, 2], missing_values="NA", delimiter=",") # drop the NA values challenger_data = challenger_data[~np.isnan(challenger_data[:, 1])] # plot it, as a function of temperature (the first column) print("Temp (F), O-Ring failure?") print(challenger_data) plt.scatter(challenger_data[:, 0], challenger_data[:, 1], s=75, color="k", alpha=0.5) plt.yticks([0, 1]) plt.ylabel("Damage Incident?") plt.xlabel("Outside temperature (Fahrenheit)") plt.title("Defects of the Space Shuttle O-Rings vs temperature") temperature = challenger_data[:, 0] D = challenger_data[:, 1] # defect or not? # notice the`value` here. We explain why below. beta = pm.Normal("beta", 0, 0.001, value=0) alpha = pm.Normal("alpha", 0, 0.001, value=0) @pm.deterministic def p(t=temperature, alpha=alpha, beta=beta): return 1.0 / (1. + np.exp(beta * t + alpha)) def logistic(x, beta, alpha=0): return 1.0 / (1.0 + np.exp(np.dot(beta, x) + alpha)) # connect the probabilities in `p` with our observations through a # Bernoulli random variable. observed = pm.Bernoulli("bernoulli_obs", p, value=D, observed=True) model = pm.Model([observed, beta, alpha]) # Mysterious code to be explained in Chapter 3 map_ = pm.MAP(model) map_.fit() mcmc = pm.MCMC(model) mcmc.sample(120000, 100000, 2) alpha_samples = mcmc.trace('alpha')[:, None] # best to make them 1d beta_samples = mcmc.trace('beta')[:, None] # 檢視對既定的溫度取值的期望概率，即我們隊所有後驗樣本取均值，得到一個最大後驗 p(t) t = np.linspace(temperature.min() - 5, temperature.max() + 5, 50)[:, None] p_t = logistic(t.T, beta_samples, alpha_samples) print "p_t" print p_t mean_prob_t = p_t.mean(axis=0) plt.plot(t, mean_prob_t, lw=3, label="average posterior \nprobability of defect") # 求均值即求對特定溫度的下事件發生的期望概率 plt.plot(t, p_t[0, :], ls="--", label="realization from posterior") # 在不同的模型引數後驗估計下，事件發生概率的後驗估計 plt.plot(t, p_t[-2, :], ls="--", label="realization from posterior") # 在不同的模型引數後驗估計下，事件發生概率的後驗估計 plt.scatter(temperature, D, color="k", s=50, alpha=0.5) plt.title("Posterior expected value of probability of defect; \ plus realizations") plt.legend(loc="lower left") plt.ylim(-0.1, 1.1) plt.xlim(t.min(), t.max()) plt.ylabel("probability") plt.xlabel("temperature") plt.show()
1. 95%置信區間（CI）
一個有趣的問題是：在哪個溫度時，我們對缺陷發生的概率是最不確定的？
我們通過編碼來看看，期望會的曲線和每個點對應的95%的置信區間（CI）
# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt import pymc as pm from scipy.stats.mstats import mquantiles np.set_printoptions(precision=3, suppress=True) challenger_data = np.genfromtxt("data/challenger_data.csv", skip_header=1, usecols=[1, 2], missing_values="NA", delimiter=",") # drop the NA values challenger_data = challenger_data[~np.isnan(challenger_data[:, 1])] # plot it, as a function of temperature (the first column) print("Temp (F), O-Ring failure?") print(challenger_data) plt.scatter(challenger_data[:, 0], challenger_data[:, 1], s=75, color="k", alpha=0.5) plt.yticks([0, 1]) plt.ylabel("Damage Incident?") plt.xlabel("Outside temperature (Fahrenheit)") plt.title("Defects of the Space Shuttle O-Rings vs temperature") temperature = challenger_data[:, 0] D = challenger_data[:, 1] # defect or not? # notice the`value` here. We explain why below. beta = pm.Normal("beta", 0, 0.001, value=0) alpha = pm.Normal("alpha", 0, 0.001, value=0) @pm.deterministic def p(t=temperature, alpha=alpha, beta=beta): return 1.0 / (1. + np.exp(beta * t + alpha)) def logistic(x, beta, alpha=0): return 1.0 / (1.0 + np.exp(np.dot(beta, x) + alpha)) # connect the probabilities in `p` with our observations through a # Bernoulli random variable. observed = pm.Bernoulli("bernoulli_obs", p, value=D, observed=True) model = pm.Model([observed, beta, alpha]) # Mysterious code to be explained in Chapter 3 map_ = pm.MAP(model) map_.fit() mcmc = pm.MCMC(model) mcmc.sample(120000, 100000, 2) alpha_samples = mcmc.trace('alpha')[:, None] # best to make them 1d beta_samples = mcmc.trace('beta')[:, None] # 檢視對既定的溫度取值的期望概率，即我們隊所有後驗樣本取均值，得到一個最大後驗 p(t) t = np.linspace(temperature.min() - 5, temperature.max() + 5, 50)[:, None] p_t = logistic(t.T, beta_samples, alpha_samples) print "p_t" print p_t mean_prob_t = p_t.mean(axis=0) # vectorized bottom and top 2.5% quantiles for "confidence interval" qs = mquantiles(p_t, [0.025, 0.975], axis=0) plt.fill_between(t[:, 0], *qs, alpha=0.7, color="#7A68A6") plt.plot(t[:, 0], qs[0], label="95% CI", color="#7A68A6", alpha=0.7) plt.plot(t, mean_prob_t, lw=1, ls="--", color="k", label="average posterior \nprobability of defect") plt.xlim(t.min(), t.max()) plt.ylim(-0.02, 1.02) plt.legend(loc="lower left") plt.scatter(temperature, D, color="k", s=50, alpha=0.5) plt.xlabel("temp, $t$") plt.ylabel("probability estimate") plt.title("Posterior probability estimates given temp. $t$") plt.show()
95%置信區間，或者稱95% CI，在圖中用紫色區塊顯示。
我們知道，pymc得到模型引數 ?, ? 並不是一個固定的值，而是一個概率分佈，而不同的 ?, ? 得到不同的邏輯函式，因此，對每一個溫度 t，會有不同的事故發生概率後驗值。
上圖紫色區域的含義是：對每一個溫度值，它都包含了95%的分佈，例如在 65 度時，有 95%的確信度相信事故發生的概率範圍為 0.25 到 0.75 之間。
更一般情況，當溫度到達 60 度時，CI’s 快速變窄；當溫度超過 70 度時 CI’s 也有類似的情況。這將給我提供接下來如何分析的思路。
1. 首先，在低溫區間，尤其是50°以下，95%置信區間是很窄的，說明在低溫時候發生事故的概率是非常大的，且這個事實非常確定； 2. 我們可能在溫度區間 60-65 上做更多 的 O 形密封圈的測試，這樣就能得到更好的估計概率。
2. 挑戰者號事故當天發生了什麼
挑戰者號發生災難的那天，外部的溫度是 31 華氏度。在此時的溫度下，發生災難的後驗分佈如何?
還是一樣注意，貝葉斯推斷得到的後驗分佈是一個概率分佈（這點在剛開始學習概率程式設計時候特別不適應），因此分佈圖如下。
# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt import pymc as pm from scipy.stats.mstats import mquantiles np.set_printoptions(precision=3, suppress=True) challenger_data = np.genfromtxt("data/challenger_data.csv", skip_header=1, usecols=[1, 2], missing_values="NA", delimiter=",") # drop the NA values challenger_data = challenger_data[~np.isnan(challenger_data[:, 1])] # plot it, as a function of temperature (the first column) print("Temp (F), O-Ring failure?") print(challenger_data) plt.scatter(challenger_data[:, 0], challenger_data[:, 1], s=75, color="k", alpha=0.5) plt.yticks([0, 1]) plt.ylabel("Damage Incident?") plt.xlabel("Outside temperature (Fahrenheit)") plt.title("Defects of the Space Shuttle O-Rings vs temperature") temperature = challenger_data[:, 0] D = challenger_data[:, 1] # defect or not? # notice the`value` here. We explain why below. beta = pm.Normal("beta", 0, 0.001, value=0) alpha = pm.Normal("alpha", 0, 0.001, value=0) @pm.deterministic def p(t=temperature, alpha=alpha, beta=beta): return 1.0 / (1. + np.exp(beta * t + alpha)) def logistic(x, beta, alpha=0): return 1.0 / (1.0 + np.exp(np.dot(beta, x) + alpha)) # connect the probabilities in `p` with our observations through a # Bernoulli random variable. observed = pm.Bernoulli("bernoulli_obs", p, value=D, observed=True) model = pm.Model([observed, beta, alpha]) # Mysterious code to be explained in Chapter 3 map_ = pm.MAP(model) map_.fit() mcmc = pm.MCMC(model) mcmc.sample(120000, 100000, 2) alpha_samples = mcmc.trace('alpha')[:, None] # best to make them 1d beta_samples = mcmc.trace('beta')[:, None] prob_31 = logistic(31, beta_samples, alpha_samples) plt.xlim(0.995, 1) plt.hist(prob_31, bins=1000, normed=True, histtype='stepfilled') plt.title("Posterior distribution of probability of defect, given $t = 31$") plt.xlabel("probability of defect occurring in O-ring"); plt.show()
看起來挑戰者號必定會發生由於 O 形密封環導致的災難。
0x4：我們的模型適用嗎？
接下來的問題是，如何驗證我們的模型的吻合度是否不好？
# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt import pymc as pm from scipy.stats.mstats import mquantiles np.set_printoptions(precision=3, suppress=True) challenger_data = np.genfromtxt("data/challenger_data.csv", skip_header=1, usecols=[1, 2], missing_values="NA", delimiter=",") # drop the NA values challenger_data = challenger_data[~np.isnan(challenger_data[:, 1])] # plot it, as a function of temperature (the first column) print("Temp (F), O-Ring failure?") print(challenger_data) plt.scatter(challenger_data[:, 0], challenger_data[:, 1], s=75, color="k", alpha=0.5) plt.yticks([0, 1]) plt.ylabel("Damage Incident?") plt.xlabel("Outside temperature (Fahrenheit)") plt.title("Defects of the Space Shuttle O-Rings vs temperature") temperature = challenger_data[:, 0] D = challenger_data[:, 1] # defect or not? # notice the`value` here. We explain why below. beta = pm.Normal("beta", 0, 0.001, value=0) alpha = pm.Normal("alpha", 0, 0.001, value=0) @pm.deterministic def p(t=temperature, alpha=alpha, beta=beta): return 1.0 / (1. + np.exp(beta * t + alpha)) def logistic(x, beta, alpha=0): return 1.0 / (1.0 + np.exp(np.dot(beta, x) + alpha)) # connect the probabilities in `p` with our observations through a # Bernoulli random variable. observed = pm.Bernoulli("bernoulli_obs", p, value=D, observed=True) model = pm.Model([observed, beta, alpha]) # Mysterious code to be explained in Chapter 3 map_ = pm.MAP(model) map_.fit() mcmc = pm.MCMC(model) mcmc.sample(120000, 100000, 2) alpha_samples = mcmc.trace('alpha')[:, None] # best to make them 1d beta_samples = mcmc.trace('beta')[:, None] simulated = pm.Bernoulli("bernoulli_sim", p) N = 10000 mcmc = pm.MCMC([simulated, alpha, beta, observed]) mcmc.sample(N) simulations = mcmc.trace("bernoulli_sim")[:] print(simulations.shape) plt.title("Simulated dataset using posterior parameters") for i in range(4): ax = plt.subplot(4, 1, i + 1) plt.scatter(temperature, simulations[1000 * i, :], color="k", s=50, alpha=0.6) plt.show()
我們希望從統計上評估上面的模擬資料和我們的原始觀測資料是否真的很相似，即我們希望評估出模型模擬得究竟有多好。需要找一種量化的方式來量化衡量這裡所謂的“好”。
1. 分離圖
我們使用作圖的方式完成這件事，下面的圖形測試是一種用於邏輯迴歸擬合程度的資料視覺化方法。這種圖被稱為分離圖。分離圖可以讓使用者用一種圖形化的方法對比不同的模型並從中選出最適合的。
# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt import pymc as pm from scipy.stats.mstats import mquantiles np.set_printoptions(precision=3, suppress=True) challenger_data = np.genfromtxt("data/challenger_data.csv", skip_header=1, usecols=[1, 2], missing_values="NA", delimiter=",") # drop the NA values challenger_data = challenger_data[~np.isnan(challenger_data[:, 1])] temperature = challenger_data[:, 0] D = challenger_data[:, 1] # defect or not? # notice the`value` here. We explain why below. beta = pm.Normal("beta", 0, 0.001, value=0) alpha = pm.Normal("alpha", 0, 0.001, value=0) @pm.deterministic def p(t=temperature, alpha=alpha, beta=beta): return 1.0 / (1. + np.exp(beta * t + alpha)) def logistic(x, beta, alpha=0): return 1.0 / (1.0 + np.exp(np.dot(beta, x) + alpha)) # connect the probabilities in `p` with our observations through a # Bernoulli random variable. observed = pm.Bernoulli("bernoulli_obs", p, value=D, observed=True) model = pm.Model([observed, beta, alpha]) # Mysterious code to be explained in Chapter 3 map_ = pm.MAP(model) map_.fit() mcmc = pm.MCMC(model) mcmc.sample(120000, 100000, 2) alpha_samples = mcmc.trace('alpha')[:, None] # best to make them 1d beta_samples = mcmc.trace('beta')[:, None] simulated = pm.Bernoulli("bernoulli_sim", p) N = 10000 mcmc = pm.MCMC([simulated, alpha, beta, observed]) mcmc.sample(N) simulations = mcmc.trace("bernoulli_sim")[:] print(simulations.shape) # 對每一個模型，給定溫度，計算後驗模擬產生值1的次數比例，並對所有返回的模擬值取均值，這樣我們得到在每一個資料點上發生缺陷的後驗可能 posterior_probability = simulations.mean(axis=0) print("posterior prob of defect | realized defect ") for i in range(len(D)): print("%.2f | %d" % (posterior_probability[i], D[i])) # 接下來我們根據後驗概率對每一列進行排序 ix = np.argsort(posterior_probability) print("probb | defect ") for i in range(len(D)): print("%.2f | %d" % (posterior_probability[ix[i]], D[ix[i]]))
我們可以在一個圖中更好地視覺化展示這些資料。
# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt import pymc as pm from scipy.stats.mstats import mquantiles from separation_plot import separation_plot np.set_printoptions(precision=3, suppress=True) challenger_data = np.genfromtxt("data/challenger_data.csv", skip_header=1, usecols=[1, 2], missing_values="NA", delimiter=",") # drop the NA values challenger_data = challenger_data[~np.isnan(challenger_data[:, 1])] temperature = challenger_data[:, 0] D = challenger_data[:, 1] # defect or not? # notice the`value` here. We explain why below. beta = pm.Normal("beta", 0, 0.001, value=0) alpha = pm.Normal("alpha", 0, 0.001, value=0) @pm.deterministic def p(t=temperature, alpha=alpha, beta=beta): return 1.0 / (1. + np.exp(beta * t + alpha)) def logistic(x, beta, alpha=0): return 1.0 / (1.0 + np.exp(np.dot(beta, x) + alpha)) # connect the probabilities in `p` with our observations through a # Bernoulli random variable. observed = pm.Bernoulli("bernoulli_obs", p, value=D, observed=True) model = pm.Model([observed, beta, alpha]) # Mysterious code to be explained in Chapter 3 map_ = pm.MAP(model) map_.fit() mcmc = pm.MCMC(model) mcmc.sample(120000, 100000, 2) alpha_samples = mcmc.trace('alpha')[:, None] # best to make them 1d beta_samples = mcmc.trace('beta')[:, None] simulated = pm.Bernoulli("bernoulli_sim", p) N = 10000 mcmc = pm.MCMC([simulated, alpha, beta, observed]) mcmc.sample(N) simulations = mcmc.trace("bernoulli_sim")[:] print(simulations.shape) # 對每一個模型，給定溫度，計算後驗模擬產生值1的次數比例，並對所有返回的模擬值取均值，這樣我們得到在每一個資料點上發生缺陷的後驗可能 posterior_probability = simulations.mean(axis=0) # 視覺化展示 separation_plot(posterior_probability, D) plt.show()
蛇形線表示排序後的後驗概率，藍色柱子表示真實發生的缺陷（觀測結果），空的地方表示沒有發生缺陷。
隨著概率增加，可以看到事故發生的概率越來越大，在圖中的右邊，說明後驗概率越來越大（因為蛇形線接近 1），表示更多的事故發生。理想情況所有的藍色柱子都要靠近右邊，這裡面會存在一些偏差，偏差存在的地方預測也就缺失了。
垂直的線表示在這個模型下，所能觀測到的事故次數。這根垂直的線可以讓我們比較模型預測出的事故次數和真實事故次數各是多少。
將這個模型的分離圖和其它模型的分離圖進行對比，會發現更多的有用資訊。下面將上面的模型與其它三個進行對比：
# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt import pymc as pm from scipy.stats.mstats import mquantiles from separation_plot import separation_plot np.set_printoptions(precision=3, suppress=True) challenger_data = np.genfromtxt("data/challenger_data.csv", skip_header=1, usecols=[1, 2], missing_values="NA", delimiter=",") # drop the NA values challenger_data = challenger_data[~np.isnan(challenger_data[:, 1])] temperature = challenger_data[:, 0] D = challenger_data[:, 1] # defect or not? # notice the`value` here. We explain why below. beta = pm.Normal("beta", 0, 0.001, value=0) alpha = pm.Normal("alpha", 0, 0.001, value=0) @pm.deterministic def p(t=temperature, alpha=alpha, beta=beta): return 1.0 / (1. + np.exp(beta * t + alpha)) def logistic(x, beta, alpha=0): return 1.0 / (1.0 + np.exp(np.dot(beta, x) + alpha)) # connect the probabilities in `p` with our observations through a # Bernoulli random variable. observed = pm.Bernoulli("bernoulli_obs", p, value=D, observed=True) model = pm.Model([observed, beta, alpha]) # Mysterious code to be explained in Chapter 3 map_ = pm.MAP(model) map_.fit() mcmc = pm.MCMC(model) mcmc.sample(120000, 100000, 2) alpha_samples = mcmc.trace('alpha')[:, None] # best to make them 1d beta_samples = mcmc.trace('beta')[:, None] simulated = pm.Bernoulli("bernoulli_sim", p) N = 10000 mcmc = pm.MCMC([simulated, alpha, beta, observed]) mcmc.sample(N) simulations = mcmc.trace("bernoulli_sim")[:] print(simulations.shape) # 對每一個模型，給定溫度，計算後驗模擬產生值1的次數比例，並對所有返回的模擬值取均值，這樣我們得到在每一個資料點上發生缺陷的後驗可能 posterior_probability = simulations.mean(axis=0) # Our temperature-dependent model separation_plot(posterior_probability, D) plt.title("Temperature-dependent model") plt.show() # Perfect model # i.e. the probability of defect is equal to if a defect occurred or not. p = D separation_plot(p, D) plt.title("Perfect model") plt.show() # random predictions p = np.random.rand(23) separation_plot(p, D) plt.title("Random model") plt.show() # constant model constant_prob = 7. / 23 * np.ones(23) separation_plot(constant_prob, D) plt.title("Constant-prediction model") plt.show()
1. 如果模型很完美，當事故確實發生時，它預測的後驗概率是 1 2. 一個隨機模型預測的結果完全是隨機的，將和溫度沒有任何關係 3. 對於一個常數模型，即P(? = 1|?) = ?, ∀?。?的最優值是觀測到的事故發生的概率，即? = 7⁄23
在隨機模型中，可以看到，隨著概率的增加，事故並沒有向右側聚集。常數模型也是這個情況。
完美模型中沒有顯示概率曲線。因為它一直位於柱子的頂端。此處的完美模型只是用於演示，並不能從中推斷出什麼科學依據。
在實際專案中，我們要評估的是我們的後驗概率估計的分離圖是如何介於隨機模型和完美模擬之間的，自然，越靠近完美模型意味著我們的模型擬合優度越好。