開啟MCMC（馬爾科夫蒙特卡洛）的黑盒子 - Pymc貝葉斯推理底層實現原理初探

我們在這篇文章裡有嘗試討論三個重點。第一，討論的 MCMC。第二，學習 MCMC 的實現過程，學習 MCMC 演算法如何收斂，收斂到何處。第三，將會介紹為什麼從後驗分佈中能返回成千上萬的樣本，也許讀者和我一樣，剛開始學習時，面對這種取樣過程看起來有點奇怪。

1. 貝葉斯景象圖

當構造一個有?個未知變數的貝葉斯推斷問題時，首先要隱式的建立 N 維空間（可以理解為 N 個隨機變數）的先驗分佈。

這 N 個隨機變數的分佈，是一個曲面或者曲線， 曲線或者曲面反映了一個點的概率，即概率密度函式。我們隱式地建立了一個 N 維空間。

0x1：先驗分佈的視覺化圖景像

我們這裡選擇2維的，即包含2個隨機變數的貝葉斯推斷問題，進行視覺化展示，選擇2維是因為可以方便進行視覺化，高維空間是很難想象的。

1. 二維均勻分佈的先驗圖景像

例如有兩個未知的隨機變數 ?1 和 ?2，這兩個隨機變數的先驗分佈為 Uniform(0,5)，?1和?2構成了一個5 × 5的平面空間， 代表概率密度的曲面位於這個5 × 5的平面空間上（由於假定了均勻分佈，因此每一點的概率都相同）

# -*- coding: utf-8 -*-

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_y[:, None], uni_x[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")

plt.show()

2. 二維指數分佈的先驗圖景像

換一個例子，如果兩個先驗分佈是 Exp(3)和 Exp(10)，則構成的空間是在二維平面上的正數，同時表示先驗概率的曲面就像從(0,0)點開始的瀑布。

# -*- coding: utf-8 -*-

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 為[0,5]內的均勻分佈
X, Y = np.meshgrid(x, y)

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_y[:, None], exp_x[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")

plt.show()

圖中顏色越是趨向於暗紅的位置，其先驗概率越高。反過來，顏色越是趨向於深藍的位置，其先驗概率越低。

0x2：觀測樣本是如何影響未知變數的先驗分佈的？

概率面描述了未知變數的先驗分佈，而觀測樣本的作用我們可以形象地理解為一隻手，每次來一個樣本，這隻手就根據觀測樣本的情況，將先驗分佈的曲面向“符合”觀測樣本的方向拉伸一點。

在MCMC過程中，觀測樣本只是通過拉伸和壓縮先驗分佈的曲面，讓曲面更符合實際的引數分佈，以表明引數的真實值最可能在哪裡。

資料?越多拉伸和壓縮就越厲害，這樣後驗分佈就變化的越厲害，可能完全看不出和先驗分佈的曲面有什麼關係，或者說隨著資料的增加先驗分佈對於後驗分佈的影響越來越小。這也體現了貝葉斯推斷的核心思想：你的先驗應該儘可能合理，但是即使不是那麼的合理也沒有太大關係，MCMC會通過觀測樣本的擬合，儘可能將先驗分佈調整為符合觀測樣本的後驗分佈。

但是如果資料?較小，那麼後驗分佈的形狀會更多的反映出先驗分佈的形狀。在小樣本的情況下，MCMC對初始的先驗分佈會非常敏感。

1. 不同的先驗概率對觀測樣本調整後驗分佈的阻力是不同的

需要再次說明的是，在高維空間上，拉伸和擠壓的變化難以視覺化。在二維空間上，這些拉伸、擠壓的結果是形成了幾座山峰。而形成這些區域性山峰的作用力會受到先驗分佈的阻撓。

先驗分佈越小意味著阻力越大；先驗分佈越大阻力越小。

有一點要特別注意，如果某處的先驗分佈為0，那麼在這一點上也推不出後驗概率。

我們通過兩個引數為 λ 的泊松分佈進行估計，它們分別用均勻分佈作為先驗，以及指數分佈作為先驗，我們來觀察在獲得一個觀察值前後的不同景象。

# -*- 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
import scipy.stats as stats

jet = plt.cm.jet

# 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, :])


# 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_y[:, None], uni_x[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_y[:, None], exp_x[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)

plt.show()

四張圖裡的黑點代表引數的真實取值。每列的上圖表示先驗分佈圖形，下圖表示後驗分佈圖形。

我們主要到，雖然觀測值相同，兩種先驗假設下得到的後驗分佈卻有不同的圖形。

我們從上圖裡注意到2個細節：

1. 右下方的指數先驗對應的後驗分佈圖形中，右上角區域的取值很低，原因是假設的指數先驗在這一區域的取值也較低。
2. 另一方面，左下方的均勻分佈對應的後驗分佈圖形中，右上角區域的取值相對較高，這也是因為均勻先驗在該區域的取值相比指數先驗更高。
3. 在右下角指數分佈的後驗圖形中，最高的山峰，也就是紅色最深的地方，向（0，0）點偏斜，原因就是指數先驗在這個角落的取值更高。

筆者的思考：這個現象其實和深度學習裡sigmoid函式的調整過程是類似的，sigmoid在越靠近0或1概率的區域中，調整的速率會越來越慢，即死衚衕效應。因為這時候sigmoid會認為收斂已經接近尾聲，要減緩調整的幅度，以穩定已經收斂到的結果。

0x3：使用MCMC來搜尋圖景像

要想找到後驗分佈上的那些山峰，我們需要對整個後驗空間進行探索，這一個空間是由先驗分佈的概率面以及觀測值一起形成的。但是，很遺憾的是，理論和實踐之間常常存在一定的鴻溝，遍歷一個N維空間的複雜度將隨著N呈指數增長，即隨著N增加，N維空間的大小將迅速爆發。毫無疑問，我們不能靠蠻力去窮舉搜尋，不過話說回來，如果未來量子計算技術面世，也許計算複雜度的問題將不復存在，那麼當今的很多所謂的啟發式搜素演算法都將退出歷史舞臺。

其實，MCMC背後的思想就是如何聰明地對空間進行搜尋，搜尋什麼呢？一個點嗎？肯定不是，貝葉斯思想的核心就是世界上沒有100%確定的東西，所有的推斷都是一個概率分佈。

實際上，MCMC的返回值是後驗分佈上的“一些樣本”，而非後驗分佈本身。我們可以把MCMC的過程想象成不斷重複地問一塊石頭：“你是不是來自我要找的那座山？”，並試圖用上千個肯定的答案來重塑那座山。最後將它們返回並大功告成。

在MCMC和PyMC的術語裡，這些返回序列裡的“石頭”就是觀測樣本，累計起來稱之為“跡”。

筆者思考：“跡”實際上就是模型當前啟發式探索到的向量空間位置的軌跡，類似於SGD優化過程中的路徑。

我們希望MCMC搜尋的位置能收斂到後驗概率最高的區域（注意不是一個確定的點，是一個區域）。為此，MCMC每次都會探索附近位置上的概率值，並朝著概率值增加的方向前進。MCMC朝著空間裡的一大塊區域移動，並繞著它隨機遊走，順便從區域中採集樣本。

1. 為何是上千的樣本？

我們可能會說，演算法模型訓練的目的不就是為了獲得我們對隨機變數的最優估計嗎？畢竟在很多時候，我們訓練得到的模型會用於之後的預測任務。但是貝葉斯學派不這麼做，貝葉斯推斷的結果更像是一個參謀，它只提供一個建議，而最後的決策需要我們自己來完成（例如我們通過取後驗估計的均值作為最大後驗估計的結果）

回到MCMC的訓練返回結果，它返回成千上萬的樣本讓人覺得這是一種低效的描述後驗概率的方式。實際上這是一種非常高效的方法。下面是其它可能用於計算後驗概率的方法

1. 用解析表示式描述“山峰區域”(後驗分佈)，這需要描述帶有山峰和山谷的 N 維曲面。在數學上是比較困難的。
2. 也可以返回圖形裡的頂峰，這種方法在數學上是可行的，也很好理解（這對應了關於未知量的估計裡最可能的取值），但是這種做法忽略了圖形的形狀，而我們知道，在一些場景下，這些形狀對於判定未知變數的後驗概率來說，是非常關鍵的
3. 除了計算原因，另一個主要原因是，利用返回的大量樣本可以利用大數定理解決非常棘手問題。有了成千上萬的樣本，就可以利用直方圖技術，重構後驗分佈曲面。

2. MCMC的演算法實現

有很多演算法實現 MCMC。大部分的演算法可以表述如下

1. 選定初始位置。即選定先驗分佈的初始位置。
2. 移動到一個新的位置(探索附件的點)。
3. 根據新資料位置和先驗分佈的緊密程度拒絕或接受新的資料點(石頭是否來自於期望的那座山上)。
4. 
    A. 如果接受，移動到新的點，返回到步驟 1。
    B. 否則不移動到新的點返回到步驟 1。
5. 在經過大量的迭代之後，返回所有接受的點。

這樣，我們就從整體上向著後驗分佈所在的方向前進，並沿途謹慎地收集樣本。而一旦我們達到了後驗分佈所在的區域，我們就可以輕鬆地採集更多的樣本，因為那裡的點幾乎都位於後驗分佈的區域裡。

我們注意到，在MCMC的演算法步驟中，每一次移動僅與當前位置有關（即總是在當前位置的附近進行嘗試）。我們將這一特性稱之為“無記憶性”，即演算法並不在乎如何達到當前位置，只要達到即可。

0x4：一個MCMC搜尋圖景像的實際的例子 - 使用混合模型進行無監督聚類

假定使用以下資料集：

# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
import pymc as pm
import scipy.stats as stats

jet = plt.cm.jet

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], "...")

fig = plt.figure()
ax = fig.add_subplot(111)
ax.hist(data[:], bins=20)

plt.show()

從直觀上觀察，這些資料說明了什麼？似乎可以看出資料具有雙模的形式，也就是說，圖中看起來有兩個峰值，一個在120附近，另一個在200附近。那麼，資料裡可能存在兩個聚類簇。

筆者思考：聚類是一個很寬泛的概念，不一定僅限於我們所熟悉的歐式空間的kmeans，實際上，聚類不一定都是幾何意義上的聚類。通過對原始資料集進行概率分佈的擬合，從而獲得資料集中每個點所屬類別的概率分佈，這也是一種聚類。

1. 選擇符合資料觀測分佈的數學模型

首先，根據數學建模的思想，我們選擇正態分佈作為擬合資料的數學模型。

下面介紹資料是如何產生的。生成的資料方法如下

1. 對於每個資料點，讓它以概率?屬於類別 1，以概率1-?屬於類別 2。
2. 畫出均值 ?? 和方差 ?? 的正態分佈的隨機變數，?是 1 中的類別 id，可以取 1 或者 2。
3. 重複 1，2 步驟。 

這個演算法可以產生與觀測資料相似的效果。所以選擇這個演算法作為模型。

但是現在的問題是我們不知道引數 ? 和正態分佈的引數。所以要學習或者推斷出這些未知變數。

用Nor0，Nor1分別表示正態分佈。兩個正態分佈的引數都是未知的，引數分別表示為??,??，? = 0，1。

2. 對模型的引數進行先驗建模

1）所屬類別分佈先驗

對於某一個具體的資料點來說，它可能來自Nor0也可能來自Nor1， 假設資料點來自Nor0的概率為?。 這是一個先驗，由於我們並不知道來自 Nor1 的實際概率，因此我們只能用 0-1 上的均勻分佈來進行建模假設（最大熵原理）。我們稱該先驗為 p。

有一種近似的方法，可以使用 PyMC 的類別(Categorical)隨機變數將資料點分配給某一類別。PyMC 類別隨機變數有一個?維概率陣列變數，必須對?維概率陣列變數進行 求和使其和變成 1，PyMC 類別隨機變數的 value 屬性是一個 0 到? − 1的值，該值如何選 擇由概率陣列中的元素決定(在本例中? = 2)。

目前還不知道將資料分配給類別 1 的 先驗概率是多少，所以選擇 0 到 1 的均勻隨機變數作為先驗分佈。此時輸入類別變數的 概率陣列為[?, 1 − ?]。

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], "..."

2）單個聚類中的正態分佈的引數分佈先驗

從對資料的觀察來看，我們會猜測兩個正態分佈具有不同的標準差。同樣，我們並不知道兩個標準差分別是多少，我們依然使用0-100上的均勻分佈對其進行建模。

此外，我們還需要兩個聚類簇中心點的先驗分佈，這兩個中心點的位置其實就是正態分佈的引數 ?。

雖然對肉眼的估計不是那麼有信心，但從資料形狀上看，我們還是猜測在?0 = 120，?1 = 190附近。不要忘了，MCMC會幫我們修正先驗中不是那麼精確的部分，所以不要擔心我們的先驗估計會存在誤差。

 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], "..."

3. MCMC搜尋過程 - 跡

PyMC 有個 MCMC 類，MCMC 位於 PyMC 名稱空間中，其實現了 MCMC 演算法。用 Model 物件例項化 MCMC 類，如下

mcmc = pm.MCMC( model )

用 MCMC 的 sample(iterations)方法可以用於探究隨機變數的空間，iteration 是演算法執行的步數。下面將演算法設定為執行 50000 步數

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

注意這裡的步數和樣本數不一定是相等的，事實上，步數一般是遠大於樣本數的，因為MCMC的剛開始的時候是在“試探性前進搜尋的”，最開始的搜尋一般都是抖動很劇烈的，這種抖動直到搜尋的後期會逐漸穩定下來。

下面畫出變數未知引數(中心點，精度和 ?)的路徑(“跡”)，要想得到跡，可以通過向 MCMC 物件中的 trace方法傳入想要獲取的PyMC變數，例如：

mcmc.trace("centers")
或者可以使用[:]或者.gettrance()得到所有的跡

程式碼如下：

# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
import pymc as pm
import scipy.stats as stats

jet = plt.cm.jet

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

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], "...")

stds = pm.Uniform("stds", 0, 100, size=2)
taus = 1.0 / stds ** 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], "...")

# and to combine it with the observations:
observations = pm.Normal("obs", center_i, tau_i, value=data, observed=True)

# below we create a model class
model = pm.Model([p, assignment, observations, taus, centers])


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

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()

plt.show()

從上圖我們可以看出什麼？

1. 這些跡並非收斂到某一點，而是收斂到一定分佈下，概率較大的點集。這就是MCMC演算法裡收斂的涵義。
2. 最初的幾千個點（訓練輪）與最終的目標分佈關係不大，所以使用這些點參與估計並不明智。一個聰明的做法是剔除這些點之後再執行估計，產生這些遺棄點的過程稱為預熱期。
3. 這些跡看起來像是在圍繞空間中某一區域隨機遊走。這就是說它總是在基於之前的位置移動。這樣的好處是確保了當前位置與之前位置之間存在直接、確定的聯絡。然而壞處就是太過於限制探索空間的效率。

4. 如何估計各個未知變數的最佳後驗估計值

上一小節討論到了跡和收斂，很顯然，我們會問一個問題：so what？

別忘了我們最大的挑戰仍然是識別各個聚類，而此刻我們的估計只能告訴我們各個未知變數的後驗概率分佈，我們把中心點和標準差的後驗分化畫出來。

# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
import pymc as pm
import scipy.stats as stats

jet = plt.cm.jet

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

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], "...")

stds = pm.Uniform("stds", 0, 100, size=2)
taus = 1.0 / stds ** 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], "...")

# and to combine it with the observations:
observations = pm.Normal("obs", center_i, tau_i, value=data, observed=True)

# below we create a model class
model = pm.Model([p, assignment, observations, taus, centers])

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

# 獲取整條跡
std_trace = mcmc.trace("stds")[:]
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"]

_i = [1, 2, 3, 4]
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()

plt.show()

可以看到，MCMC 演算法已經給出聚類中心最可能的地方分別在 120 和 200 處。同理，標準方差也可以有類似的推斷過程。
同時也給出同一類別下的資料的後驗分佈（資料點屬於哪一類）

# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
import pymc as pm
import scipy.stats as stats
import matplotlib as mpl

jet = plt.cm.jet

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

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], "...")

stds = pm.Uniform("stds", 0, 100, size=2)
taus = 1.0 / stds ** 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], "...")

# and to combine it with the observations:
observations = pm.Normal("obs", center_i, tau_i, value=data, observed=True)

# below we create a model class
model = pm.Model([p, assignment, observations, taus, centers])

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

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.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")

plt.show()

x 軸表示先驗資料的排序值。 紅色的地方表示類別 1，藍色地方代表類別 0。

在下圖中估計出了每個資料點屬於 0 和屬於 1 的頻。

# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
import pymc as pm
import scipy.stats as stats
import matplotlib as mpl

jet = plt.cm.jet

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

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], "...")

stds = pm.Uniform("stds", 0, 100, size=2)
taus = 1.0 / stds ** 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], "...")

# and to combine it with the observations:
observations = pm.Normal("obs", center_i, tau_i, value=data, observed=True)

# below we create a model class
model = pm.Model([p, assignment, observations, taus, centers])

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

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"]


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")

plt.show()

可以看到，雖然對正態分佈對兩類資料進行了建模，MCMC也根據觀測樣本得到了未知變數的後驗概率分佈。但是我們仍然沒有得到能夠最佳匹配資料的正態分佈，而僅僅是得到了關於正態分佈各引數的分佈。當然，這也體現了貝葉斯推斷的一個特點，貝葉斯推斷並不直接作出決策，它更多地是提供線索和證據，決策還是需要統計學家來完成。

那接下來一個很自然的問題是，我們如何能夠選擇能夠滿足最佳匹配的引數 - 均值、方差呢？

一個簡單粗暴的方法是選擇後驗分佈的均值（當然，這非常合理且有堅實的理論支撐）。在下圖中，我們以後驗分佈的均值作為正態分佈的各引數值，並將得到的正態分佈於觀測資料形狀疊加到一起。

# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
import pymc as pm
import scipy.stats as stats
import matplotlib as mpl

jet = plt.cm.jet

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

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], "...")

stds = pm.Uniform("stds", 0, 100, size=2)
taus = 1.0 / stds ** 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], "...")

# and to combine it with the observations:
observations = pm.Normal("obs", center_i, tau_i, value=data, observed=True)

# below we create a model class
model = pm.Model([p, assignment, observations, taus, centers])

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

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"]

std_trace = mcmc.trace("stds")[:]

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")

plt.show()

可以看到，取均值作為後驗比較好的“擬合”了觀測資料

5. 回到聚類：預測問題 - 到了該決策的時候了！

前面說了那麼多，假設我們使用均值作為最大後驗估計，好！我們現在得到了模型引數的估計了。接下來我們要來解決最後一步問題，即預測。

假設有個新的點? = 175，我們想知道這個點屬於哪個類別。

更一般的情況是：我們更感興趣? = 175這個點屬於類別 1 的概率是多大。將?屬於某一類別表示成??，我 們感興趣的是?(??|? = 175)。顯然，預測問題被轉化成了概率計算以及比大小問題。

下面將使用貝葉斯定理解決後驗推斷問題。貝葉斯定理如下

在此時的例子中?用?? = 1表示，?是觀測到的資料，即? = 175。對於後驗分佈的 引數(?0, σ0, ?1, σ1, ?)，我們感興趣的是“x 屬於 1 的概率是不是比 x 屬於 0 的概率要大?”，概率的大小和選擇的引數有關

因為分母都是一個樣的所以可以將其忽略掉：

# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
import pymc as pm
import scipy.stats as stats
import matplotlib as mpl

jet = plt.cm.jet

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

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], "...")

stds = pm.Uniform("stds", 0, 100, size=2)
taus = 1.0 / stds ** 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], "...")

# and to combine it with the observations:
observations = pm.Normal("obs", center_i, tau_i, value=data, observed=True)

# below we create a model class
model = pm.Model([p, assignment, observations, taus, centers])

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

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"]

std_trace = mcmc.trace("stds")[:]


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.06006)

 

2. MCMC收斂性討論

0x1：使用MAP來改進收斂性

如果我們重複多次執行上面的程式，我們會發現每次得到的結果都不盡相同。問題在於，MCMC得到的跡實際上是演算法起始值的函式，即MCMC是初始值敏感的。這也很自然，MCMC的搜尋過程是在做啟發式搜尋，類似於“盲人摸象”的過程，所以很自然地，不用的起點，其之後走的跡自然也是不同的。

從數學上可以證實，如果讓MCMC通過更多的取樣執行得足夠久，就可以忽略起始值的位置，其實這也就是MCMC收斂的定義所在。

因此，如果我們看到不同的後驗分佈結果，那可能是因為MCMC還沒有充分地收斂，此時的樣本還不適合用作分析（我們應該加大預熱期）。實際上，錯誤的起始位置可能阻礙任何的收斂，或者使之遲緩。

理想情況下，我們希望其實位置就分佈在概率圖形的山峰處，因為這其實就是後驗分佈的所在區域，如果以山峰為起點，就能避免很長的預熱期以及錯誤的的估計結果。通常，我們將山峰位置稱為最大後驗，或簡稱為MAP。

筆者思考：我們常常會在深度學習專案中，直接基於resNet、googleNet這種已經經過訓練優化後的模型，其背後的思想也有一些減少預熱期的意思，在resNet、googleNet的基礎上，在繼續進行訓練，可以更快地收斂到我們的目標概率分佈上。

當然，MAP的真實位置是未知的。PyMC提供了一個用於估計該位置的物件。MAP 物件位於 PyMC 的主名稱空間中，它接受一 個 PyMC 的 Model 例項作為引數。呼叫 MAP 例項的 fit()方法，將模型中的引數設定成 MAP 的值

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

理論上，MAP也能直接當做問題的解，因為從數學上說它是未知量最可能的取值。實際上，MAP最大後驗估計在工程專案中被使用到的頻率很高。

但是另一方面，我們要明白：這種單個點的解會忽略未執行，也無法得到分佈形式的返回結果。

常常在呼叫mcmc前呼叫方法 MAP(model).fit()。fit（）函式本身開銷並不大，但是卻能顯著減少預熱時間。

0x2：如何判斷收斂？

1. 自相關 - 序列遞迴推演性

要回答如何判斷收斂，我們需要先來討論一個概念：自相關！

自相關性用於衡量一串數字與自身的相關程度。1 表示和自身完全正相關，0 表示沒有自相關性，-1 表示完全負相關。如果你熟悉標準方法，自相關就是序列??在時刻 ?和 ? − ?的相關程度。

例如有如下兩個序列

??~Normal(0,1)，?0 =1
??~Normal(??−1,1)，?0 = 1

這兩個序列的影象如下

# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
import pymc as pm
import scipy.stats as stats
import matplotlib as mpl


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()

plt.show()

可以這樣來理解自相關，“如果知道在時刻?序列的位置，它能否幫助我們瞭解序列在時刻?時的位置在哪裡”。

在序列?? 中，是無法通過?時刻的序列位置判斷?時刻的序列位置的。如果知道?2 = 0.5，能否有利於猜測?3的位置在哪裡，答案是不能。

另一方面??是自相關的。如果知道?2 = 10，可以很自信的說?3離 10 的位置不會太遠。對?4可以進行類似的猜測(猜測的結果的可信度應該比?3要小些)，它不可能離 0 或者 20 比較近，離 5 近的可能性較大，對?5的猜測的可信度又比?4小。

所以可以得出如下結論，隨著?的減小，即資料點之間的間隔變小，相關性會增加。這個結論如下圖

# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
import pymc as pm
import scipy.stats as stats
import matplotlib as mpl


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)

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.")

plt.show()

從圖中可以看出，隨著?的增加，?? 的自相關性從最高點快隨減小。?? 的自相關性看起來就像噪聲一樣(實際上它就是噪聲，白噪聲)，因此可以得出?? 沒有自相關性。

筆者思考：讀者還記得HMM隱馬爾科夫假設嗎？即當前的節點狀態只和之前有限步驟（例如1步）的節點狀態有關，雖然理論上應該是和歷史上所有的節點狀態有相關，但是其實越往前，相關性越小，甚至小到可以忽略，因為HMM的假設實際上並沒有丟失太多的概率資訊。

2. 自相關性和MCM的收斂有何關係？

MCMC演算法會天然地返回具有自相關性的取樣結果，這是因為MCMC的“摸索性行走”演算法：總是從當前位置，移動到附近的某個位置。

從影象上看，如果取樣的跡看起來像蜿蜒緩慢、流淌不停地河流，那麼該過程的自相關性會很高。想象我們是河流裡的一粒水分子，如果我知道你現在在什麼位置，那麼我可以較為精確地估計你下一步會在哪。

而另一方面，如果過程的自相關性很低，那麼此時你很難估計你下一步的位置。

 

3. MCMC的一些小技巧

0x1：聰明的初始值

初始選擇在後驗概率附近，這樣花很少的時間就可以計算出正確結果。我們可以在建立隨機過程變數的時候，通過指定value引數來告訴演算法我們猜測後驗分佈會在哪裡。

在許多情況下我們可以為引數猜測一個合理的結果。例如，如果有資料符合正態分佈，希望估計引數?，最優的初值就是資料的均值：

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

對於大多數的模型引數，都可以以頻率論進行猜測，以得到良好的MCMC演算法處置。

當然了，有些引數是無法估計出頻率值的，不過還是要儘量的近似的初始值，這樣有利於計算。即使你的猜測是錯誤的，MCMC 也會收斂到合適的分佈，所以即使估計錯誤也沒有什麼損失。
在實際的工程專案中，我們應該儘可能結合我們的業務領域經驗，將領域經驗融合到先驗初始值中。

0x2：先驗

另外重要的一點是，不合適的初值是 PyMC 出 bug 的主要原因，同時也不利於收斂。

如果先驗分佈選擇的不好，MCMC 演算法可能不會收斂，至少收斂比較困難。考慮下：如果先驗分佈沒有包含真的引數，類別 0 的先驗概率將是未知的，因此屬於類別 0 的後 驗概率也將未知。這可能會導致病態的結果。
所以要小心選擇先驗分佈。

一般來說，如果收斂性不好或者樣本擁擠的分佈在某一個邊界，說明先驗分佈選擇的有問題。

0x3：統計計算的無名定理

如果你的計算遇到了問題，你的模型有可能是錯誤的。即先驗有問題，或者選擇的分佈函式有問題。