概率分布及 PyMC 初探

概率分布

首先，讓我們回憶一下扮饶，什么是概率分布具练？設(shè)Z是一個(gè)隨機(jī)變量，則甜无，必定存在一個(gè)與Z有關(guān)的概率分布函數(shù)扛点。 給定Z的任意取值，我們都可以得到該Z對(duì)應(yīng)的概率值岂丘。

我們把隨機(jī)變量分為以下 3 種類(lèi)別：
Z 為離散的陵究。離散隨機(jī)變量的取值只能存在于某個(gè)特定的列表中。像人民幣（只有5毛奥帘、1元铜邮、20元、100元等面值）寨蹋、投票數(shù)等都屬于離散的隨機(jī)變量松蒜。
Z 為連續(xù)的。連續(xù)型隨機(jī)變量的值可以是任意精確度的數(shù)值已旧。像溫度和時(shí)間等都屬于連續(xù)型變量秸苗。因此對(duì)于這些數(shù)值，我們可以將其精確到任意程度运褪。連續(xù)型變量和離散型變量是一組相對(duì)的變量惊楼。
Z 為混合的〗斩铮混合型隨機(jī)變量的值檀咙，可以為以上兩種形式。也就是結(jié)合了上面兩種隨機(jī)變量形式的變量嗦枢。

離散情況

如果 Z 是離散的攀芯，則它的分布就是概率質(zhì)量函數(shù)。該函數(shù)度量的是文虏， Z 取值為 k 時(shí)的概率侣诺，即 P(Z = k)殖演。換句話說(shuō)，概率分布函數(shù)描述了隨機(jī)變量 Z年鸳。也就是說(shuō)趴久，如果知道了 Z 的概率質(zhì)量方程，我們就能夠完全的掌握 Z 的表現(xiàn)情況搔确。下面我們會(huì)介紹一些常見(jiàn)的概率質(zhì)量方程彼棍。

image.png

λ 被稱(chēng)為分布的一個(gè)參數(shù)，可以為任意正數(shù)膳算，它決定了我們這個(gè)分布的強(qiáng)度座硕。
k不同于λ， 可以為任意非負(fù)整數(shù)涕蜂，即k必須為 0,1,2 之類(lèi)的值勺鸦。這是非常重要的妆距，比如，當(dāng)我們?cè)谀M人口分布時(shí)，我們是不可以假設(shè)有 3.42 個(gè)或者 5.212 個(gè)人的肾档。
讓我們來(lái)使用一下這個(gè)分布函數(shù)妄壶，并且觀察一下λ對(duì)該概率函數(shù)的影響吊奢。這里我們會(huì)使用到上一個(gè)試驗(yàn)中用到的 scipy.stats 庫(kù)央勒，該庫(kù)已經(jīng)為我們定義了這些分布函數(shù)。讓我們先導(dǎo)入它們葱跋。

import scipy.stats as stats
from IPython.core.pylabtools import figsize
import numpy as np
from matplotlib import pyplot as plt
%matplotlib inline
# 定義 possion 函數(shù)
poi = stats.poisson
# 定義 兩個(gè) lambda 值
lambda_ = [1.5, 4.25]
# 定義 k值為0-15
k = np.arange(16)
k

接下來(lái)持寄，讓我們把 k 值傳入概率質(zhì)量函數(shù)中，得到每個(gè) k 對(duì)應(yīng)的概率值：

# 所有的分布函數(shù)類(lèi)中都存在 pmf 函數(shù)用于計(jì)算相應(yīng)的概率值
# 通過(guò)傳入不同的 labmda 值年局，求出兩組概率值
P_lambda_0 = poi.pmf(k, lambda_[0])
P_lambda_1 = poi.pmf(k, lambda_[1])
P_lambda_0, P_lambda_1

如上际看，我們已經(jīng)得到了每個(gè)k相應(yīng)的概率值。為了更好的了解λ 對(duì)于整個(gè)函數(shù)的影響矢否，這里仲闽，我們對(duì)我們求到的值進(jìn)行可視化。
首先僵朗，讓我們導(dǎo)入 Matplotlib 庫(kù)赖欣，并畫(huà)出相關(guān)圖像：

import matplotlib as mpl
from matplotlib import pyplot as plt
%matplotlib inline
# 設(shè)置所畫(huà)圖像的大小
figsize(12.5, 4)
# 將兩個(gè)lambda對(duì)應(yīng)的概率設(shè)置為兩種，不同的值
colours = ["#348ABD", "#A60628"]


plt.bar(k, P_lambda_0, color=colours[0],
        label="$\lambda = %.1f$" % lambda_[0], alpha=0.60,
        edgecolor=colours[0], lw="3")

plt.bar(k, P_lambda_1, color=colours[1],
        label="$\lambda = %.1f$" % lambda_[1], alpha=0.60,
        edgecolor=colours[1], lw="3")

plt.xticks(k + 0.4, k)
plt.legend()
plt.ylabel("probability of $k$")
plt.xlabel("$k$")
plt.title("Probability mass function of a Poisson random variable; differing \
$\lambda$ values")

從上圖中可以看到验庙，對(duì)于Poisson分布來(lái)說(shuō)顶吮，隨著λ 的增大，得到較大值的概率會(huì)增大粪薛。相反地悴了，當(dāng) λ 減小時(shí)，得到較小值的概率會(huì)增大。其次湃交，雖然x軸在 15 的時(shí)候就停止了熟空，但是分布并沒(méi)有在這里截止。他可以延伸到任意非負(fù)的整數(shù)搞莺。

image.png

image.png

image.png

a = np.linspace(0, 4, 100)
expo = stats.expon
lambda_ = [0.5, 1]

for l, c in zip(lambda_, colours):
    plt.plot(a, expo.pdf(a, scale=1./l), lw=3,
             color=c, label="$\lambda = %.1f$" % l)
    plt.fill_between(a, expo.pdf(a, scale=1./l), color=c, alpha=.33)

plt.legend()
plt.ylabel("PDF at $z$")
plt.xlabel("$z$")
# 設(shè)置y軸范圍
plt.ylim(0, 1.2)
plt.title("Probability density function of an Exponential random variable;\
 differing $\lambda$")

值得注意的是坤学，概率密度函數(shù)方程的某一個(gè)點(diǎn)的值并不等于它在這個(gè)點(diǎn)的概率疯坤。
什么是λ
這個(gè)問(wèn)題，我們可以理解為統(tǒng)計(jì)的動(dòng)機(jī)是什么深浮。在現(xiàn)實(shí)世界中压怠，我們并不知道λ 的存在，我們能直觀感受的就是變量Z飞苇。因此菌瘫，為了確定參數(shù)λ，我們就必須深入到整個(gè)事件的背景中去布卡。這個(gè)問(wèn)題雨让，其實(shí)很難，因?yàn)椴⒉淮嬖赯到λ 的對(duì)應(yīng)公式忿等。
對(duì)于λ的估計(jì)有很多設(shè)計(jì)好的方法栖忠，但因?yàn)棣瞬皇且粋€(gè)可以真正觀察到的東西。因此，誰(shuí)也不能說(shuō)哪一種方法是最好的庵寞。
貝葉斯推斷提出虚汛，就是為了對(duì)λ取值進(jìn)行估計(jì)。與其不斷的猜測(cè) λ的精確取值皇帮，不如用一個(gè)概率分布來(lái)描述λ 的可能取值卷哩。
這看起來(lái)或許有些奇怪。畢竟属拾，λ是一個(gè)定值将谊，它不一定是隨機(jī)的！我們?cè)趺茨軐?duì)一個(gè)非隨機(jī)變量值賦予一個(gè)概率呢渐白？不尊浓，這樣的思維方式其實(shí)是老舊的頻率派的思考方式。在貝葉斯的哲學(xué)體系下纯衍，我們可以通過(guò)所統(tǒng)計(jì)的數(shù)據(jù)栋齿，對(duì)λ賦予概率值（后面會(huì)詳細(xì)闡述）。因此襟诸，對(duì)參數(shù)λ估計(jì)是完全可以接受的瓦堵。
總結(jié)一下，貝葉斯推斷的主要思路歌亲，就是先給λ賦上一個(gè)先驗(yàn)分布菇用，然后找到我們能觀測(cè)到的數(shù)據(jù)和 λ之間的關(guān)系式，進(jìn)而建立模型陷揪。然后將我們觀測(cè)到的數(shù)據(jù)集放入模型中進(jìn)行訓(xùn)練惋鸥，最后得到λ的后驗(yàn)分布。

實(shí)例：從短信數(shù)據(jù)推斷行為

接下來(lái)悍缠，我們來(lái)模擬一個(gè)有趣的實(shí)例卦绣。這是一個(gè)關(guān)于用戶發(fā)送和收到短信的例子。
這里我們?yōu)槟闾峁┝艘粋€(gè)用戶收發(fā)短信條數(shù)的數(shù)據(jù)集合飞蚓，讓我們將它加載到本地：

!wget -nc "https://labfile.oss.aliyuncs.com/courses/1520/lab2-1.csv"

接下來(lái)滤港，我們還是利用 Python 先對(duì)這些數(shù)據(jù)進(jìn)行可視化，再來(lái)對(duì)其行為進(jìn)行推斷玷坠。

from matplotlib import pyplot as plt
from IPython.core.pylabtools import figsize
import numpy as np
%matplotlib inline
# 利用 numpy 加載數(shù)據(jù)蜗搔，該文件夾里存的就是每天短信的條數(shù) count_data
count_data = np.loadtxt("lab2-1.csv")
n_count_data = len(count_data)
n_count_data, count_data

根據(jù)上面的結(jié)果，該文件共存儲(chǔ)了該用戶 74 天內(nèi)發(fā)送和收到的短信條數(shù)八堡。接下來(lái)樟凄，讓我們來(lái)觀察用戶的短信使用行為是否隨著時(shí)間有所改變。短信的條數(shù)是循序漸進(jìn)兄渺？還是突然的變化缝龄？

figsize(12.5, 3.5)
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)

從上圖中，我們可以很好的發(fā)現(xiàn)，用戶在后半段時(shí)間收發(fā)短信的條數(shù)明顯變多叔壤，即他的收發(fā)短信行為可能發(fā)生了改變瞎饲。但是，從圖中炼绘，我們很難判斷出他到底是何時(shí)發(fā)生的行為改變嗅战。因此，我們需要利用貝葉斯推斷進(jìn)行計(jì)算俺亮。

image.png

image.png

image.png

import pymc3 as pm
# 定義一個(gè)變量模型幌氮，后面的所有隨機(jī)變量都會(huì)寫(xiě)入其中
model = pm.Model()

接下來(lái)缭受，讓我們根據(jù)公式，將(τ,λ1,λ2)等隨機(jī)變量寫(xiě)入 model 中该互。

with model:  # 利用with定義一個(gè)上下文管理器米者，以 model 作為上下文。
           # 在這個(gè)上下文中定義的變量都會(huì)被添加到這個(gè)模型中去

    # α 變量為樣本中計(jì)算平均值的逆（即倒數(shù)的意思）
    alpha = 1.0/count_data.mean()

    # 兩個(gè) λ的值都服從指數(shù)密度函數(shù)分布宇智，函數(shù)參數(shù)為α
    lambda_1 = pm.Exponential("lambda_1", alpha)
    lambda_2 = pm.Exponential("lambda_2", alpha)

    # τ蔓搞，即λ發(fā)生改變的時(shí)間，為 0 - 73 中的任意一天(總共74天)
    # 因?yàn)棣邮谴淼奶鞌?shù)普筹，是離散的败明，這里通過(guò) DiscreteUniform 設(shè)置它為離散型變量
    # 即產(chǎn)生隨機(jī)數(shù)
    tau = pm.DiscreteUniform("tau", lower=0, upper=n_count_data - 1)
model

在上面代碼中，我們創(chuàng)建了對(duì)應(yīng)于參數(shù)λ1和λ2的 PyMC 變量太防。并且設(shè)置它們?yōu)?PyMC 的隨機(jī)變量，這些變量服從于它們自己的分布函數(shù)酸员。

image.png

我們用變量 lambda_ 來(lái)存每個(gè)時(shí)刻的λ值蜒车，代碼如下：

with model:
    idx = np.arange(n_count_data)  # idf 表示天數(shù)
    # lambda_ 表示每天的 λ的值
    # 當(dāng)小于τ 時(shí)為lambda_1 ，大于為lambda_2
    lambda_ = pm.math.switch(tau >= idx, lambda_1, lambda_2)
type(lambda_)

這里是無(wú)法將 lambda_ 的具體的值顯示出來(lái)的幔嗦。因?yàn)?lambda_1酿愧、lambda_2、tau 是隨機(jī)的邀泉，所以 lambda_ 也會(huì)是隨機(jī)的嬉挡。它們只有在模型編譯時(shí)钝鸽，才會(huì)出現(xiàn)具體的值。目前我們只是定義了它們的分布方式庞钢，還未進(jìn)行具體計(jì)算拔恰。
讓我們總結(jié)一下思路，我們觀察了我們得到的用戶收發(fā)短信的數(shù)據(jù)基括，并且假設(shè)了這些數(shù)據(jù)服從于 Poisson 分布颜懊。而對(duì)于 Poisson 分布，最關(guān)鍵的便是λ (即變量 lambda_ )的求解风皿。為了更加準(zhǔn)確的求取 lambda_ 河爹，根據(jù)圖像，我們將 lambda_ 分成了兩個(gè)階段桐款，分別為 lambda_1,lambda_2 這兩個(gè)值咸这。再加上參數(shù)變化的時(shí)刻 tau。整個(gè)模型需要求解的便是 (lambda_1魔眨，lambda_2媳维，tau) 這三個(gè)變量。
接下來(lái)我們就來(lái)將定義模型與實(shí)際統(tǒng)計(jì)的數(shù)據(jù)相結(jié)合冰沙。

with model:
    # 將統(tǒng)計(jì)數(shù)據(jù) 與 參數(shù)為 lambda_ 的 Poisson 分布相結(jié)合侨艾。
    observation = pm.Poisson("obs", lambda_, observed=count_data)
type(observation)

我們假設(shè)的觀察數(shù)據(jù)是通過(guò) Poisson 分布產(chǎn)生，因此我們需要把這些參數(shù)和 poisson 分布結(jié)合拓挥，得到觀測(cè)數(shù)據(jù)的產(chǎn)生模型唠梨。然后將真實(shí)的觀測(cè)數(shù)據(jù)放入模型中進(jìn)行訓(xùn)練，進(jìn)而得到那幾個(gè)參數(shù)的具體后驗(yàn)侥啤。
目標(biāo)數(shù)據(jù)已經(jīng)傳入了 model当叭，需要求解的參數(shù)也已經(jīng)傳入了 model。接下來(lái)盖灸，其實(shí)我們只需要按下“學(xué)習(xí)”的按鈕蚁鳖，讓計(jì)算機(jī)自行學(xué)習(xí)，就可以得到參數(shù) (τ,λ1,λ2)的值了赁炎。
接下來(lái)醉箕，我們就來(lái)編寫(xiě)按下這個(gè)“按鈕”的代碼。但是這里會(huì)使用到一種叫做 馬爾科夫鏈蒙特卡洛(MCMC) 的復(fù)雜理論徙垫。這個(gè)復(fù)雜理論和代碼將在后面的實(shí)驗(yàn)進(jìn)行闡述讥裤。為了展示結(jié)果，這里姻报，我們只需要運(yùn)行下面代碼即可己英。無(wú)需自己手動(dòng)敲。

# 下面代碼將在之后的實(shí)驗(yàn)中吴旋，進(jìn)行闡述
with model:
    step = pm.Metropolis()
    trace = pm.sample(10000, tune=5000, step=step)

通過(guò)上面的學(xué)習(xí)损肛，我們從λ1厢破、λ2和τ各自所對(duì)應(yīng)的后驗(yàn)分布函數(shù)中，得到了數(shù)千個(gè)隨機(jī)樣本治拿。樣本如下摩泪，這些變量都是通過(guò)學(xué)習(xí)到的各自的分布函數(shù)所產(chǎn)生。因?yàn)檫@些變量的具體后驗(yàn)分布公式很能表示忍啤，因此返回的是由后驗(yàn)分布產(chǎn)生的隨機(jī)樣本加勤。

lambda_1_samples = trace['lambda_1']
lambda_2_samples = trace['lambda_2']
tau_samples = trace['tau']
lambda_1_samples, lambda_2_samples, tau_samples

我們只需要將這些值進(jìn)行統(tǒng)計(jì)即可。例如同波，計(jì)算 lambda_1 的樣例中每個(gè)數(shù)字所出現(xiàn)的頻率鳄梅。通過(guò)這些頻率，我們就能夠畫(huà)出每個(gè)參數(shù)的后驗(yàn)分布的樣子未檩。下面戴尸，則是我們通過(guò)直方圖，來(lái)展示每種參數(shù)的后驗(yàn)概率冤狡。

figsize(12.5, 10)
# 下面的代碼全為圖像可視化代碼
# 將圖像分為三份孙蒙，現(xiàn)在ax代表的是第一行第一張
ax = plt.subplot(311)
ax.set_autoscaley_on(False)
# 畫(huà)出lambda_1 的后驗(yàn)分布
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""" $\lambda_1,\;\lambda_2,\;\tau$ """)
plt.xlim([15, 30])
plt.xlabel("$\lambda_1$ value")

# 畫(huà)出lambda_2 的后驗(yàn)分布
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")

# 畫(huà)出lambda_3的后驗(yàn)分布
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$ (day)")
plt.ylabel("probability")

說(shuō)明
回想一下，貝葉斯方法返回一個(gè)分布悲雳。因此挎峦，現(xiàn)在我們可以用分布來(lái)描述未知的λ和τ。我們也可以從上面的分布中合瓢，找到參數(shù)的合理值：λ1大概為 18 和λ2大概為23坦胶。這兩個(gè)\lambdaλ 的后驗(yàn)分布明顯不同，這表明用戶接受短信的行為確實(shí)發(fā)生了變化晴楔。
我們還可以注意到λ的后驗(yàn)分布并不像是指數(shù)分布顿苇。事實(shí)上，后驗(yàn)分布并不是我們可以從原始模型中辨別的任何分布税弃。正因如此纪岁，我們利用貝葉斯推斷出來(lái)的不是一個(gè)分布的函數(shù)值，而是一批分布中的樣本點(diǎn)则果。如果我們改用數(shù)學(xué)方式處理這個(gè)問(wèn)題幔翰，它就會(huì)變得棘手和混亂。
再來(lái)看看τ的一個(gè)分布西壮。由于它是一個(gè)離散的變量导匣，所以它的后驗(yàn)看起來(lái)和其他的兩個(gè)參數(shù)有點(diǎn)不同，它不存在概率區(qū)間茸时。我們可以看到在第3天，有50%的把握可以說(shuō)用戶的行為是有所改變的赋访。并且可都，從上面結(jié)果可以看出缓待，僅僅有三到四天可以認(rèn)為是潛在的轉(zhuǎn)折點(diǎn)。
后驗(yàn)樣本的作用
在本課程的后面章節(jié)渠牲，我們會(huì)著重闡述這個(gè)問(wèn)題⌒矗現(xiàn)在我們用另外一個(gè)實(shí)例對(duì)它進(jìn)行簡(jiǎn)單的闡述。
我們需要用我們剛才得到的后驗(yàn)樣本來(lái)回答下面問(wèn)題：
在第t(0≤t≤73)天中签杈，期望收發(fā)所少條短信呢瘫镇？
由于Poisson分布的期望值等于它的參數(shù)λ。因此問(wèn)題相當(dāng)于：在時(shí)間 t中答姥，參數(shù)λ的期望值是多少铣除。

image.png

N = tau_samples.shape[0]
expected_texts_per_day = np.zeros(n_count_data)
for day in range(0, n_count_data):
    ix = day < tau_samples
    # 求第 t 天的期望尚粘，即所有可能的 lambda 的平均值
    # 當(dāng)t<tau時(shí)，所有l(wèi)ambda2 會(huì)為0敲长，則不會(huì)將 lambda_2_samples計(jì)算進(jìn)去
    # lambda_1_samples[ix].sum() 和 lambda_2_samples[~ix].sum() 是互斥的
    expected_texts_per_day[day] = (lambda_1_samples[ix].sum()
                                   + lambda_2_samples[~ix].sum()) / N
expected_texts_per_day

最后郎嫁，我們將這些期望展示到原圖上，觀察它與實(shí)際收發(fā)信息的關(guān)系祈噪。

figsize(12.5, 5)
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")

上圖的結(jié)果很明顯的說(shuō)明了轉(zhuǎn)折點(diǎn)的重要影響泽铛。我們的分析結(jié)果非常符合之前的估計(jì)：用戶行為確實(shí)發(fā)生了改變，而且這是一個(gè)突然的變化辑鲤，而非循序漸進(jìn)的變化盔腔。我們可以推測(cè)這種情況的產(chǎn)生原因是：短信費(fèi)用降低，天氣提醒短信的訂閱或者一段新的感情等遂填。
兩個(gè)λ值是否真的一樣
在短信接收例子中铲觉，我們只管地觀測(cè)了λ1和 λ2的先驗(yàn)信息。并認(rèn)為他們是不同的吓坚。這很公平撵幽，畢竟先驗(yàn)的位置基本離得非常遠(yuǎn)。但是這種觀察的主管因素過(guò)強(qiáng)礁击，下面盐杂，我們將會(huì)介紹一種較為正式的方法，來(lái)判斷兩個(gè) λ 值是否相等哆窿。
一種方法是計(jì)算出 P(λ1<λ2∣data)链烈。即在獲得觀察數(shù)據(jù)的前提下，λ1的真實(shí)值比λ2小的概率挚躯。如果這個(gè)概率接近 50% ,那就相當(dāng)于是投硬幣的結(jié)果强衡，我們則不能確定這兩個(gè)λ 值是不同的。如果码荔，這個(gè)概率接近 100% ,那么我們就可以很自信的說(shuō)漩勤，這兩個(gè)值必定不同感挥。
利用上我們得到的后驗(yàn)樣本，這種計(jì)算是非常簡(jiǎn)單的越败。我們只需計(jì)算λ1的后驗(yàn)樣本中比λ2的后驗(yàn)樣本小的次數(shù)的比例即可触幼。

print((lambda_1_samples < lambda_2_samples).mean())

結(jié)果很顯然，現(xiàn)在我們幾乎 100% 的把握可以說(shuō)這兩個(gè)值是不等的究飞。
當(dāng)然置谦，我們還可以問(wèn)的更加詳細(xì)一點(diǎn)。比如亿傅，兩個(gè)值之間相差至少 1媒峡、2、5袱蜡、10 的概率有多大丝蹭？

for d in [1, 2, 5, 10]:
    v = (abs(lambda_1_samples-lambda_2_samples) >= d).mean()
    print("兩個(gè) lambda 之間相差至少{}的概率為 : %{}".format(d, v*100))

多行為變化推斷

數(shù)據(jù)集的下載

同樣，讓我們先從藍(lán)橋云課的網(wǎng)站上下載這些數(shù)據(jù)集：

!wget -nc "https://labfile.oss.aliyuncs.com/courses/1520/lab2-1.csv"

接下來(lái)坪蚁，從這個(gè)數(shù)據(jù)集合讀入到模型中：

import scipy.stats as stats
from IPython.core.pylabtools import figsize
import numpy as np
from matplotlib import pyplot as plt
%matplotlib inline

count_data = np.loadtxt("lab2-1.csv")
n_count_data = len(count_data)
n_count_data, count_data

模型的建立

image.png

挑戰(zhàn)：將上面幾個(gè)參數(shù)加入我們定義的模型中

import pymc3 as pm
import theano.tensor as tt
#定義一個(gè)變量模型奔穿，后面的所有隨機(jī)變量都會(huì)寫(xiě)入其中
model = pm.Model()
with model:#利用with定義一個(gè)上下文管理器，以 model 作為上下文敏晤。
           #在這個(gè)上下文中定義的變量都會(huì)被添加到這個(gè)模型中去
        
    # α 變量為樣本中計(jì)算平均值的逆（即倒數(shù)的意思）
    alpha = 1.0/count_data.mean()  
    
    #三個(gè) λ的值都服從指數(shù)密度函數(shù)分布贱田，函數(shù)參數(shù)為α
    lambda_1 = pm.Exponential("lambda_1", alpha)
    lambda_2 = pm.Exponential("lambda_2", alpha)
    lambda_3 = pm.Exponential("lambda_3", alpha)
                              
    #τ，即λ發(fā)生改變的時(shí)間嘴脾，為 0 -70 中的任意一天(總共71天)
    #因?yàn)棣邮谴淼奶鞌?shù)男摧，是離散的，這里通過(guò) DiscreteUniform 設(shè)置它為離散型變量
    #即產(chǎn)生隨機(jī)數(shù)
    tau_1 = pm.DiscreteUniform("tau_1", lower=0, upper=n_count_data - 1)
    tau_2 = pm.DiscreteUniform("tau_2", lower=0, upper=n_count_data - 1)
model

挑戰(zhàn)：完成分段函數(shù) lambda_ 的賦值译打。

with model:
    idx = np.arange(n_count_data) # idf 表示天數(shù)
    #lambda_ 表示每天的 λ的值
    #使用where 條件語(yǔ)句去達(dá)到 3 分段函數(shù)賦值的效果
    lambda_ =  pm.math.where (idx<tau_1,lambda_1, 
                              pm.math.where (idx<tau_2,lambda_2,lambda_3))
type(lambda_)  

接下來(lái)耗拓，就是將這我們定義的參數(shù)和實(shí)際數(shù)據(jù)相結(jié)合，代碼和上一個(gè)試驗(yàn)相同：

with model:
    # 將統(tǒng)計(jì)數(shù)據(jù) 與 參數(shù)為 lambda_ 的 Poisson 分布相結(jié)合奏司。
    observation = pm.Poisson("obs", lambda_, observed=count_data)
type(observation)

最后就是摁下“學(xué)習(xí)”的按鈕乔询，讓模型進(jìn)行訓(xùn)練與學(xué)習(xí)：

with model:
    step = pm.Metropolis()
    trace = pm.sample(1000, tune=500, step=step)

通過(guò) trace 函數(shù)得到這些參數(shù)的后驗(yàn)分布函數(shù)下的幾千個(gè)隨機(jī)變量：

lambda_1_samples = trace['lambda_1']
lambda_2_samples = trace['lambda_2']
lambda_3_samples = trace['lambda_3']
tau1_samples = trace['tau_1']
tau2_samples = trace['tau_2']
print("lambda1:", lambda_1_samples)
print("lambda2:", lambda_2_samples)
print("lambda3:", lambda_3_samples)
print("tau_1:", tau1_samples)
print("tau_2:", tau2_samples)

數(shù)據(jù)可視化

挑戰(zhàn)：分別畫(huà)出以上 5 個(gè)參數(shù)的分布圖像

figsize(12.5, 10)
#下面的代碼全為圖像可視化代碼
#將圖像分為三份，現(xiàn)在ax代表的是第一行第一張
ax = plt.subplot(511)
ax.set_autoscaley_on(False)
#畫(huà)出lambda_1 的后驗(yàn)分布
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")

#畫(huà)出lambda_2 的后驗(yàn)分布
ax = plt.subplot(512)
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([50, 65])
plt.xlabel("$\lambda_2$ value")

#畫(huà)出lambda_3 的后驗(yàn)分布
ax = plt.subplot(513)
ax.set_autoscaley_on(False)
plt.hist(lambda_3_samples, histtype='stepfilled', bins=30, alpha=0.85,
         label="posterior of $\lambda_3$ ", color="#7A68A6", normed=True)
plt.legend(loc="upper left")
plt.xlim([15, 30])
plt.xlabel("$\lambda_2$ value")

#畫(huà)出 tau_1 的后驗(yàn)分布
plt.subplot(514)
w = 1.0 / tau1_samples.shape[0] * np.ones_like(tau1_samples)
plt.hist(tau1_samples, bins=n_count_data, alpha=1,
         label=r" posterior of $\tau_1$ ",
         color="#467821", weights=w, rwidth=2.)
plt.xticks(np.arange(n_count_data))

plt.legend(loc="upper left")
plt.ylim([0, .75])
plt.xlim([20, 35])
plt.xlabel(r"$\tau$ (day)")
plt.ylabel("probability");

##畫(huà)出 tau_2 的后驗(yàn)分布
plt.subplot(515)
w = 1.0 / tau2_samples.shape[0] * np.ones_like(tau2_samples)
plt.hist(tau2_samples, bins=n_count_data, alpha=1,
         label=r"posterior of $\tau_2$ ",
         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$ (day)")
plt.ylabel("probability");