pymc3 MCMC应用【转载】-CSDN博客

转载：
https://blog.csdn.net/qq_44009891/article/details/106069563

第一部分编程准备
贝叶斯思维：和更传统的统计推断不同，贝叶斯推断会保留不确定性，在贝叶斯派的世界观中，概率是被解释为我们对一件事情发生的相信程度或者说信心（飞机事故，总统选举）。

需注意的是，我们每个人都可以给事件赋概率值，而不是存在某个唯一的概率值，因为不同的人拥有不同的信息，因此他们对同一事件发生的信心也可以有不同的值，但这些不同并不说明其他人是错误的。

飞机事故：综合某航空公司过去十年出现飞机事故的次数，来推断该公司出现飞机事故的概率；
总统选举：根据候选人的声望、民意等等，推断他能当上总统的概率。
贝叶斯推断的工作方式：我们会随着新的证据不断更新之前的信念，但很少做出绝对的判断，除非所有其他的可能都被一一排除（代码 bug）。

代码 bug：“我的代码通常会有 bug”，这是我的先验信息，我现在要来验证是否有 bug，于是我测试了一个简单的案例，这个案例通过后，我再用一个稍微复杂的测试案例，再次通过了，接下来更难的测试案例也通过了，通过这一连串的测试，我开始觉得我的这段代码很可能已经没有 bug 了，于是我得到我的后验概率，“我的这段代码有 80% 的可能没有 bug”。
相关的 Python 第三方库安装：

Ipython
Numpy
Scipy
Pandas
Matplotlib
Seaborn
PyMC3
Jupyter Notebook 安装及配置使用：

Jupyter Notebook 安装及常见问题：https://www.cnblogs.com/justaman/p/11282655.html
Jupyter Notebook 更改默认浏览器：https://cloud.tencent.com/developer/article/1420759
Jupyter Notebook 更改默认目录：https://blog.csdn.net/nico2333/article/details/84186063
Jupyter Notebook 更改默认目录常见问题：https://blog.csdn.net/lmshan123/article/details/92109976
第二部分概率编程入门实例——贝叶斯点估计
我们来以《贝叶斯统计》（韦来生著）这本书上的一个习题作为我们入门的实例，因为学习 PyMC3 必然是要落地的，通过书上的问题能帮我们迅速拉进和 PyMC3 的距离，这里选出第四章贝叶斯统计推断中的关于贝叶斯点估计的习题来看一下。

1 观察问题

这是第四章习题中的第 10 题，这里把题目截图放在这里，我们看一看这个题目，它有五个样本 X1 到 X5，而且已经有了观测值然后这些样本服从指数分布 Exp(1 / θ)，而 θ 的先验分布是逆伽马分布 Γ-1(10，100)，然后我们要求什么呢，我们要求参数 θ 的广义最大似然估计（这里后验均值估计以及后验方差就不求了），好，题目是这个样子，我们现在使用编程来解决这个问题。

2 确定观测数据
第二步，导入第三方库，以及确定观测数据：

import pymc3 as pm

X = [5, 12, 14, 10, 12]
1
2
3
3 定义模型与变量
接着第三步，定义模型，我们使用 pm.Model() 方法，以及 with 上下文管理器将该代码块中的变量都添加到 my_first_model 模型，注意，这个 with pm.Model() as my_first_model 不同于其他的有关 with 的使用方法，my_first_model 这个模型不会在退出 with 时自动关闭，而是会一直保存，作为连接整个代码程序的枢纽：

with pm.Model() as my_first_model:

theta = pm.InverseGamma('theta', alpha=10, beta=100)
X_obs = pm.Exponential('X_obs', lam=1 / theta, observed=X)

1
2
3
4
然后我们定义 θ 随机变量，它是逆伽马分布，PyMC3 也提供了对应的 API —— InverseGamma 类，这是对应的网址：https://docs.pymc.io/api/distributions/continuous.html#pymc3.distributions.continuous.InverseGamma，我们点进去看一下，界面如下面截图所示，我们可以发现 InverseGamma 类有很多参数，但是只有前两个参数 alpha 和 beta 是我们需要传入的，它们分别对应分布中的 α 和 β，而后面的这些参数都是可选或可替代的，其中，mu 参数表示逆伽马分布的期望，sigma 参数对应分布的标准差，因为 α 和 β 与分布的期望和标准差之间可以通过计算进行转换，因此当传入参数 alpha 和 beta 后，就不必再传入参数 mu 和 sigma，同理，反过来也是一样，所以称参数 mu 和 sigma 是可供替代的参数，还有一个参数 sd，这个参数其实也代表标准差，不过当它有传入进来的值时即不为 None 时，优先选择 sd 而不是 sigma 作为标准差。

还有一点需要注意的是，InverseGamma 类还有一个隐藏的必填参数，那就是 name，我们可以看到如果仅仅传入 alpha 和 beta 参数运行时会出现报错：

报错信息告诉我们 new() 缺少一个参数，但是 InverseGamma 类的 init() 函数明明没有 name 这样一个必选参数啊，为什么会出现这样的报错信息，其实这是因为 InverseGamma 类是一个多重继承的类，它最开始的父类是 Distribution 类，这个类有一个 new 函数，需要一个必填的字符串型的 name 参数，而 new 函数会在 InverseGamma 类实例创建之前被调用，它的任务就是创建实例，此时缺少 name 参数 new 函数就会报错，关于 new 函数的具体细节，大家可以参照这一篇博文：https://blog.csdn.net/qq_41020281/article/details/79638370。

因此，我们需要传入一个合适的名字，一般选择变量名作为 name，这样将方便从输出结果中进行查找对应的变量，如下所示：

刚才说了，InverseGamma 类是一个继承类，但是它也有着自己的类函数，分别是 logp() 和 random()，实例化类时我们传入了 alpha 参数和 beta 参数就有了分布，而 logp() 函数可以帮助我们计算分布的对数似然，这里的似然就是我们的分布密度，它的参数 value 就是分布密度公式中的 x，还有 random() 函数可以帮助我们从分布中抽取随机的样本。

不过在这里我还遇到了一个很大的坑，之前关于这两个函数我还没有理解透，但在和老师和同学交流之后，稍稍有了一点感悟，大家如果有兴趣，可以移步至知识扩展。

刚刚讲了一下关于 InverseGamma 类的内容，现在来重新看一下写在上面的这一步的代码，让我们回到初始问题上来：

with pm.Model() as my_first_model:

theta = pm.InverseGamma('theta', alpha=10, beta=100)
X_obs = pm.Exponential('X_obs', lam=1 / theta, observed=X)

1
2
3
4
确定了 θ 之后，我们要来确定需要拟合的变量，即 X_obs，它服从指数分布 Exp(1 / θ)，因此传入的参数是 1 / theta，而参数 observed=X 是告诉模型这个变量的值是已经被观测到了的，不会被拟合算法改变，而观测值就是 X，在这里，X 是一开始我们定义的列表。

这里也给出 PyMC3 中的指数分布类 Exponential 所对应的网址：
https://docs.pymc.io/api/distributions/continuous.html#pymc3.distributions.continuous.Exponential，该类和刚才讲的 InverseGamma 类差不多，没有太大变化，可以说 PyMC3 所提供的分布类基本上都可以套用上述的理解方式去学习，不过这个现象可能是表面上的，其中还是有许多坑的，大家慎踩。

4 统计推断
定义了模型之后，我们就要开始对模型进行统计推断了，这里选择 MAP 作为推断方法，那什么是 MAP 呢，MAP（maximum a posteriori probability estimate）中文名称叫极大似然估计，其实它就是《贝叶斯统计》这本书上所讲的广义最大似然估计，即后验众数估计，PyMC3 也提供了对应的函数 find_MAP()，不过它的计算方式和我们书上说的不同，因为书上是通过分析计算的方式推断出参数的估计值，但是这种方式一般只适合于标准分布，而实际上我们碰到的很多分布都是非标准分布，很难或无法通过分析计算进行推断，这个时候就需要通过 find_MAP() 所采用的数值优化的方法得出参数的估计值。

map_estimate = pm.find_MAP(model=my_first_model, start={‘theta’:10})
1
实现极大似然估计这一步的代码如上所示，非常简单，不过大家可能对于 find_MAP() 还有一些疑问，请大家移步至我的另一篇博文：PyMC3 API 解读（一）—— find_MAP() 函数。

5 代码整理与输出
到了这里，问题的基本代码讲解完毕，我们现在将它们整合到一起，然后输出最终的结果，这也顺便帮助我们回顾一下刚才讲解的整个过程：

import pymc3 as pm

X = [5, 12, 14, 10, 12]

with pm.Model() as my_first_model:
theta = pm.InverseGamma(‘theta’, alpha=10, beta=100)
X_obs = pm.Exponential(‘X_obs’, lam=1 / theta, observed=X)

map_estimate = pm.find_MAP(model=my_first_model, start={‘theta’:10})

print(map_estimate)
1
2
3
4
5
6
7
8
9
10
11
最后的输出结果如下，theta 的点估计值是 9.56249988：

我们也会发现除 theta 外还有一个输出值—— theta_log__，它就是 theta 经过代数变换后的变量，我们也可以验证一下，发现确实是这样，9.56249988 经过 log 变换，得到了 theta_log__的值—— 2.2578491866036345：

然后我们再比对一下《贝叶斯统计》这本书上该题的答案，我们可以看到答案是 9.653，这和我们的点估计结果非常接近，说明我们的代码运行无误，问题得到解决：

到这里，这一部分的代码就讲解完毕了，我们通过编程的手段来解决了一个实际问题，虽然整个部分的代码只有七行，但是它们内含的门道特别绕，需要去一点点地推敲琢磨，敲完这段代码或许只需要几分钟，但理解它们可能需要几个小时甚至更多时间。

6 知识扩展
关于 InverseGamma 类的这两个类函数—— logp() 和 random()，之前我一直觉得不太懂，因为它们和其他的分布类比如说 Normal 正态分布类的同名类函数有着一定的差别，当时我怀疑这可能是因为逆伽马分布有着形状参数，比如说 α 就是形状参数（关于形状参数大家可以参照这篇博文：https://blog.csdn.net/weixin_30782871/article/details/95108032），因此它在 PyMC3 中的存在和运作形式和正态分布这样没有形状参数的分布有一些不同，但是在经过与老师和同学的讨论后，我发现我的猜想是错误的，在这里我将讲一讲我的疑惑和感悟。

我们先通过打印 InverseGamma 类和 Normal 类的实例化对象的类型来观察一下，发现它们的实例化对象在模型中被定义为不同的类型：
实际上，为什么会产生这样的不同是因为逆伽马分布对应的函数是有界函数，而正态分布对应的函数是无界函数，有界的意思是指分布密度函数中的 x 被限定再一个范围内，比如说逆伽马分布要求 x 大于 0，因此它对应的是有界函数，而正态分布中的 x 可以取任何一个值，从负无穷取到正无穷都没有问题，因此它对应的是无界函数，所以它对应的类型为 pymc3.model.FreeRV，表示不受限制的随机变量的意思（大家需要注意的是，这里的有界函数和我们高中学的有界函数并不相同，官网的说法是 Automatic transforms of bounded RVs ——有界随机变量的自动（代数）变换，所以只是在这里使用有界函数的称呼，该称呼可能对大家造成一定程度的误解，请大家谅解）。

一个例子可能还不能说明这个现象，我们再选取两个有界函数对应的分布——指数分布和贝塔分布，指数分布要求 x 大于 0，而贝塔分布要求 x 大于 0 小于 1，我们使用它们分别初始化随机变量 RV 然后打印它们的类型，发现它们确实是 pymc3.model.TransformedRV 类型，这进一步佐证了我们的论述：

PyMC3 在生成变量的过程中为了方便以后的抽样，它会将有界变量自动转换为无界变量，而可以做到这一切的幕后黑手其实就是参数 transform，那些有界函数的对应的分布在 PyMC3 中自带有一个 transform 参数，在模型初始化有界随机变量时，它就会被自动调用，将有界函数变换为无界函数，而对应的类型也会变成 pymc3.model.TransformedRV，表示被代数变换后的变量的意思，我们来尝试一下，看到底是不是这样：

从上图我们可以看到，当 transform 参数被设置为 None 时，将不会发生从有界到无界的变换，打印出的类型也还是 pymc3.model.FreeRV，也可以通过调用 logp() 函数来得到对应的对数似然，然而，当不对 transform 参数进行设置时，对应的随机变量就会自动转变为 pymc3.model.TransformedRV，而且还会出现 ‘TransformedRV’ object has no attribute ‘logp’ 的报错信息，因为在这种情况下，虽然我们表面上没有传参给 transform，但实际上它自己会自动向 transform 参数传入一个函数，接着调用该函数帮助逆伽马分布完成从有界到无界的华丽逆转。

那 transform 参数接收到的这个函数到底是什么呢，我们可以通过两种方式查看：一种是打印模型的 FreeRV 列表，因为从无界变换过来的变量会自动变成 FreeRV，在这里就是 x_inverse_gamma 是 TransformedRV，而经过 log 变换后的 x_inverse_gamma_log__ 就是 FreeRV（注意，x_inverse_gamma 依然会保存在模型当中）；另一种方法就是直接打印 transformation.name 从而得到进行代数变换的函数名（在这里，我又产生了新的疑惑，为什么逆伽马分布经过 log 变换后会从有界变成无界呢）。

第三部分概率编程基础实例——线性回归
1 观察问题
设 Y 1 , Y 2 , ⋯ , Y n 为从正态总体 N ( μ , σ 2 ) 中抽取的随机样本 , σ 代表观察误差 , 而 μ = α + β 1 X 1 + β 2 X 2 , α 是偏置项 , β i 是随机变量 X i 的系数 , 其中 X 1 和 X 2 都服从标准正态分布 N ( 0 , 1 ) , 且已经分别拥有 N 个观测值 : { x 11 , x 12 , ⋯ , x 1 n } , { x 21 , x 22 , ⋯ , x 2 n } , 并且我们已经知道所有参数的先验分布 , 其中 , α 和 β 的先验为零均值的正态先验 , 而 σ 的先验分布为半正态分布 , 如下所示 , 现在我们根据 X 1 和 X 2 的观测值和未知参数得到了 Y 的 N 个观测值 : { y 1 , y 2 , ⋯ , y n } , 求各个未知参数的点估计值。
设 Y1,Y2,⋯,Yn 为从正态总体 N(μ,σ2) 中抽取的随机样本, σ 代表观察误差, 而 μ=α+β1X1+β2X2, α 是偏置项, βi 是随机变量 Xi 的系数, 其中 X1 和 X2 都服从标准正态分布 N(0,1), 且已经分别拥有 N 个观测值:{x11,x12,⋯,x1n}, {x21,x22,⋯,x2n}, 并且我们已经知道所有参数的先验分布, 其中, α 和 β 的先验为零均值的正态先验, 而 σ 的先验分布为半正态分布, 如下所示, 现在我们根据 X1 和 X2 的观测值和未知参数得到了 Y 的 N 个观测值:{y1,y2,⋯,yn}, 求各个未知参数的点估计值。
设 Y1,Y2,⋯,Yn 为从正态总体 N(μ,σ2) 中抽取的随机样本, σ 代表观察误差, 而 μ=α+β1X1+β2X2, α 是偏置项, βi 是随机变量 Xi 的系数, 其中 X1 和 X2 都服从标准正态分布 N(0,1), 且已经分别拥有 N 个观测值:{x11,x12,⋯,x1n}, {x21,x22,⋯,x2n}, 并且我们已经知道所有参数的先验分布, 其中, α 和 β 的先验为零均值的正态先验, 而 σ 的先验分布为半正态分布, 如下所示, 现在我们根据 X1 和 X2 的观测值和未知参数得到了 Y 的 N 个观测值:{y1,y2,⋯,yn}, 求各个未知参数的点估计值。

设 Y
1

,Y
2

,⋯,Y
n

为从正态总体 N(μ,σ
2
) 中抽取的随机样本, σ 代表观察误差, 而 μ=α+β
1

X
1

+β
2

X
2

, α 是偏
置项, β
i

是随机变量 X
i

的系数, 其中 X
1

和 X
2

都服从标准正态分布 N(0,1), 且已经分别拥有 N 个观测值:{x
11

,
x
12

,⋯,x
1n

}, {x
21

,x
22

,⋯,x
2n

}, 并且我们已经知道所有参数的先验分布, 其中, α 和 β 的先验为零均值的正态先
验, 而 σ 的先验分布为半正态分布, 如下所示, 现在我们根据 X
1

和 X
2

的观测值和未知参数得到了 Y 的 N 个观测
值:{y
1

,y
2

,⋯,y
n

}, 求各个未知参数的点估计值。

α ∼ N ( 0 , 100 ) β i ∼ N ( 0 , 100 ) σ ∼ ∣ N ( 0 , 1 ) ∣ \alpha \sim N(0, 100) \ \beta_i \sim N(0, 100) \ \sigma \sim |N(0, 1)|
α∼N(0,100)
β
i

∼N(0,100)
σ∼∣N(0,1)∣

该实例是 PyMC3 官网上 getting_started 页面的一个例子，同时也参考了一位大佬的博文，网址分别是：https://docs.pymc.io/notebooks/getting_started 和 https://blog.csdn.net/jackxu8/article/details/71080865，而本题是我针对该实例模仿《贝叶斯统计》这本书上的习题格式抽象出来的问题，以此帮助大家理解该实例。

观察问题，我们发现，其实这就是一个线性回归预测的问题，简单来说，就是现在有函数 y = α + β1x1 + β2x2 + σx3，我们从四维空间（X1, X2, X3, Y）中抽取了 N 个符合该函数的样本点，然后我们要根据这 N 个样本点进行训练，最后得到 α、β1、β2 和 σ 这四个未知参数的估计值，这样就得到了一个关于函数 y = α + β1x1 + β2x2 + σx3 的预测模型，然后我们就可以拿该模型去预测其他样本点的 Y 了。

这其实和平面上的线性回归预测问题没啥区别，平面上的函数是 y = b + wx，然后根据已知的（x, y) 样本点预测 b 和 w，而我们的函数是 y = α + β1x1 + β2x2 + σx3，然后根据已知的样本点预测这四个未知参数。

讲到这里，大家可能对刚才的突然冒出来的函数 y = α + β1x1 + β2x2 + σx3 比较奇怪，会疑惑这个函数是如何得到的，其中的 x3 又是怎么来的，接下来解释一下，因为 Y 是服从正态分布 N(μ, σ2) 的，因此可以得出下面的公式：

( Y − μ ) σ ∼ N ( 0 , 1 ) \frac {(Y - \mu)} \sigma \sim N(0, 1) \
σ
(Y−μ)

∼N(0,1)

而我们假设有一个随机变量 X3，它服从标准正态分布 N(0, 1)，则上面的公式可以写成：

( Y − μ ) σ = X 3 ∼ N ( 0 , 1 ) \frac {(Y - \mu)} \sigma = X_3 \sim N(0, 1) \
σ
(Y−μ)

=X
3

∼N(0,1)

接着进行变换，就可以得到在四维空间（X1, X2, X3, Y）上的线性函数：

Y = μ + σ X 3 = α + β 1 X 1 + β 2 X 2 + σ X 3 Y = \mu + \sigma X_3 = \alpha + \beta_1 X_1 + \beta_2 X_2 + \sigma X_3
Y=μ+σX
3

=α+β
1

X
1

+β
2

X
2

+σX
3

因此，我们的问题表面上只有 X1 和 X2 两个自变量，但实际上其中还隐藏着一个自变量 X3，找出了这个 X3，等于我们就得到了一个标准的线性函数：y = α + β1x1 + β2x2 + σx3，然后对这个函数进行建模和训练就可以解决本题了。

不过，需注意的是，这里的 X3 提出来是帮助大家理解，在实际使用 PyMC3 进行编程的时候，这个 X3 是隐式的，把它用于模拟一下样本观测值就可以了。

注：半正态分布就是正态分布的一半，即当 x < 某一阈值时，概率密度为零，而当 x >= 某一阈值时，概率密度又符合以该阈值为均值，方差不定的某一正态分布，这里设置阈值为 0，具体信息可以参照知乎上的这个回答：https://www.zhihu.com/question/264514748。

正态分布：
半正态分布：

2 模拟观测数据
因为题目没有直接提供样本的观测值，所以我们需要自己模拟观测数据，那么我们首先就要模拟参数 α, βi 和 σ，并根据参数对分布进行抽样，用抽样出来的样本来模拟观测值，如下：

import numpy as np

设置随机数种子

np.random.seed(123)

设置参数

true_alpha = 1
true_sigma = 1
true_beta = [1, 2.5]

设置样本数

N = 100

从标准正态分布中抽取 N 个样本，组成列表

X1 = np.random.randn(N)
X2 = np.random.randn(N)
X3 = np.random.randn(N)

通过对标准正态分布的观测值进行变换得到 Y 的观测值

Y = true_alpha + true_beta[0]*X1 + true_beta[1]X2 + true_sigmaX3
print(“X1:”, X1)
print(“X2:”, X2)
print(“Y :”, Y)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
打印模拟出来的观测数据：

X1: [-1.0856306 0.99734545 … 0.37940061 -0.37917643]
X2: [ 0.64205469 -1.97788793 … -0.34126172 -0.21794626]
Y : [ 2.22281624 -3.54547971 … 0.68679074 0.89571852]
1
2
3
得到样本的观测值后，通过绘制散射图来观察随机变量 X1、X2 和 Y之间的关系：

import matplotlib.pyplot as plt
%matplotlib inline

绘制散射图，其中参数 1 和 2 分别代表子图的行数和列数，一共有 1 * 2 个子图像，通过 figszie 可以设置子图的宽度和高度

fig1, ax1 = plt.subplots(1, 2, figsize=(10,5))
ax1[0].scatter(X1, Y)
ax1[0].set_xlabel(‘X1’)
ax1[0].set_ylabel(‘Y’)
ax1[1].scatter(X2, Y)
ax1[1].set_xlabel(‘X2’)
ax1[1].set_ylabel(‘Y’)
1
2
3
4
5
6
7
8
9
10
11
我们发现 X2 和 Y 呈现出正相关的关系，并且相关度比较高，但是 X1 和 Y 的相关度较低，这是因为我们刚才设置它们的系数 β1 和 β2 时，β1 设置为了 1，而 β2 设置为了 2.5，把他们看作权重的话，那很明显，X2 的权重相对更大一点，因此对 Y 的影响也更大一点。

到这里，得到了样本观测值，并且也大概了解了随机变量 X1、X2 和 Y之间的关系，那么现在就请立马忘掉刚才设置的四个参数，假装不知道它们是多少，假装我们的样本观测值是题目给出来的，因为现在我们有了观测数据，题目完整了，之后就需要估计未知参数们的值了。

3 定义模型与变量
这一步在前面的内容中讲的应该比较清楚了，先定义模型，然后再在该模型定义对应的变量，不清楚可以看看前面的或者通过注释了解一下：

import pymc3 as pm

定义一个新的模型对象，这个对象是模型中随机变量的容器

with pm.Model() as my_second_model:

# 定义两个具有正态分布先验和一个具有半正态分布先验的随机性随机变量
alpha = pm.Normal('alpha', mu=0, sd=10) # 正态分布
beta = pm.Normal('beta', mu=0, sd=10, shape=2) # 正态分布
sigma = pm.HalfNormal('sigma', sd=1) # 半正态分布

# 定义一个确定性随机变量，确定性指这个变量的值完全由右端值确定
mu = alpha + beta[0]*X1 + beta[1]*X2

# 定义 Y 的观测值变量
Y_obs = pm.Normal('Y_obs', mu=mu, sd=sigma, observed=Y)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
4 模型拟合
模型拟合这一步其实就等价于贝叶斯的统计推断，至于为什么把这一步命名为模型拟合，一是因为官网上的称呼是 Model fitting，我们应该尊重权威，二是因为模型拟合确实能够指出接下来讲的一系列代码的关键点—— PyMC3 使用迭代优化的方式对未知参数进行拟合来得到各个参数的最优估计值，而不是通过分析计算的方式得到。

第一种方法 MAP 方法
关于 MAP 极大似然估计方法，前面也讲过了，具体的解读我也独立出来到另一篇从博文 PyMC3 API 解读（一）—— find_MAP() 函数当中了，这里只稍微点一下：

使用不同的优化方法找到参数的极大后验估计

map_estimate1 = pm.find_MAP(model=my_second_model) # 默认使用 L-BFGS—B 优化算法
map_estimate2 = pm.find_MAP(model=my_second_model, method=“BFGS”) # 使用 BFGS 优化算法
1
2
3
打印不同的优化方法找到的 MAP 估计值：

打印估计值

print(map_estimate1)
print(map_estimate2)
1
2
3
得到的输出结果为：

可以看到通过两种算法得到的估计值都和我们原本设定的值较为接近，说明使用 MAP 点估计的方法能够解决本问题。

不过需注意的是，MAP 常常找到的最优估计值只是局部最优的，我们一般使用接下来要讲的第二种方法—— MCMC 采用算法来拟合模型的全局最优，当然，可能在这一题，我们暂时看不出它们效果的差异，因为它们的估计值都和初始设定的较为接近，难以判断孰优孰劣，不过在现实世界的很多后验实验中，它们的差距是比较大的。

第二种方法 MCMC 采样方法
对于 MCMC 采用算法，大家学过贝叶斯的应该都比较熟悉了，但如果没有深入理解的话，终究只是盲人摸象，所以希望大家在使用 PyMC3 的过程中，通过实践不断地加深自己对 MCMC 算法的理解，同时这句话也送给我自己，因为我也没有理解清楚 MCMC 算法，今后希望和大家一起交流。

PyMC3 中提供了一些内置的 MCMC 采样算法（官方把采样算法称为 step method），比如 NUTS、Metropolis、Slice、HamiltonianMC 和 BinaryMetropolis 等，采样算法可以由 PyMC3 自动指定也可以手动进行指定，自动指定是根据模型中的变量类型决定的，手动指定优先，可覆盖自动指定，如下为自动指定的采样算法：

二值变量：自动指定 BinaryMetropolis；
离散变量：自动指定 Metropolis；
连续变量：自动指定 NUTS。
（1）MCMC 采样
在这一小步中，我们使用 NUTS 和 Slice 两种不同的采样算法对诸参数 α, βi 和 σ 的后验分布进行采样，关于 NUTS 采样算法可参照：https://zhuanlan.zhihu.com/p/59473302，关于 Slice 采样算法可参照：https://blog.csdn.net/baihaoli123/article/details/53886077，对于这两个采用算法不了解也没关系，我们可以在日后的实践过程中慢慢了解，我反而觉得现在掌握整个过程的脉络是更为重要的，我们先看一下这一小步的代码，然后再对它进行解读：

with my_second_model:
# 自动指定 NUTS 采样算法，并对后验分布进行 6000 次采样，其中，预烧期为 1000 次，正式采样 5000 次
trace1 = pm.sample(draws=5000, tune=1000)

# 使用 MAP 获得初始值
start = pm.find_MAP(method="BFGS")
# 对于 sigma 指定使用 Slice 采样算法，而其他两个参数 alpha 和 beta 的采样方法由模型自动指定为 NUTS 采样算法
step = pm.Slice(vars=[sigma])
# 对后验分布进行 6000 次采样，其中，预烧期为 1000 次，正式采样 5000 次
trace2 = pm.sample(draws=5000, tune=1000, step=step, start=start)

1
2
3
4
5
6
7
8
9
10
上面的代码中，my_second_model 是之前已经定义好的模型，也是连接整个代码程序的枢纽，然后，就可以在 my_second_model 这个模型中对诸参数 α, βi 和 σ 的后验分布进行采样，这会用到 PYMC3 当中的一个关键函数 sample()，在我的另一篇博文 PyMC3 API 解读（二）—— sample() 函数当中，我对它进行了较为详细的解释。

这段代码运行时的进度条如下：

（2）打印采样结果
我们在上一小步通过 MCMC 采样方法对诸参数 α, βi 和 σ 的后验分布进行了采样，而收集到的采样值都存储在 pymc3.backends.base.MultiTrace 类的对象中，按照采样值获得的先后顺序存储，pymc3.backends.base.MultiTrace 类的对象是一个类似于字典的对象，大家可以直接把它当作字典来用，现在我们来打印采样结果：

print(“trace1:”)
print(“alpha”, trace1[‘alpha’].shape)
print(trace1[‘alpha’], “\n”)
print(“beta”, trace1[‘beta’].shape)
print(trace1[‘beta’], “\n”)
print(“sigma”, trace1[‘sigma’].shape)
print(trace1[‘sigma’], “\n\n”)

print(“trace2:”)
print(“alpha”, trace2[‘alpha’].shape)
print(trace2[‘alpha’], “\n”)
print(“beta”, trace2[‘beta’].shape)
print(trace2[‘beta’], “\n”)
print(“sigma”, trace2[‘sigma’].shape)
print(trace2[‘sigma’], “\n\n”)

采样结果如下：

（3）绘制后验采样的趋势图
在这一小步，我们通过 traceplot() 函数来绘制 MCMC 采样的趋势图，左侧的图是诸参数 α, βi 和 σ 的边缘后验分布的直方图，并使用核密度估计进行了平滑处理，右侧是针对通过 MCMC 链采样出来的值按顺序绘制的趋势图，其中，对于参数 beta 会有两条边缘后验分布直方图和 MCMC 采样趋势图：

pm.traceplot(trace1)

pm.traceplot(trace2)

（4）统计总结
使用 summary 函数获得诸参数 α, βi 和 σ 的统计总结，我们可以看到 mean 那一列，这一列就是我们得到的最终结果了，大家可以发现它们和我们开始设定的值相当接近了，虽然我们使用了两个不同的采样算法，但现在看来，它们都能对参数们的后验分布进行非常好的采样：

pm.summary(trace1)
1

pm.summary(trace2)
1

好了本文到此就快结束了，以上就是我分享的关于 PyMC3 概率编程入门的内容，大家参照一下即可，重要的还是经过自己的实践来理解 PyMC3 这个优秀的第三方库。