还在抛硬币找答案？还在为抛硬币不理想而烦恼？现在让Python来帮助你抛硬币吧，又准又公道。

PyMC3是一个用于概率编程的Python库，当前最新的版本号是2016年10月4号发布的3.0.rc2。PyMC3提供了一套非常简洁直观的语法，非常接近统计学中描述概率模型的语法，可读性很高。PyMC3是用Python写的，其中的核心计算部分基于NumPy和Theano。Theano是一个用于深度学习的Python库，可以高效地定义、优化和求解多维数组的数学表达式。PyMC3使用Theano的主要原因是某些采样算法（如NUTS）需要计算梯度，而Theano可以很方便地进行自动求导。而且，Theano将Python代码转化成了C代码，因而PyMC3的速度相当快。

用计算的方法解决抛硬币问题

让我们重新回顾下抛硬币问题，这次我们使用PyMC3。首先我们需要获取数据，这里我们使用手动构造的数据。由于数据是我们自己生成，所以知道真实的参数

，以下代码中用theta_real变量表示。显然，在真实数据中，我们并不知道参数的真实值，而是要将其估计出来。

np.random.seed(123)
n_experiments = 4
theta_real = 0.35
data = stats.bernoulli.rvs(p=theta_real, size=n_experiments)
print(data)
array([1, 0, 0, 0])

模型描述

现在有了数据，需要再指定模型。回想一下，模型可以通过指定似然和先验的概率分布完成。对于似然，我们可以用参数分别为

和

的二项分布来描述，对于先验，我们可以用参数为

的beta分布描述。这个beta分布与[0,1]区间内的均匀分布是一样的。我们可以用数学表达式描述如下：

这个统计模型与PyMC3的语法几乎一一对应。第1行代码先构建了一个模型的容器，PyMC3使用with语法将所有位于该语法块内的代码都指向同一个模型，你可以把它看作是简化模型描述的“语法糖”，这里将模型命名为our_first_model。第2行代码指定了先验，可以看到，语法与数学表示很接近。我们把随机变量命名为

，需要注意的是，这里变量名与Beta函数的第1个参数名一样；保持相同的名字是个好习惯，这样能避免混淆。然后，我们通过变量名从后验采样中提取信息。这里变量

是一个随机变量，我们可以将该变量看做是从某个分布（在这里是beta分布）中生成数值的方法而不是某个具体的值。第3行代码用跟先验相同的语法描述了似然，唯一不同的是我们用observed变量传递了观测到的数据，这样就告诉了PyMC3我们的似然。其中，data可以是一个Python列表或者Numpy数组或者Pandas的DataFrame。这样我们就完成了模型的描述。

with pm.Model() as our_first_model:
 theta = pm.Beta('theta', alpha=1, beta=1)
 y = pm.Bernoulli('y', p=theta, observed=data)

按下推断按钮

对于抛硬币这个问题，后验分布既可以从分析的角度计算出来，也可以通过PyMC3用几行代码从后验的采样中得到。代码中的第1行，调用了find_MAP函数，该函数调用SciPy中常用的优化函数尝试返回最大后验（Maximum a Posteriori，MAP）。调用find_MAP是可选的，有时候其返回值能够为采样方法提供一个不错的初始点，不过有时候却并没有多大用，因此大多数时候会避免使用它。然后，下一行定义了采样方法。这里用的是Metropolis-Hastings算法，函数名取了简写Metropolis。PyMC3可以让我们将不同的采样器赋给不同的随机变量；眼下我们的模型只有一个参数，不过后面我们会有多个参数。我们也可以省略该行，PyMC3会根据不同参数的特性自动地赋予一个采样器，例如，NUTS算法只对连续变量有用，因而不能用于离散的变量，Metropolis算法能够处理离散的变量，而另外一些类型的变量有专门的采样方法。总的来说，我们可以让PyMC3为我们选一个采样方法。最后一行是执行推断，其中第1个参数是采样次数，第2个和第3个参数分别是采样方法和初始点，可以看到这两个参数是可选的。

 start = pm.find_MAP()
 step = pm.Metropolis()
 trace = pm.sample(1000, step=step, start=start)

这样，只需要几行代码我们就完成了整个模型的描述和推断。感谢PyMC3的开发者们为我们提供了这么棒的库。

诊断采样过程

现在我们根据有限数量的样本对后验做出了近似，接下来要做的第一件事就是检查我们的近似是否合理。我们可以做一些测试，有些是可视化的，有些是定量的。这些测试会尝试从样本中发现问题，但并不能证明我们得到的分布是正确的，它们只能提供证据证明样本看起来是合理的。如果我们通过样本发现了问题，解决办法有如下几种。

• 增加样本次数。
• 从样本链（迹）的前面部分去掉一定数量的样本，称为 老化 （Burn-in）。在实践中，MCMC方法通常需要经过一段时间的采样之后，才得到真正的目标分布。老化在无限多次的采样中并不是必须的，因为这部分并没有包含在马尔科夫理论中。事实上，去掉前面部分的样本只不过是在有限次采样中提升结果的一个小技巧。需要注意，不要被数学对象和数学对象的近似弄糊涂了，球体、高斯分布以及马尔科夫链等数学对象只存在于柏拉图式的想象世界中，并不存在于不完美但却真实的世界中。
• 重新参数化你的模型，也就是说换一种不同但却等价的方式描述模型。
• 转换数据。这么做有可能得到更高效的采样。转换数据的时候需要注意对结果在转换后的空间内进行解释，或者再重新转换回去，然后再解释结果。

收敛性

通常，我们要做的第一件事就是查看结果长什么样，traceplot函数非常适合该任务：

burnin = 100
chain = trace[burnin:]
pm.traceplot(chain, lines={'theta':theta_real});

对于未观测到的变量，我们得到了两幅图。左图是一个 核密度估计 （Kernel Density Estimation，KDE）图，可以看做是平滑之后的直方图。右图描绘的是每一步采样过程中得到的采样值。注意图中红色的线表示变量theta_real的值。

在得到这些图之后，我们需要观察什么呢？首先，KDE图看起来应该是光滑的曲线。通常，随着数据的增加，根据中心极限定理

，参数的分布会趋近于高斯分布。当然，这并不一定是正确的。右侧的图看起来应该像白噪声，也就是说有很好的 混合度（mixing） ，我们看不到任何可以识别的模式，也看不到向上或者向下的曲线，相反，我们希望看到曲线在某个值附近震荡。对于多峰分布或者离散分布，我们希望曲线不要在某个值或区域停留过多时间，我们希望看到采样值在多个区间自由移动。此外，我们希望迹表现出稳定的相似性，也就是说，前10%看起来跟后50%或者10%差不多。再次强调，我们不希望看到规律的模式，相反我们期望看到的是噪声。下图展示了一些迹呈现较好混合度（右侧）与较差混合度（左侧）的对比。

如果迹的前面部分跟其他部分看起来不太一样，那就意味着需要进行老化处理，如果其他部分没有呈现稳定的相似性或者可以看到某种模式，那这意味着需要更多的采样，或者需要更换采样方法或者参数化方法。对于一些复杂的模型，我们可能需要结合使用前面所有的策略。

PyMC3可以实现并行地运行一个模型多次，因而对同一个参数可以得到多条并行的迹。这可以通过在采样函数中指定njobs参数实现。此时使用traceplot函数，便可在同一幅图中得到同一个参数的所有迹。由于每组迹都是相互独立的，所有的迹看起来都应该差不多。除了检查收敛性之外，这些并行的迹也可以用于推断，我们可以将这些并行的迹组合起来提升采样大小而不是扔掉多余的迹：

with our_first_model:
 step = pm.Metropolis()
 multi_trace = pm.sample(1000, step=step, njobs=4)

burnin = 0
multi_chain = multi_trace[burnin:]
pm.traceplot(multi_chain, lines={'theta':theta_real});

一种定量地检测收敛性的方法是 Gelman-Rubin 检验。该检验的思想是比较不同迹之间的差异和迹内部的差异，因此，需要多组迹来进行该检验。理想状态下，我们希望得到

。根据经验，我们认为如果得到的值低于1.1，那么可以认为是收敛的了，更高的值则意味着没有收敛：

pm.gelman_rubin(multi_chain)
{'theta': 1.0074579751170656, 'theta_logodds': 1.009770031607315}

我们还可以用forestplot函数将

和每个参数的均值、50%HPD和95%HPD可视化地表示出来：

pm.forestplot(multi_chain, varnames=['theta']);

函数summary提供了对后验的文字描述，它可以提供后验的均值、标准差和HPD区间：

pm.summary(multi_chain)
theta:
 Mean SD MC Error 95% HPD interval
 -------------------------------------------------------------------
 0.339 0.173 0.006 [0.037, 0.659]
Posterior quantiles:
 2.5 25 50
75 97.5
 |--------------|==============|==============|--------------|
 0.063 0.206 0.318
0.455 0.709

此外，df_summary函数会返回类似的结果，不过类型是Pandas中的DataFrame：

pm.df_summary(multi_chain)

其中，返回值之一是mc_error，这是对采样引入误差的估计值，该值考虑的是所有的采样值并非真的彼此独立。mc_error是迹中不同块的均值的标准差，每一块是迹中的一部分：

该误差值显然低于我们结果的准确度。由于采样方法是随机的，每次重跑的时候，summary或者df_summary返回的值都会不同，不过没关系，mc_error的值应该是相似的，如果返回的值有很大不同，则说明我们可能需要更多的样本。

自相关

最理想的采样应该不会是自相关的，也就是说，某一点的值应该与其他点的值是相互独立的。在实际中，从MCMC方法（特别是Metropolis-Hastings）中得到的采样值是自相关的。由于参数之间的相互依赖关系，可能模型会导致更多的自相关采样。PyMC3有一个很方便的函数用来描述自相关。

pm.autocorrplot(chain)

该图显示了采样值与相邻连续点（最多100个）之间的平均相关性。理想状态下，我们不会看到自相关性，实际中我们希望看到自相关性降低到较低水平。参数越自相关，要达到指定精度的采样次数就需要越多，也就是说，自相关性不利于降低采样次数。

有效采样大小

一个有自相关性的采样要比没有自相关性的采样所包含的信息量更少，因此，给定采样大小和采样的自相关性之后，我们可以尝试估计出该采样的大小为多少时，该采样没有自相关性而且包含的信息量不变，该值称为有效采样大小。理想情况下，两个值是一模一样的；二者越接近，我们的采样效率越高。有效采样大小可以作为我们的一个参考，如果我们想要估计出一个分布的均值，我们需要的最小采样数至少为100；如果想要估计出依赖于尾部分布的量，比如可信区间的边界，那么我们可能需要1000到10000次采样。

pm.effective_n(multi_chain)['theta']
667

显然，提高采样效率的一个方法是换一个更好的采样方法；另一个办法是转换数据或者对模型重新设计参数，此外，还有一个常用的办法是对采样链压缩。所谓压缩其实就是每隔 k 个观测值取一个，在Python中我们称为切片。压缩会降低自相关性，但代价是同时降低了样本量。因此，实际使用中通常更倾向于增加样本量而不是切片。不过有时候，压缩会很有用，比如降低存储空间。如果仍不能避免高自相关性，我们就只能算出更长的采样链，模型中的参数很多的话，存储量会是个问题。而且，我们可能还会对后验做一些计算量很大的后处理，此时在自相关性尽可能小的前提下，采样数量的大小就显得尤为重要。

总结

目前为止，所有的诊断测试都是经验性而非绝对的。实际使用中，我们会先运行一些测试，如果看起来没什么问题，我们就继续往下分析。如果发现了一些问题，就需要回过头解决它们，这也是建模过程的一部分。要知道，进行收敛性检查并非贝叶斯理论的一部分，由于我们是用数值的方式在计算后验，因而这只是贝叶斯实践过程中的一部分。