黑客級概率程序設計和貝葉斯方法

揭開MCMC的神秘面紗

前面兩章對讀者隱藏了PyMC的內部機制蚕苇，也就是Markov chain Monte Carlo（MCMC）童社。我們引入本章目的有三翘地。第一點涡上，任何關于貝葉斯推斷的書籍肯定會討論MCMC猩系。我無法抗拒此點。要抱怨就去找那些統(tǒng)計學家吧兆蕉。：）第二點羽戒，知道MCMC的過程會給你關于算法收斂的理解。（收斂什么虎韵？等會我們再介紹易稠。）第三點，我們會弄清楚為何我們可以從后驗分布中產生的上千個樣本作為最終的結果包蓝，而這些剛開始我們會覺得相當困惑的驶社。

貝葉斯理論地貌（landscape）

landscape 不知道怎么翻譯，太 tm 糾結了养晋。

當我們開始一個包含N個未知量的貝葉斯推斷問題衬吆，我們隱式地創(chuàng)建了一個N維空間來包含先驗分布。與這個空間關聯的是一個額外的維度绳泉，這個位于空間之上的東西被稱作表面surface或者曲線curve，該維度也反應了一個特定點的先驗概率姆泻×憷遥空間上的表面由我們的先驗分布定義。例如拇勃，若我們有兩個未知量p_1和p_2四苇，兩者的先驗分布都是Uniform(0, 5)，那么創(chuàng)建出來的空間就是一個長度為5的正方形方咆，而表面就是一個扁平的平面月腋，位于這個正方形紙上，代表了每個點瓣赂。

%matplotlib inline
import scipy.stats as stats
from IPython.core.pylabtools import figsize
import numpy as np
figsize(12.5, 4)

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

jet = plt.cm.jet
fig = plt.figure()
x = y = np.linspace(0, 5, 100)
X, Y = np.meshgrid(x, y)

plt.subplot(121)
uni_x = stats.uniform.pdf(x, loc=0, scale=5)
uni_y = stats.uniform.pdf(y, loc=0, scale=5)
M = np.dot(uni_x[:, None], uni_y[None, :])
im = plt.imshow(M, interpolation='none', origin='lower',
                cmap=jet, vmax=1, vmin=-.15, extent=(0, 5, 0, 5))

plt.xlim(0, 5)
plt.ylim(0, 5)
plt.title("Landscape formed by Uniform priors.")

ax = fig.add_subplot(122, projection='3d')
ax.plot_surface(X, Y, M, cmap=plt.cm.jet, vmax=1, vmin=-.15)
ax.view_init(azim=390)
plt.title("Uniform prior landscape; alternate view");

Uniform prior.png

換個例子榆骚，如果兩個先驗分布是Exp(3)和Exp(10)，那么空間是在2-D平面上的所有正數煌集，并且由兩個先驗分布產生的表面看起來就像是瀑布一樣從(0, 0)開始妓肢，流經所有的正數。

下面的圖將這個場景圖形化出來了苫纤。顏色越是深紅碉钠，先驗概率就分配得越高纲缓。相反地，越是深藍則表示分配的概率越小喊废。

figsize(12.5, 5)
fig = plt.figure()
plt.subplot(121)

exp_x = stats.expon.pdf(x, scale=3)
exp_y = stats.expon.pdf(x, scale=10)
M = np.dot(exp_x[:, None], exp_y[None, :])
CS = plt.contour(X, Y, M)
im = plt.imshow(M, interpolation='none', origin='lower',
                cmap=jet, extent=(0, 5, 0, 5))
#plt.xlabel("prior on $p_1$")
#plt.ylabel("prior on $p_2$")
plt.title("$Exp(3), Exp(10)$ prior landscape")

ax = fig.add_subplot(122, projection='3d')
ax.plot_surface(X, Y, M, cmap=jet)
ax.view_init(azim=390)
plt.title("$Exp(3), Exp(10)$ prior landscape; \nalternate view")

exponential prior.png

這些例子都是很簡單的2D空間上的祝高，我們的大腦還是能夠理解這些表面的。實際應用中污筷，由先驗分布產生出來的空間和表面都是很高維的工闺。

如果這些表面描述了關于未知量的先驗分布，那么在我們獲得了觀測值X之后颓屑，這個空間將會發(fā)生怎樣的變化呢斤寂？數據不會對空間產生變化，但是卻會對空間的表面產生影響揪惦，拉伸或者壓縮先驗表面的特性來反應哪些地方是真的參數遍搞。更多的數據意味著更多的拉伸和壓縮，原始的形狀會和最終形成的形狀有一定的差異器腋，有時候會完全被顛覆掉溪猿。少的數據，則原始形狀和最終結果的形狀差不多纫塌。不管怎樣诊县，最終的表面描述了后驗分布。

不幸的是措左，我必須強調的是依痊，高維空間不可能可視化。對2D空間怎披，數據常常會上推表面來產生高的山峰胸嘁。觀察道的數據在某些區(qū)域上推后驗概率的傾向（tendency）由先驗概率分布檢查院溺，使得低先驗概率意味著更多的阻抗（resistance）瘟芝。因此，在上面的雙指數先驗例子中玛痊，山峰可能在靠近(0, 0)爆發(fā)得比靠近(5, 5)處爆發(fā)得更高状飞，因為在(5, 5)處附近有更多的阻抗毫胜。這個峰值反映了后驗概率，哪里真的參數更可能被發(fā)現诬辈。需要重視的是酵使，如果先驗已經被分配概率為0，那么沒有后驗概率被分配在那個地方自晰。

假設上面提到的先驗概率表示不同的參數\lambda給定的泊松分布凝化。我們觀察到一些數據點并可視化新的地貌（landscape）：

# create the observed data

# sample size of data we observe, trying varying this (keep it less than 100 ;)
N = 1

# the true parameters, but of course we do not see these values...
lambda_1_true = 1
lambda_2_true = 3

#...we see the data generated, dependent on the above two values.
data = np.concatenate([
    stats.poisson.rvs(lambda_1_true, size=(N, 1)),
    stats.poisson.rvs(lambda_2_true, size=(N, 1))
], axis=1)
print "observed (2-dimensional,sample size = %d):" % N, data

# plotting details.
x = y = np.linspace(.01, 5, 100)
likelihood_x = np.array([stats.poisson.pmf(data[:, 0], _x)
                        for _x in x]).prod(axis=1)
likelihood_y = np.array([stats.poisson.pmf(data[:, 1], _y)
                        for _y in y]).prod(axis=1)
L = np.dot(likelihood_x[:, None], likelihood_y[None, :])

觀察到（2D，樣本大谐贶瘛＝1）：[[2 2]]

figsize(12.5, 12)
# matplotlib heavy lifting below, beware!
plt.subplot(221)
uni_x = stats.uniform.pdf(x, loc=0, scale=5)
uni_y = stats.uniform.pdf(x, loc=0, scale=5)
M = np.dot(uni_x[:, None], uni_y[None, :])
im = plt.imshow(M, interpolation='none', origin='lower',
                cmap=jet, vmax=1, vmin=-.15, extent=(0, 5, 0, 5))
plt.scatter(lambda_2_true, lambda_1_true, c="k", s=50, edgecolor="none")
plt.xlim(0, 5)
plt.ylim(0, 5)
plt.title("Landscape formed by Uniform priors on $p_1, p_2$.")

plt.subplot(223)
plt.contour(x, y, M * L)
im = plt.imshow(M * L, interpolation='none', origin='lower',
                cmap=jet, extent=(0, 5, 0, 5))
plt.title("Landscape warped by %d data observation;\n Uniform priors on $p_1, p_2$." % N)
plt.scatter(lambda_2_true, lambda_1_true, c="k", s=50, edgecolor="none")
plt.xlim(0, 5)
plt.ylim(0, 5)

plt.subplot(222)
exp_x = stats.expon.pdf(x, loc=0, scale=3)
exp_y = stats.expon.pdf(x, loc=0, scale=10)
M = np.dot(exp_x[:, None], exp_y[None, :])

plt.contour(x, y, M)
im = plt.imshow(M, interpolation='none', origin='lower',
                cmap=jet, extent=(0, 5, 0, 5))
plt.scatter(lambda_2_true, lambda_1_true, c="k", s=50, edgecolor="none")
plt.xlim(0, 5)
plt.ylim(0, 5)
plt.title("Landscape formed by Exponential priors on $p_1, p_2$.")

plt.subplot(224)
# This is the likelihood times prior, that results in the posterior.
plt.contour(x, y, M * L)
im = plt.imshow(M * L, interpolation='none', origin='lower',
                cmap=jet, extent=(0, 5, 0, 5))

plt.scatter(lambda_2_true, lambda_1_true, c="k", s=50, edgecolor="none")
plt.title("Landscape warped by %d data observation;\n Exponential priors on \
$p_1, p_2$." % N)
plt.xlim(0, 5)
plt.ylim(0, 5)

landscape.png

還有要注意的是，那些對應于深紅顏色的地方的那些最高的點倦畅，在指數情形下更加偏向在(0, 0)角落遮糖，這里是指數先驗分布給出更多權重的地方。黑色的點告訴我們真實的參數叠赐，即使有一個樣本點欲账，山峰嘗試包含真實參數。當然了芭概，以樣本大小為1進行的推斷是相當簡易的赛不，選擇這樣一個小的樣本作為例子只是為了更好地解釋原理。

當我們嘗試更大的樣本集合罢洲，可以觀測到山峰后驗分布如何改變踢故。

使用MCMC來探索地貌

我們應該去研究變形后的后驗空間來尋找后驗山峰。然而惹苗，我們不能夠天真地直接對空間進行搜索：任何計算機科學家都會告訴你遍歷N維空間是N的指數級困難的：空間的大小會在N增長時畴椰，發(fā)生快速的膨脹（維度詛咒）。我們如何找尋這些隱藏的山峰呢鸽粉？MCMC背后的思想是執(zhí)行一個空間上的智能搜索。這里的search指的是我們在找尋一個特定的點抓艳，這個可能b

MCMC從后驗分布中返回樣本而不是分布本身触机。MCMC執(zhí)行了一種類似于詢問“我發(fā)現的這個pebble來自我尋找的山峰的可能性有多大？”并且完成這個任務通過返回可以接受的pebble來重新創(chuàng)造原始的山峰玷或。用MCMC和PyMC的語言就是儡首，返回的pebble序列也就是那些樣本，合起來稱為軌跡(traces)偏友。

當我說MCMC智能地搜索時蔬胯，我實際上在說MCMC會有希望收斂到高后驗概率的地方。MCMC通過探索周圍的位置并移動到具有更高概率的地方位他。并且氛濒，收斂不是一個精確的空間产场，但是MCMC移動到一個空間中更加廣闊的地方并且在那里進行隨機游走，從那片地方中選出樣本舞竿。

為何需要成千個樣本京景？

首先，返回成千個樣本給用戶可能聽起來不是一種可以高效地描述后驗分布的方式骗奖。我將論證這種方式是特別高效的确徙。考慮下面替代的方式：

1. 返回一個對山峰范圍一個數學公式执桌，描述了帶有任意峰頂和山谷的N維表面
2. 返回地貌的峰頂鄙皇，數學上可能和合理的東西往往是最高點，對應于最可能的對未知量的估計仰挣，忽略地貌的形狀伴逸，而這個東西是我們之前認為是相當重要的可以確定后驗置信度。

除了計算方面的原因椎木，返回樣本的很可能最強的原因是我們可以很容易地使用大數定律來解決難解問題违柏。我們將這個討論推后到下一章節(jié)。擁有了上千個樣本香椎，我們可以將他們放在一個直方圖中重新構造后驗分布的表面漱竖。

執(zhí)行MCMC的算法

有一大群算法都可以完成MCMC。大多數算法粗略的過程如下：

1. 從當前位置開始
2. 提議一個到新位置的移動（找到靠近你的pebble）
3. 基于這個位置對數據的信任度和先驗分布（詢問pebble是否可能來自山峰）
4. 如果你接受：移動到新的位置畜伐，返回step 1馍惹；否則不移動，返回step 1.
5. 在大量迭代后玛界，返回所有接受的位置

這個方式讓我們在宏觀上可以移動到后驗分布存在的區(qū)域万矾，并且謹慎地收集到在這個過程中產生的樣本。
如果MCMC算法當前的位置是在一個極其低的概率的地方慎框，常常是當算法開始的時候（典型滴就是空間中隨機獲得一個位置）良狈，算法會移動到那些不太可能從后驗分布中產生出來的位置，不過比周圍的地方要略好笨枯。因此算法首先的移動不會對后驗分布產生影響薪丁。
在上面的偽代碼中，只有當前位置起到作用馅精。我們可以將這種特性稱為無記憶性严嗜，算法不關注如何到達當前位置的，只關心在當前位置上是什么洲敢。

其余的一些求解后驗分布的近似解

除了MCMC漫玄，還有一些方法可以確定后驗分布。Laplace近似使用簡單的函數來近似后驗分布压彭。更加高級的方法叫做變分貝葉斯法睦优。這三種方法都有各自的優(yōu)缺點渗常。本書中我們將只關注MCMC法。

That being said, my friend Imri Sofar likes to classify MCMC algorithms as either "they suck", or "they really suck". He classifies the particular flavour of MCMC used by PyMC as just sucks ;)

例子：使用mixture模型來進行無監(jiān)督聚類

假設我們有下面的數據集：

figsize(12.5, 4)
data = np.loadtxt("data/mixture_data.csv", delimiter=",")

plt.hist(data, bins=20, color="k", histtype="stepfilled", alpha=0.8)
plt.title("Histogram of the dataset")
plt.ylim([0, None])
print data[:10], "..."

[ 115.85679142  152.26153716  178.87449059  162.93500815  107.02820697
  105.19141146  118.38288501  125.3769803   102.88054011  206.71326136] ...

example.png

這個數據暗示了什么刨秆？看起來是一個雙模型凳谦，有兩個峰值，一個接近120衡未，另一個200尸执。可能這個數據集合有兩個聚類了缓醋。

這個數據集是上一章里面的數據生成技術的很好的例子如失。我們可以假設數據是如何被產生出來的。下面就是數據產生算法：

1. 對每個數據點送粱，以概率p選擇cluster 1褪贵，否則就選擇cluster 2.
2. 然后按照1中選出來的cluster i，從對應的參數為\mu_i和\sigma_i的正態(tài)分布中產生一個隨機值
3. 重復這個步驟抗俄，知道所有的點都被產生脆丁。

這個算法可以創(chuàng)建一個與真實觀測數據集合相似的數據集合，所以我們選擇它作為自己的模型动雹。當然槽卫，我們并不知道p或者正態(tài)分布的參數。因此我們肯定要推斷胰蝠，或者學習這些未知量歼培。

將正態(tài)分布表示為Nor_0和Nor_1（從0開始就是pythonic而已）。兩個都是不知道對應的均值和標準方差茸塞，分別表示為\mu_i和\sigma_i躲庄，i=0,1。特定的數據點可以從Nor_0或者Nor_1中產生钾虐，我們假設數據點以概率p分配給Nor_0噪窘。

使用一個PyMC Categorical隨機變量可以用來將數據點分給相應的cluster。它的參數是一個k長度的概率的數組效扫，其中項的和應該為1效览，它的值是在0到k-1之間的整數，取的概率就是對應于數組中的概率荡短。在我們這里的情形k=2。先驗地哆键，我們并不知道分配給cluster 1的概率是多少掘托，所以我們創(chuàng)建了一個0，1上的均勻分布變量對此進行建模籍嘹。設為p闪盔。因此這個輸入到Categorical變量的概率數組就是[p, 1-p]弯院。

import pymc as pm

p = pm.Uniform("p", 0, 1)

assignment = pm.Categorical("assignment", [p, 1 - p], size=data.shape[0])
print "prior assignment, with p = %.2f:" % p.value
print assignment.value[:10], "..."

觀察一下上面的數據集，我們猜測這兩個正態(tài)分布的標準方差都是不同的泪掀。為了保持對標準方差的未知听绳，我們首先會將他們建模成0~100上的均勻分布。實際上我們在談論的是\tau异赫，也就是正態(tài)分布的準確性椅挣，但是可以從標準方差的角度來衡量這樣的準確性。我們PyMC的代碼需要轉換我們的標準方差稱準確性塔拳，他們有下面的關系：

\tau = \frac{1}{\sigma^2}

在PyMC中鼠证，我們可以使用下面這句代碼完成轉換：

\tau = 1.0 / pm.Uniform( "stds", 0, 100, size=2 ) ** 2

注意我們將size設置為2：我們將這兩個\tau均建模稱一個單個的PyMC變量。注意這里不包含兩個\tau之間的關系靠抑。

我們同樣需要指定cluster的中心的先驗量九。中心實際上就是正態(tài)分布的參數\mu。他們的先驗可以建模成一個正態(tài)分布颂碧≤校看一下數據，我猜測這兩個中心分別是120和190载城，盡管我們對這兩個肉眼得到的結果并不是很確定肌似。因此我將\mu_0=120, \mu_1=190,\sigma_{0,1}=10（然后在PyMC變量中的\tau = frac{1}{\sigma^2} = 0.01）

taus = 1.0 / pm.Uniform("stds", 0, 100, size=2) ** 2
centers = pm.Normal("centers", [120, 190], [0.01, 0.01], size=2)

"""
The below deterministic functions map an assignment, in this case 0 or 1,
to a set of parameters, located in the (1,2) arrays `taus` and `centers`.
"""

@pm.deterministic
def center_i(assignment=assignment, centers=centers):
    return centers[assignment]

@pm.deterministic
def tau_i(assignment=assignment, taus=taus):
    return taus[assignment]

print "Random assignments: ", assignment.value[:4], "..."
print "Assigned center: ", center_i.value[:4], "..."
print "Assigned precision: ", tau_i.value[:4], "..."

PyMC有一個MCMC的類，MCMC在PyMC的主命名空間中个曙，實現了MCMC的搜索算法锈嫩。我們通過下面的方法初始化一個Model實例：

mcmc=pm.MCMC( model )

讓MCMC搜索空間的方式是采樣（迭代），其中迭代就是你期望算法執(zhí)行的步數垦搬。我們嘗試50000步：

mcmc=pm.MCMC(model)
mcmc.sample(50000)
_____________
[****************100%******************]  50000 of 50000 complete

下面我們畫出來未知參數（中心呼寸，準確度和p）的path(或者軌跡)。我們這里可以使用MCMC對象中trace方法猴贰，這個接受的參數是PyMC變量名稱对雪。例如，mcmc.trace("centers")會檢索到一個可以被索引的Trace對象（使用[:]或者.gettrace()來查詢所有的路徑或者fancy-indexing[1000:]）

figsize(12.5, 9)
plt.subplot(311)
lw = 1
center_trace = mcmc.trace("centers")[:]

# for pretty colors later in the book.
colors = ["#348ABD", "#A60628"] \
if center_trace[-1, 0] > center_trace[-1, 1] \
    else ["#A60628", "#348ABD"]

plt.plot(center_trace[:, 0], label="trace of center 0", c=colors[0], lw=lw)
plt.plot(center_trace[:, 1], label="trace of center 1", c=colors[1], lw=lw)
plt.title("Traces of unknown parameters")
leg = plt.legend(loc="upper right")
leg.get_frame().set_alpha(0.7)

plt.subplot(312)
std_trace = mcmc.trace("stds")[:]
plt.plot(std_trace[:, 0], label="trace of standard deviation of cluster 0",
     c=colors[0], lw=lw)
plt.plot(std_trace[:, 1], label="trace of standard deviation of cluster 1",
     c=colors[1], lw=lw)
plt.legend(loc="upper left")

plt.subplot(313)
p_trace = mcmc.trace("p")[:]
plt.plot(p_trace, label="$p$: frequency of assignment to cluster 0",
     color="#467821", lw=lw)
plt.xlabel("Steps")
plt.ylim(0, 1)
plt.legend()

MCMC4.png

注意下面的特征：

1. 路徑收斂到米绕，不是到一個單一的點瑟捣，但是是到一個可能點的分布。這就是MCMC算法中的收斂性栅干。
2. 使用前面幾千個點作為推斷的依據不是很靠譜迈套，因為他們與最終的分布沒有太大的關聯。所以最好的方法是將前面的樣本去除掉碱鳞。我們成這個階段為燃燒階段（burn-in period）
3. 軌跡看起來像是在空間上進行隨機游走桑李，路徑表現出來的關聯性只是和前一個位置相關。這個有好處也有壞處：我們總是可以讓當前位置僅和下一位置相關；而無法很好的探索整個空間贵白。在本章最后我們會詳細討論率拒。

為了達到最終的收斂，我們將這行更多的MCMC步數禁荒。在MCMC已經被調用后重新開始并不是指重新整個算法猬膨。在上面給出的偽代碼中，只有當前的位置才是起作用的呛伴，這個被存放在PyMC變量value中勃痴。因此停止一個MCMC算法和觀察它的整個過程，然后重新開始是OK的磷蜀。value屬性不會被重寫召耘。

我們將會進行100000次取樣并且可視化這個過程：

mcmc.sample(100000)
[****************100%******************]  100000 of 100000 complete

figsize(12.5, 4)
center_trace = mcmc.trace("centers", chain=1)[:]
prev_center_trace = mcmc.trace("centers", chain=0)[:]

x = np.arange(50000)
plt.plot(x, prev_center_trace[:, 0], label="previous trace of center 0",
     lw=lw, alpha=0.4, c=colors[1])
plt.plot(x, prev_center_trace[:, 1], label="previous trace of center 1",
     lw=lw, alpha=0.4, c=colors[0])

x = np.arange(50000, 150000)
plt.plot(x, center_trace[:, 0], label="new trace of center 0", lw=lw, c="#348ABD")
plt.plot(x, center_trace[:, 1], label="new trace of center 1", lw=lw, c="#A60628")

plt.title("Traces of unknown center parameters")
leg = plt.legend(loc="upper right")
leg.get_frame().set_alpha(0.8)
plt.xlabel("Steps")

MCMC5.png

聚類分析

我們還未忘記自己的使命——識別其中的聚類衫贬。我們已經確定了未知量的后驗分布。我們在下面畫出了中心和標準方差的變量：

figsize(11.0, 4)
std_trace = mcmc.trace("stds")[:]

_i = [1, 2, 3, 0]
for i in range(2):
    plt.subplot(2, 2, _i[2 * i])
    plt.title("Posterior of center of cluster %d" % i)
    plt.hist(center_trace[:, i], color=colors[i], bins=30,
             histtype="stepfilled")

    plt.subplot(2, 2, _i[2 * i + 1])
    plt.title("Posterior of standard deviation of cluster %d" % i)
    plt.hist(std_trace[:, i], color=colors[i], bins=30,
             histtype="stepfilled")
    # plt.autoscale(tight=True)

plt.tight_layout()

MCMC6.png

import matplotlib as mpl
figsize(12.5, 4.5)
plt.cmap = mpl.colors.ListedColormap(colors)
plt.imshow(mcmc.trace("assignment")[::400, np.argsort(data)],
       cmap=plt.cmap, aspect=.4, alpha=.9)
plt.xticks(np.arange(0, data.shape[0], 40),
       ["%.2f" % s for s in np.sort(data)[::40]])
plt.ylabel("posterior sample")
plt.xlabel("value of $i$th data point")
plt.title("Posterior labels of data points")

MCMC7.png

看看上面的圖村砂，看起來最為不確定的是在150到170之間的烂斋。上面的圖可能會產生些許的不準確，因為xaxis不是真實的標量（它展示了第i個排序后的數據點的值）础废。更加清晰的圖我們在下面給出汛骂，其中我們已經估計了屬于label0和1的每個數據點的頻率。

cmap = mpl.colors.LinearSegmentedColormap.from_list("BMH", colors)
assign_trace = mcmc.trace("assignment")[:]
plt.scatter(data, 1 - assign_trace.mean(axis=0), cmap=cmap,
        c=assign_trace.mean(axis=0), s=50)
plt.ylim(-0.05, 1.05)
plt.xlim(35, 300)
plt.title("Probability of data point belonging to cluster 0")
plt.ylabel("probability")
plt.xlabel("value of data point")

MCMC9.png

即使我們建模cluster使用的是正態(tài)分布评腺，我們也沒有得到一個簡單的正態(tài)分布來很好的擬合數據帘瞭，但是獲得了一個正態(tài)分布的參數的值的分布。那么我們如何選擇一個均值和方差的對來確定一個擬合最后的正態(tài)分布呢蒿讥？

一個簡單粗暴的方式是使用后驗分布的均值图张。下面我們使用觀測的數據來更新了正態(tài)密度函數锋拖，使用后驗分布的均值作為被選擇的參數：

norm = stats.norm
x = np.linspace(20, 300, 500)
posterior_center_means = center_trace.mean(axis=0)
posterior_std_means = std_trace.mean(axis=0)
posterior_p_mean = mcmc.trace("p")[:].mean()

plt.hist(data, bins=20, histtype="step", normed=True, color="k",
     lw=2, label="histogram of data")
y = posterior_p_mean * norm.pdf(x, loc=posterior_center_means[0],
                                scale=posterior_std_means[0])
plt.plot(x, y, label="Cluster 0 (using posterior-mean parameters)", lw=3)
plt.fill_between(x, y, color=colors[1], alpha=0.3)

y = (1 - posterior_p_mean) * norm.pdf(x, loc=posterior_center_means[1],
                                      scale=posterior_std_means[1])
plt.plot(x, y, label="Cluster 1 (using posterior-mean parameters)", lw=3)
plt.fill_between(x, y, color=colors[0], alpha=0.3)

plt.legend(loc="upper left")
plt.title("Visualizing Clusters using posterior-mean parameters")

MCMC10.png

不要混淆后驗樣本

在上面的例子中，可能場景就是cluster 0有一個非常大的標準方差祸轮，而cluster1有個小一些的標準方差。這個可與你滿足這個證據侥钳，盡管會比我們初始時的推斷略少适袜。另外，它會變的不可能對兩個分布都有一個小的標準方差舷夺，因為數據不可能支持這樣的假設苦酱。因此兩個標準方差時依賴于彼此的：如果一個小，另外一個就會表達给猾。實際上疫萤，所有這些未知量都是以一個相似的方式關聯在一起的。例如敢伸，如果一個標準方差很大扯饶，均值就會在一個寬泛的空間中取值。相反地池颈，小的標準方差限制了均值落在一嬌小的區(qū)間內尾序。
在MCMC的過程中，我們返回表示從后驗分布中產生的樣本的向量躯砰。不同向量的元素不可以在一起使用每币，因為這個會打破上述邏輯：可能一個樣本已經返回了cluster 1 有一個小的標準方差，因此所有其他的變量在這個樣本中將會被影響并且被做出相應的調整琢歇。當然避免這樣的問題倒是很簡單兰怠，只需要你確保正確地索引了軌跡。

另一個小的例子可以解釋這點李茫。假設有兩個變量x和y揭保，其關聯是x+y=10。我們將x建模為一個正態(tài)隨機變量其均值是4涌矢，探索500個樣本掖举。

import pymc as pm

x = pm.Normal("x", 4, 10)
y = pm.Lambda("y", lambda x=x: 10 - x, trace=True)

ex_mcmc = pm.MCMC(pm.Model([x, y]))
ex_mcmc.sample(500)

plt.plot(ex_mcmc.trace("x")[:])
plt.plot(ex_mcmc.trace("y")[:])
plt.title("Displaying (extreme) case of dependence between unknowns")

MCMC11.png

回到聚類分析：預測

上面的聚類可以推廣至 k 個簇。選擇 k = 2 我們可以更好地對 MCMC 進行可視化名秀，并看看某些非常有趣的畫圖励负。
那關于預測呢？假設我們遇到一個新的數據點匕得，比方說 x = 175继榆。我們希望給它分配一個 簇巾表。直接將其分配到最近的簇中心是很傻的，因為這忽略了簇本身的標準差略吨。更加形式化的說法是集币，我們對分配cluster 1 給 x = 175的概率感興趣。記 x 的分配為 L_x（取值為 0 或者 1）翠忠，我們關注的是概率 P(L_x =1 | x = 175)鞠苟。
一種原始的方法就是使用額外的數據點來重新運行上面的 MCMC 。這個方法的缺點就是對每個新來的數據點的推斷會相當緩慢秽之。所以当娱，我們嘗試一個略微不精確，但是更快捷的方法考榨。
這里我們使用貝葉斯定理跨细。回想一下貝葉斯定理是這樣的：

貝葉斯定理

分配概率 1

因為分母都是一樣的溃槐，所以可以忽略分母（總算解脫了，因為計算這個分母的概率還是非常困難的）

分配概率 2

norm_pdf = stats.norm.pdf
p_trace = mcmc.trace("p")[:]
x = 175

v = p_trace * norm_pdf(x, loc=center_trace[:, 0], scale=std_trace[:, 0]>(1 - p_trace) * norm_pdf(x, loc=center_trace[:,1], scale=std_trace[:,1])

print "Probability of belonging to cluster 1: ", v.mean()
# Probability of belonging to cluster 1: 0.025

給出概率而非特定標簽其實是很有用的科吭。不再需要 L = 1 if prob > 0.5 else 0昏滴，我們可以使用損失函數來優(yōu)化我們的學長，這也是第五章我們會介紹的內容对人。

使用 MAP 來提升收斂效果

如果你自己運行上面的例子谣殊，可能會注意到我們的結果并不是一致的：可能你獲得的劃分更加分散，也可能是比較集中牺弄。問題在于我們的 trace 是 MCMC 算法初始值的函數姻几。
可以用數學的方式證明，讓 MCMC 運行足夠長時間势告，算法會忘記其初始位置蛇捌。實際上，這就是 MCMC 收斂的真正含義（在實踐中咱台，我們可能難以到達完全的收斂）络拌。因此如果我們觀察不同的后驗分析結果，可能就是因為我們的 MCMC 還沒有完全收斂回溺，因此我們還不能從這種狀態(tài)下的分布中進行采樣（也就是說春贸，需要更長的 burn-in 階段）混萝。
實際上，很差的初始值可以影響收斂的效果萍恕，或者極大地減緩收斂過程逸嘀。理想地看，我們希望讓 chain 從一個地貌的峰頂開始運行允粤，因為這正是后驗分布所在之處厘熟。因此，如果我們從峰頂開始维哈，就可以避免長期 burn-in 過程和錯誤的推斷慈省。一般地论巍，我們就叫這樣的峰頂為 MAP(maximum a posterior)。
當然炕柔，我們不知道 MAP 在那里脑蠕。 PyMC 提供了一個對象用來近似（如果不是確切找到） MAP 位置购撼。在 PyMC 的主命名空間中， MAP 對象接受一個 PyMC 模型實例谴仙。調用 MAP 實例中的 .fit() 方法可以設置模型中變量為他們對應的 MAP 值迂求。

map_ = pm.MAP( model )
map_.fit()

MAP.fit() 方法給擁有允許用戶選擇使用哪一種優(yōu)化算法的靈活性（總之，這是一個優(yōu)化問題：我們在尋找那些最大化地貌的值）晃跺，因為不是所有優(yōu)化算法都是適用的揩局。fit 適用的默認優(yōu)化算法就是 scipy 的 fmin 算法（嘗試最小化地貌的負值）。另一個可以使用算法是 Powell Method掀虎，PyMC 的 blogger Abraham FlaxmanAbraham Flaxman最喜歡的一種方法凌盯，可以使用 fit(method='fmin_powell) 得到。從我個人的經驗看來烹玉，默認的用得更多驰怎，但是如果收斂很慢的話，我也會嘗試 Powell method二打。
MAP 也可以看成是推斷問題的解县忌，因為數學上它就是未知量的最可能的值。但是本章前面也提到過继效，這種結果忽略了不確定性症杏，并沒有返回一個分布。
一般來說在調用 mcmc 之前使用 MAP(model).fit() 都是一個很好的想法莲趣。對 fit 的調用計算量不會非常重鸳慈，也能夠在一定能夠程度上節(jié)約花在 burn-in 階段的時間。

Burn-in 階段的信息

即使我們使用 MAP 先驗來調用 MCMC.sample喧伞，提供 Burn-in 階段還是很好的做法走芋，這樣可以更加安全绩郎。我們可以使用 PyMC 設置 burn 的參數來自動地丟棄前 n 個樣本。因為我們不知道 chain 什么時候已經完全收斂翁逞，我傾向于丟棄前面一半的樣本肋杖，甚至有時達到 90% 的樣本。接著前面的聚類例子挖函，我現在的代碼就像這樣：

model = pm.Model( [p, assignment, taus, centers ] )

map_ = pm.MAP( model )
map_.fit() #stores the fitted variables' values in foo.value

mcmc = pm.MCMC( model )
mcmc.sample( 100000, 50000 )

Diagnosing Convergence

Autocorrelation

自相關是一系列數字其本身相關性的度量状植。度量值為 1.0 則是完美的正自相關，0.0 則是沒有自相關怨喘，而 - 1.0 表示完美的負自相關津畸。如果你熟悉標準相關性定義，那么自相關就是一個序列 x_{\tau} 在時間 t 和 在時間 t-k 的相關性：

Paste_Image.png

例如必怜，下面兩個序列：

Paste_Image.png

擁有對應的樣本路徑：

figsize(12.5, 4)
import pymc as pm
x_t = pm.rnormal(0, 1, 200)
x_t[0] = 0
y_t = np.zeros(200)
for i in range(1, 200): 
    y_t[i] = pm.rnormal(y_t[i - 1], 1)
plt.plot(y_t, label="$y_t$", lw=3)
plt.plot(x_t, label="$x_t$", lw=3)
plt.xlabel("time, $t$")
plt.legend()

Paste_Image.png

一種考慮自相關性的角度是“如果我知道序列在時間 s 處的位置肉拓，這個信息能不能幫助我知道在時間 t 所處的位置？”在序列 x_t 處梳庆，答案是 不能暖途。通過構造，x_t 是隨機變量膏执。如果我告訴你 x_2 = 0.5驻售，你能夠告訴我更加準確的關于 x_3 的推測么 ？也不行更米。
換句話說欺栗，y_t 是自相關的。通過構造壳快，如果我知道 y_2 = 10纸巷，我可以很確定地認為 y_3
不會太偏離 10。類似的眶痰，我甚至可以對 y_4 進行一個相對不太肯定的推測：它不太可能靠近 0 或者 20瘤旨，但是 5 卻是很有可能的。我可以對 y_5 進行推測竖伯，但是這種置信度是更低的存哲。更加嚴格的表述是，我們承認在時間點之間的距離 k 增加時七婴，自相關性會隨之下降祟偷。可視化如下：

def autocorr(x):
    # from http://tinyurl.com/afz57c4
    result = np.correlate(x, x, mode='full')
    result = result / np.max(result)
    return result[result.size / 2:]

colors = ["#348ABD", "#A60628", "#7A68A6"]

x = np.arange(1, 200)
plt.bar(x, autocorr(y_t)[1:], width=1, label="$y_t$",
        edgecolor=colors[0], color=colors[0])
plt.bar(x, autocorr(x_t)[1:], width=1, label="$x_t$",
        color=colors[1], edgecolor=colors[1])

plt.legend(title="Autocorrelation")
plt.ylabel("measured correlation \nbetween $y_t$ and $y_{t-k}$.")
plt.xlabel("k (lag)")
plt.title("Autocorrelation plot of $y_t$ and $x_t$ for differing $k$ lags.")

Paste_Image.png

注意打厘，在隨著 k 增加時修肠，y_t 的自相關性會從一個很高的點開始下降。對比 x_t 的自相關性（非常像噪聲）户盯，我們可以說自相關性存在于這些序列中嵌施。

How does this relate to MCMC convergence

根據 MCMC 算法的本質饲化，我們總是返回那些顯示出自相關性的樣本（這是因為從你當前位置的步驟，移動到靠近你的另外一個）
能夠很好滴遍歷空間 chain 會表現出很高的自相關性吗伤。從可視化圖上看吃靠，如果 trace 看起來如河流一樣蜿蜒曲折，并且沒有靜止下來足淆，這個 chain 就會有很高的自相關性巢块。
這并不是說，一個收斂的 MCMC 擁有低自相關性巧号。因此低自相關性并不是收斂的必要條件族奢。PyMC 在模塊 Matplot 中擁有內置的畫出自相關性函數。

Thinning

如果在后驗分布的樣本中存在高自相關性丹鸿，會出現另一個問題歹鱼。很多后期處理算法需要樣本之間是獨立的。通過僅僅返回第 n 個樣本來解決或者至少緩解這個問題卜高，可以去除一些自相關性。下面我們執(zhí)行對 y_t 的自相關性畫圖南片，按照不同 thinning 的層次進行：

max_x = 200 / 3 + 1
x = np.arange(1, max_x)

plt.bar(x, autocorr(y_t)[1:max_x], edgecolor=colors[0],
        label="no thinning", color=colors[0], width=1)
plt.bar(x, autocorr(y_t[::2])[1:max_x], edgecolor=colors[1],
        label="keeping every 2nd sample", color=colors[1], width=1)
plt.bar(x, autocorr(y_t[::3])[1:max_x], width=1, edgecolor=colors[2],
        label="keeping every 3rd sample", color=colors[2])

plt.autoscale(tight=True)
plt.legend(title="Autocorrelation plot for $y_t$", loc="lower left")
plt.ylabel("measured correlation \nbetween $y_t$ and $y_{t-k}$.")
plt.xlabel("k (lag)")
plt.title("Autocorrelation of $y_t$ (no thinning vs. thinning) \
at differing $k$ lags.")

Paste_Image.png

pymc.Matplot.plot()

在每次執(zhí)行 MCMC 時手動創(chuàng)建直方圖可霎、自相關畫圖和 trace 畫圖并不是很方便。PyMC 的作者已經包含了這樣的可視化工具宴杀。
正如標題所說的癣朗，pymc.Matplot 模塊包含一個命名不太好的函數 plot，我傾向于作為 mcplot 導入旺罢，這樣可以避免一些與其他的命名空間的 plot 或者 mcplot 的沖突旷余。plot 接受 MCMC 對象返回每個變量（最多 10 個變量）的后驗分布绢记，trace 和 自相關性。
下面我們使用這個工具來畫出設置 thinning = 10在25 000次采樣后簇的中心荣暮。

from pymc.Matplot import plot as mcplot

mcmc.sample(25000, 0, 10)
mcplot(mcmc.trace("centers", 2), common_scale=False)

# [****************100%******************]  25000 of 25000 completePlotting centers_0
# Plotting centers_1

Paste_Image.png

名稱為 center_0_acorr 和 center_1_acorr 是生成的自相關性的畫圖。這看起來和之前展示的那些不同诈泼，但是僅僅在于 0-lag 在不在圖的中央懂拾。

Useful tips for MCMC

如果不是MCMC的計算上造成的困難，那么貝葉斯推斷的首選肯定是這個方法铐达。實際上岖赋，MCMC 是讓很多的人都放棄了在實際問題使用它。下面給出一些好的啟發(fā)式建議來幫助快速收斂瓮孙。

Intelligent starting values

從后驗分布附近開始 MCMC 算法肯定很有效贾节，這樣只需要少量的時間就可以正確地進行采樣了。我們可以通過在創(chuàng)建 Stochastic 變量時指定 value 參數來幫助算法衷畦。在很多情況下栗涂，我們可以產生關于參數的合理的猜測。例如祈争，如果我們有來自正態(tài)分布的數據斤程，想要估計 \mu 參數，那么就可以從數據的均值來作為開始點。

 mu = pm.Uniform( "mu", 0, 100, value = data.mean() )

對大多數的模型參數忿墅，有一個頻率主義的估計扁藕。這些估計是 MCMC 算法很好的起始點。當然疚脐，這并不總是可行的亿柑，但是包含盡可能多的合適的初始值總歸是很好的想法。即使你猜測的值是錯誤的棍弄，MCMC 仍然會收斂到合適的分布望薄，很少會失敗。
這就是 MAP 試著做出的事情呼畸，通過給 MCMC 設置出好的初始值痕支。所以為何還需要指定用戶自定義的值呢？因為蛮原，給定好的 MAP 值就能夠幫助它找到 MAP卧须。
同樣重要的一點是，壞的初始值也是 PyMC 主要的問題之一儒陨，會導致收斂失敗花嘶。

Priors

如果先驗選擇得不好，MCMC 算法可能不會收斂蹦漠，至少收斂困難察绷。想想如果選擇的先驗沒有能夠包含真實的參數：先驗賦給未知量概率為 0，這樣后驗同樣得到 0 的概率津辩。這就出問題了。
所以容劳，最好精心選擇先驗喘沿。同樣，收斂的缺失或者樣本集中在邊界上的情況就是因為選錯了先驗（參看下面的 統(tǒng)計計算的 Folk Theorem）竭贩。

Covariance matrices and eliminating parameters

The Folk Theorem of Statistical Computing

如果你在計算上出了問題蚜印，可能你的模型就已經錯了。

Conclusion

PyMC 給出了執(zhí)行貝葉斯推斷的強大支持平臺留量，因為其擁有抽象的 MCMC 的內部機制窄赋。除此之外，在確保推斷本身沒有受到 MCMC 的迭代特性的影響也是需要注意的楼熄。

References

Flaxman, Abraham. "Powell's Methods for Maximization in PyMC." Healthy Algorithms. N.p., 9 02 2012. Web. 28 Feb 2013. http://healthyalgorithms.com/2012/02/09/powells-method-for-maximization-in-pymc/.