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

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：统计计算的无名定理

如果你的计算遇到了问题，你的模型有可能是错误的。即先验有问题，或者选择的分布函数有问题。