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

从图中能够看出，随着𝑘的增长，𝑦𝑡 的自相关性从最高点快随减少。𝑥𝑡 的自相关性看起来就像噪声同样(实际上它就是噪声，白噪声)，所以能够得出𝑥𝑡 没有自相关性。

笔者思考：读者还记得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：统计计算的无名定理

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