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

時間 2019-12-01

標籤打開 mcmc 蒙特卡洛黑盒子 pymc 貝葉推理底層實現原理初探简体版

原文原文鏈接

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

1. 貝葉斯景象圖

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

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

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

咱們這裏選擇2維的，即包含2個隨機變量的貝葉斯推斷問題，進行可視化展現，選擇2維是由於能夠方便進行可視化，高維空間是很難想象的。dom

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)點開始的瀑布。post

# -*- 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過程當中，觀測樣本只是經過拉伸和壓縮先驗分佈的曲面，讓曲面更符合實際的參數分佈，以代表參數的真實值最可能在哪裏。this

數據𝑋越多拉伸和壓縮就越厲害，這樣後驗分佈就變化的越厲害，可能徹底看不出和先驗分佈的曲面有什麼關係，或者說隨着數據的增長先驗分佈對於後驗分佈的影響愈來愈小。這也體現了貝葉斯推斷的核心思想：你的先驗應該儘量合理，可是即便不是那麼的合理也沒有太大關係，MCMC會經過觀測樣本的擬合，儘量將先驗分佈調整爲符合觀測樣本的後驗分佈。google

可是若是數據𝑋較小，那麼後驗分佈的形狀會更多的反映出先驗分佈的形狀。在小樣本的狀況下，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()