貝葉斯推斷 && 機率編程初探

時間 2019-12-05

標籤貝葉推斷機率編程初探简体版

原文原文鏈接

1. 寫在以前的話

0x1：貝葉斯推斷的思想

咱們從一個例子開始咱們本文的討論。小明是一個編程老手，可是依然堅信bug仍有可能在代碼中存在。因而，在實現了一段特別難的算法以後，他開始決定先來一個簡單的測試用例，這個用例經過了。接着，他用了一個稍微複雜的測試用例，再次經過了。接下來更難的測試用例也經過了，這時，小明開始以爲這段代碼出現bug的可能性大大大大下降了....git

上面這段白話文中，已經包含了最質樸的貝葉斯思想了！簡單來講，貝葉斯推斷是經過新獲得的證據不斷地更新咱們的信念。github

貝葉斯推斷不多會作出絕對的判斷，但能夠作出很是可信的判斷。在上面的例子中，小明永遠沒法100%確定代碼是完好陷的，除非他窮舉測試了每一種可能出現的情形，這在實踐中幾乎是不可能的。可是，咱們能夠對代碼進行大量的測試，若是每一次測試都經過了，咱們「更」有把握以爲這段代碼是沒問題的。算法

貝葉斯推斷的工做方式就在這裏：咱們會隨着新的證據不斷更新以前的信念，但不多作出絕對的判斷，除非全部其餘的可能都被一一排除。編程

0x2：貝葉斯思惟和頻率派思惟的區別

1. 頻率派關於機率的解釋 - 統計

頻率派是更古典的統計學派，它們認爲機率是事件在長時間內發生的機率。例如，在頻率派的哲學語境裏，飛機事故的機率指的是長期來看，飛機事故的機率值。安全

對許多事件來講，這樣解釋是有邏輯的。可是對某些沒有長期頻率的事件來講，這樣的解釋是難以理解的。例如，在總統選舉時，某個候選人的獲選的機率是難以用統計的方式來回答的，由於這個候選人的選舉可能只發生一次。頻率論爲了解決這個問題，提出了「替代現實」的說法，從而用在全部這些替代的「現實」中發生的機率定義了這個機率。服務器

2. 貝葉斯派關於機率的解釋 - 信息

貝葉斯派對機率有更直觀的解釋。它把機率解釋成是對事件發生的信心。簡單地說，機率是觀點的概述。框架

若是某人把機率 0 賦給某個事件的發生，代表他徹底肯定此事不會發生；相反，若是某人賦的機率值是 1，則代表他十分確定此事必定會發生。dom

0 和 1 之間的機率值能夠代表此事可能發生的權重。機器學習

這個機率定義能夠和飛機事故機率一致，若是除去全部外部信息，一我的對飛機事故發生的信心應該等同於他了解到的飛機事故的頻率。一樣，貝葉斯機率定義也能適用於總統選舉，而且使得機率（信心）更加有意義：即你對候選人 A 獲勝的信息有多少？ide

對咱們而言，將對一個事件發生的信息等同於機率是十分天然的，這正是咱們長期以來同世界打交道的方式。咱們只能瞭解到部分的真相，但能夠經過不斷收集證據來完善咱們對事物的觀念。

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}$

# -*- 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. 挑戰者號事故案例

1986 年 1 月 28 日，美國航天飛機計劃以災難了結。在起飛後不就，航天飛機的其中一個助推器發生爆炸，七名宇航員所有遇難。總統的調查委員會得出的結論是，這場災難是因爲火箭推動器上的一個密封圈的失效引發的，致使失效的緣由是密封圈對外界許多因素很是敏感，包括溫度。

從以前的 24 次發射中，能夠獲得這次失效的 O 形密封圈的 23 批次的數據(其中一個密封圈丟失在海上了)。在挑戰者發射的前一天晚上，對這 23 組數據進行了專門的討論，不幸的是得出只有七組數據反映出密封圈會形成損壞事故，這說明沒有明顯的趨勢代表密封圈會出問題。數據展現以下:

# -*- 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：咱們如何對數據分佈進行建模？

咱們知道，表達自變量和因變量之間關係的方式有2種，一種是函數表示，即：Y = F (X)；另外一種是機率表示，即：Y = P (X)。

兩種方法本質上都是同樣的，只是說函數表示傾向於獲得一個判別結果（判別式模型），而機率表示更傾向於獲得一個機率分佈（生成式模型）。

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）溫度和事故發生的機率如何建模？

咱們的目標在於模擬在不一樣溫度下事故發生的機率，在溫度和事故之間看起來沒有明顯的截止點。最終的問題是在溫度𝑡 時，事故發生的機率是多少。這個例子的目標就是要回答這個問題。

須要一個溫度的函數，稱之爲 p(t)，該函數的值域爲 [0，1]，這種函數有不少，可是根據觀察數據，一個可能比較適合的函數就是logistic函數。

𝛽 是這個模型的參數，是咱們須要從數據中進行推理獲得的。下面畫出𝛽 = 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：咱們的模型適用嗎？

目前爲止，咱們的機率推斷彷佛很完美，咱們完美地解釋了挑戰者號事故的發生。

可是有的讀者可能會問。本文選擇了 logistic 函數和先驗機率是合理的嗎。若是選擇了其它函數結果也許就不同了。如何可以知道咱們選擇的模型是好的？

這種懷疑徹底正確。考慮一種極端狀況，若是讓𝑝(𝑡) = 1, ∀𝑡，這保證事故必定會發生，咱們應該能夠再次預測在 1 月 28 號會發生災難。顯然這是一個不好的模型。

另外一方面，若是選擇了logistic 函數，並讓先驗機率都在 0 附近，可能會獲得不一樣的後驗機率。如何可以知道模型可以對數據進行了很好的表達？這就須要檢驗模型的吻合度問題。即度量模型的 擬合優度。
接下來的問題是，如何驗證咱們的模型的吻合度是否很差？

一種想法是將觀測數據和根據後驗參數估計獲得的模擬的數據集進行對比。這種作法的理論根據是，若是模擬出的數據和觀測數據在統計學上沒有類似的地方，那麼咱們的模型極可能沒有準確的表示觀測數據。

# -*- 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

在隨機模型中，能夠看到，隨着機率的增長，事故並無向右側彙集。常數模型也是這個狀況。

完美模型中沒有顯示機率曲線。由於它一直位於柱子的頂端。此處的完美模型只是用於演示，並不能從中推斷出什麼科學依據。

在實際項目中，咱們要評估的是咱們的後驗機率估計的分離圖是如何介於隨機模型和完美模擬之間的，天然，越靠近完美模型意味着咱們的模型擬合優度越好。

相關標籤/搜索

每日一句

每一个你不满意的现在，都有一个你没有努力的曾经。