概率编程基础实例——线性回归¶
概率编程允许在用户自定义的概率模型上进行自动贝叶斯推断。新的MCMC(Markoc chain Monte Carlo)采样方法允许在复杂模型上进行推断。这类MCMC采样方法被称为HMC(Hamliltinian Monte Carlo)，但是其推断需要的梯度信息有时候是不获得的。PyMC3是一个用Python编写的开源的概率编程框架，使用Theano通过变分推理进行梯度计算，并使用了C实现加速运算。不同于其他概率编程语言，PyMC3允许使用Python代码来定义模型。这种没有作用域限制的语言极大的方便了模型定义和直接交互。这篇文章介绍了这个软件包。
PyMC3具有先进的下一代MCMC采样算法如No-U-Turn Sampler (NUTS; Hoffman, 2014)和Hamiltonian Monte Carlo自整定变体(HMC; Duane, 1987)。这类采样算法在高维和复杂的后验分布上具有良好的效果，允许对复杂模型进行拟合而不需要对拟合算法有特殊的了解。NUTS和HMC算法从似然函数中获得梯度信息，因此其收敛速度比传统采样方法快很多，特别是针对大模型。NUTS也具有集合自整定过程，因此使用者不需要了解算法细节。
考虑一个简单的贝叶斯线性回归模型，其参数具有正态分布 (Normal) 先验。我们预测的具有正态分 布的观测值 $Y$ ，其期望 $\mu$ 是两个预测变量的线性组合， $X_1$ 和 $X_2$ :
$$
\begin{aligned}
&Y \sim \mathcal{N}\left(\mu, \sigma^2\right) \\
&\mu=\alpha+\beta_1 X_1+\beta_2 X_2
\end{aligned}
$$
其中 $\alpha$ 是截距， $\beta_i$ 是变量 $X_i$ 的系数； $\sigma$ 代表观察误差。
由于我们需要简历贝叶斯模型，模型中的末知变量需要赋予先验分布，这里选择零均值的正态先验。
其中，系数 $\alpha 、 \beta_i$ 的方差为 100 ，代表参数的弱信息。选择半正态分布作为观测误差 $\sigma$ 的先验分 布。
$$
\begin{aligned}
\alpha & \sim \mathcal{N}(0,100) \\
\beta_i & \sim \mathcal{N}(0,100) \\
\sigma & \sim|\mathcal{N}(0,1)|
\end{aligned}
$$
模拟观测数据¶
使用NumPy的随机函数 random 模块来产生模拟数据，再使用PyMC3尝试恢复相应的参数。
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
np.random.seed(123)
alpha=1
sigma=1
beta =[1, 2.5]
N=100
X1=np.random.randn(N)
X2=np.random.randn(N)
Y=alpha + beta[0]*X1 + beta[1]*X2 + np.random.randn(N)*sigma
%matplotlib inline
fig1,ax1 = plt.subplots(1, 2, figsize=(10,4));
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');
fig2 = plt.figure(2);
ax2 = Axes3D(fig2);
ax2.scatter(X1,X2,Y);
ax2.set_xlabel('X1');
ax2.set_ylabel('X2');
ax2.set_zlabel('Y');
定义模型¶
PyMC3中的模型定义语言和统计中的公式很接近，一般是一条语句对应一个公式
import pymc3 as pm
basic_model = pm.Model()
with basic_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_obs=pm.Normal('Y_obs',mu=mu,sd=sigma,observed=Y)
第一行定义了一个新的模型对象，这个对象是模型中随机变量的容器,然后用with语句定义了一个上下文管理器 context manager ，以 basic_model 作为上下文，在这个上下文中定义的变量都被添加到这个模型。上下文中定义了三个具有正态分布先验的随机性随机变量（stochastic random veriables）定义了一个确定性随机变量（deterministic random veriable）,确定性指这个变量的值完全由右端值确定（其实就和一般的变量一个意思）,最后定义了 Y 的观测值，这是一个特殊的观测随机变量（observed stochastic），表示模型数据的可能性。通过 observed 参数来告诉这个变量其值是已经被观测到了的（就是 Y 变量），不会被拟合算法改变。可以使用 numpy.ndarry 或 pandas.DataFrame 对象作为数据参数传入。
模型拟合¶
获得模型中未知变量的后验估计。考虑两个方法：
- （1）使用优化方法找到参数的极大后验估计点(maximum a posteriori(MAP))
- （2）使用MCMC采样方法获得后验分布来计算。
极大后验估计 MAP¶
使用最优化方法找到点估计，快速简单。PyMC3提供了 find_MAP 函数，返回参数的一个估计值（point）。
默认情况下，find_MAP 使用Broyden–Fletcher–Goldfarb–Shanno (BFGS) 算法进行最优化运算，找到对数后验分布的最大值。这里也可以指定 scipy.optimize 模块中的其他最优化算法完成寻优。
map_estimate = pm.find_MAP(model=basic_model)
值得注意的是，MAP估计不总是有效的，特别是在极端模型情况下。在高维后验中，可能出现某一出概率密度异常大，但整体概率很小的情况，在层次化模型中较为常见。
大多数寻找MAP极大值算法都找到局部最优解，那么这种算法这在差异非常大的多模态后验中会失效。
采样算法¶
虽然极大后验估计是一个简单快速的估计未知参数的办法，但是没有对不确定性进行估计是其缺陷。相反，比如MCMC这类基于采样的方法可以获得后验分布接近真实的采样值。
为了使用MCMC采样以获得后验分布的样本，PyMC3需要指定和特定MCMC算法相关的步进方法（采样方法）,如Metropolis, Slice sampling, or the No-U-Turn Sampler (NUTS)。PyMC3中的 step_methods 子模块提供了如下采样器： NUTS, Metropolis, Slice, HamiltonianMC, and BinaryMetropolis。可以由PyMC3自动指定也可以手动指定。自动指定是根据模型中的变量类型决定的：
- 二值变量：指定 BinaryMetropolis
- 离散变量：指定 Metropolis
- 连续变量：指定 NUTS
手动指定优先，可覆盖自动指定。
基于梯度的采样算法¶
PyMC3有许多基本采样算法，如自适应切片采样、自适应Metropolis-Hastings采样，但最厉害是的No-U-Turn采样算法（NUTS），特别是在具有大量连续型参数的模型中。NUTS基于对数概率密度的梯度，利用了高概率密度的区域信息，以获得比传统采样算法更快的收敛速度。PyMC3使用Theano的自动微分运算来进行后验分布的梯度计算。通过自整定步骤来自适应的设置Hamiltonian Monte Carlo(HMC)算法的可调参数。NUTS不可用于不可微分变量（离散变量）的采样，但是对于具有不可微分变量的模型中的可微分变量是可以使用的。
NUTS需要一个缩放矩阵作为参数，这个矩阵提供了分布的大致形状，使NUTS不至于在某些方向上的步长过大而在另外一些方向上的步长过小。设置有效的缩放矩阵有助于获得高效率的采样，特别是对于那些有很多未知的（未被观察的）随机型随机变量的模型或者具有高度非正态后验的模型。不好的缩放矩阵将导致采样速度很慢甚至完全停止。此外，在多数时候，采样的初始位置对于高效的采样也是很关键的。
幸运的是，PyMC3使用自动变分推理ADVI (auto-diff variational inference)来初始化NUTS算法，并在 step 参数没有被指定的情况下会自动指定一个合适的迭代方法（step，采样器）。
下面的例子对 basic_model 的后验分布进行了2000次采样。
with basic_model:
trace = pm.sample(2000)
sample 函数用指定的迭代器（采样算法）进行了2000次迭代，收集到的采样值存储在 Trace 对象中，按照采样值获得的先后顺序存储。Trace对象是一个字典中包含了多个map。
print("alpha",trace['alpha'].shape)
print(trace['alpha'][0:5],"\n")
print("beta",trace['beta'].shape)
print(trace['beta'],"\n")
print("sigma",trace['sigma'].shape)
print(trace['sigma'])
可以通过 step 参数指定特定的采样器（迭代器）来替换默认的迭代器NUTS。这里给 sigma 变量指定 Slice 采样器（step），那么其他两个变量（ alpha 和 beta）的采样器会自动指定（NUTS）。
with basic_model:
# 用MAP获得初始点
start = pm.find_MAP(fmin=optimize.fmin_powell)
# 实例化采样器
step = pm.Slice(vars=[sigma])
# 对后验分布进行5000次采样
trace = pm.sample(5000, step=step,start=start)
后验分析¶
PyMC3提供了 traceplot 函数来绘制后验采样的趋势图。
左侧是每个随机变量的边际后验的直方图，使用核密度估计进行了平滑处理。右侧是马尔可夫链采样值按顺序绘制。对于向量参数 beta 会有两条后验分布直方图和后验采样值。
同时也提供文字形式的后验统计总结，使用summary函数获得。
参考资料
