我們在這篇文章里有嘗試討論三個重點。第一，討論的 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：統計計算的無名定理
如果你的計算遇到了問題，你的模型有可能是錯誤的。即先驗有問題，或者選擇的分布函數有問題。