摘要

概率编程允许在用户自定义的概率模型上进行自动贝叶斯推断。新的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 <script type="math/tex" id="MathJax-Element-1">Y</script> ，其期望 μ<script type="math/tex" id="MathJax-Element-2">\mu</script> 是两个预测变量的线性组合， X1 <script type="math/tex" id="MathJax-Element-3">X_1</script> 和 X2 <script type="math/tex" id="MathJax-Element-4">X_2</script> ：

Yμ∼(μ,σ2)=α+β1X1+β2X2
<script type="math/tex; mode=display" id="MathJax-Element-5">\begin{aligned} Y &\sim \mathcal{N}(\mu, \sigma^2) \\ \mu &= \alpha + \beta_1 X_1 + \beta_2 X_2 \end{aligned}</script>

其中 α <script type="math/tex" id="MathJax-Element-6">\alpha</script> 是截距， βi <script type="math/tex" id="MathJax-Element-7">\beta_i</script> 是变量 Xi <script type="math/tex" id="MathJax-Element-8">X_i</script> 的系数； σ <script type="math/tex" id="MathJax-Element-9">\sigma</script> 代表观察误差。

由于我们需要简历贝叶斯模型，模型中的未知变量需要赋予先验分布，这里选择零均值的正态先验。其中，系数 α <script type="math/tex" id="MathJax-Element-10">\alpha</script> 、 βi <script type="math/tex" id="MathJax-Element-11">\beta_i</script> 的方差为100，代表参数的弱信息。选择半正态分布作为观测误差 σ <script type="math/tex" id="MathJax-Element-12">\sigma</script> 的先验分布。

αβiσ∼(0,100)∼(0,100)∼∣(0,1)∣
<script type="math/tex; mode=display" id="MathJax-Element-13">\begin{aligned} \alpha &\sim \mathcal{N}(0, 100) \\ \beta_i &\sim \mathcal{N}(0, 100) \\ \sigma &\sim \lvert\mathcal{N}(0, 1){\rvert} \end{aligned}</script>

模拟观测数据

使用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)
basic_model = pm.Model()

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

with basic_model:

然后用with语句定义了一个上下文管理器 context manager ，以 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)

上下文中定义了三个具有正态分布先验的随机性随机变量（stochastic random veriables）

mu=alpha+beta[0]*X1+beta[1]*X2

定义了一个确定性随机变量（deterministic random veriable）,确定性指这个变量的值完全由右端值确定（其实就和一般的变量一个意思）

Y_obs=pm.Normal('Y_obs',mu=mu,sd=sigma,observed=Y)

最后定义了 Y <script type="math/tex" id="MathJax-Element-14">Y</script> 的观测值，这是一个特殊的观测随机变量（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)
Optimization terminated successfully.
         Current function value: 149.015731
         Iterations: 13
         Function evaluations: 20
         Gradient evaluations: 20

from scipy import optimize
map_estimate2 = pm.find_MAP(model=basic_model,fmin=optimize.fmin_powell)
Optimization terminated successfully.
         Current function value: 149.015731
         Iterations: 6
         Function evaluations: 269

下面是优化的结果，即点估计结果。α<script type="math/tex" id="MathJax-Element-15">\alpha</script>, βi <script type="math/tex" id="MathJax-Element-16">\beta_i</script>, σ <script type="math/tex" id="MathJax-Element-17">\sigma</script>分别估计出其最优值。

print(map_estimate)
print(map_estimate2)
{'alpha': array(0.9066185315916999), 'beta': array([ 0.94850339,  2.52245771]), 'sigma_log_': array(-0.03278227409845171)}
{'alpha': array(0.9066106657681885), 'beta': array([ 0.94849496,  2.52246567]), 'sigma_log_': array(-0.03278217927386032)}

值得注意的是，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)
Auto-assigning NUTS sampler...
Initializing NUTS using advi...
Average ELBO = -155.08: 100%|██████████| 200000/200000 [00:19<00:00, 10128.66it/s]
Finished [100%]: Average ELBO = -155.12
100%|██████████| 2000/2000 [00:03<00:00, 600.86it/s]

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'])
alpha (2000,)
[ 0.89826662  0.89826662  0.87772947  0.82995563  1.00371943] 

beta (2000, 2)
[[ 0.92705487  2.49968019]
 [ 0.92705487  2.49968019]
 [ 0.93297739  2.50139308]
 ..., 
 [ 0.90528159  2.52904582]
 [ 0.99905303  2.51454039]
 [ 0.99905303  2.51454039]] 

sigma (2000,)
[ 1.00889559  1.00889559  0.95644403 ...,  0.9445938   1.00392469
  1.00392469]

可以通过 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)
Assigned NUTS to alpha
Assigned NUTS to beta


Optimization terminated successfully.
         Current function value: 149.015731
         Iterations: 6
         Function evaluations: 269


100%|██████████| 5000/5000 [00:07<00:00, 676.52it/s]

后验分析

PyMC3提供了 traceplot 函数来绘制后验采样的趋势图。

pm.traceplot(trace);

这里写图片描述

左侧是每个随机变量的边际后验的直方图，使用核密度估计进行了平滑处理。右侧是马尔可夫链采样值按顺序绘制。对于向量参数 beta 会有两条后验分布直方图和后验采样值。

同时也提供蚊子形式的后验统计总结，使用summary函数获得。

pm.summary(trace)
alpha:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------

  0.906            0.098            0.001            [0.713, 1.098]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|

  0.713          0.842          0.903          0.973          1.098


beta:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------

  0.948            0.088            0.001            [0.768, 1.112]
  2.522            0.100            0.002            [2.334, 2.720]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|

  0.771          0.889          0.948          1.008          1.118
  2.329          2.454          2.522          2.590          2.716


sigma:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------

  0.991            0.071            0.001            [0.855, 1.129]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|

  0.862          0.940          0.987          1.037          1.143

前面设置的参数如下

alpha=1
beta =[1, 2.5]
sigma=1

这里通过后验估计可以看出，

mean(alpha)=[0.906]
mean(beta)=[0.948,2.522]
mean(sigma)=[0.991]

得到了准确的估计。

案例学习1：随机波动

这个案例研究一种随机波动现象，随机股市波动，来说明PyMC3在现实问题上的应用。市场波动性是高度非正态的，这导致了对其波动特征的采样异常困难。这个模型有400多个参数，使用传统的采样方法（如Metropolis-Hastins）会很困难，得到高度自相关的结果。这里使用更为高校的NUTS算法。

模型

价格是随时间波动的（日回报），在某些时刻是快速变化的，某些时刻是静止的。随机波动模型对这些随时间变化的变量进行建模。下面的模型参考了NUTS文章（Hoffman 2014, p. 21）：

σνsilog(yi)∼exp(50)∼exp(.1)∼(si−1,σ−2)∼t(ν,0,exp(−2si))
<script type="math/tex; mode=display" id="MathJax-Element-18">\begin{aligned} \sigma &\sim exp(50) \\ \nu &\sim exp(.1) \\ s_i &\sim \mathcal{N}(s_{i-1}, \sigma^{-2}) \\ log(y_i) &\sim t(\nu, 0, exp(-2 s_i)) \end{aligned}</script>

其中 y <script type="math/tex" id="MathJax-Element-19">y</script> 是一个遵循Student-T分布的日回报序列，有一个未知的自由度参数 ν<script type="math/tex" id="MathJax-Element-20">\nu</script>，和一个由潜在过程 s <script type="math/tex" id="MathJax-Element-21">s</script> 决定的缩放因子。在潜在的对数波动过程中，独立的 si<script type="math/tex" id="MathJax-Element-22">s_i</script> 是独立的日对数波动。

数据

数据包含2008年金融危机期间的 S&P 500 的日回报数据。使用 pandas-datareader 从雅虎金融上获取。可以使用pip install pandas-datareader安装。

from pandas_datareader import data
returns=data.get_data_yahoo('SPY',start='2008-5-1',end='2009-12-1')['Adj Close'].pct_change()
print(len(returns))
401

returns.plot(figsize=(10, 6));
plt.ylabel('daily returns in %');

这里写图片描述

模型定义

PyMC3中模型定义直接对应了统计学上的公式定义。这个模型使用了一些新的分布： ν <script type="math/tex" id="MathJax-Element-23">\nu</script> 和 σ <script type="math/tex" id="MathJax-Element-24">\sigma</script> 先验使用了指数分布 Exponential；日回报的使用了Student-T分布 StudentT；隐含波动使用来的高斯随机走动 GaussianRandomWalk。

在PyMC3中，只有正数的先验分布（如 Exponential）使用了对数转换以使采样更加鲁棒，在处理过程中，一个经过对数变换的无约束空间变量（更名为 varName_log）被替代以用于采样。本例中，自由度参数 ν <script type="math/tex" id="MathJax-Element-25">\nu</script> 和隐含过程的缩放银子 σ <script type="math/tex" id="MathJax-Element-26">\sigma</script> 的先验都是指数分布，使用了这种内建变换。具有约束于两边的先验，如beta和Uniform分布，使用了log odds变换进行无约束处理。

不同于PyMC2，PyMC3一般不需要指定迭代初始值，不过为了避免一个无效值出现的情况下，也可以使用testval参数为任何分布提测试，覆盖默认的测试值（通常是均值、中位数）。这个测试值在默认情况下就作为采样和优化操作的起点（初始值）。

潜在波动的先验分布通过 GaussianRandomWalk 给出， GaussianRandomWalk 的值是一个向量，由随机的正态分布走动n个长度获得，sigma确定了正态分布的精度

with pm.Model() as sp500_model:
    sigma = pm.Exponential('sigma',50.0,testval=0.1)
    nu=pm.Exponential('nu',0.1,testval=5.0)
    s=pm.distributions.timeseries.GaussianRandomWalk('s',sigma**-2,shape=len(returns))
    volatility_process=pm.Deterministic('volatility_process',pm.math.exp(-2*s))
    r=pm.StudentT('r',nu,0,1.0/volatility_process,observed=returns)
??pm.StudentT

这里将对数隐含过程 s <script type="math/tex" id="MathJax-Element-27">s</script> 用 exp(-2*s) 转换成了隐含过程。 exp 是Theano的函数，不是NumPy中的函数，Theano提供了NumPy中绝大多数的数学运算函数。

另外，这里的第一句

with pm.Model() as sp500_model:

等价于下面两句

sp500_model = pm.Model()
with sp500_model:

拟合

with sp500_model:
    trace = pm.sample(2000)
Auto-assigning NUTS sampler...
Initializing NUTS using advi...
Average ELBO = 883.23: 100%|██████████| 200000/200000 [01:12<00:00, 2775.64it/s]
Finished [100%]: Average ELBO = 883.14
100%|██████████| 2000/2000 [03:51<00:00, 16.00it/s]

使用 traceplot 检查参数 nu 和 sigma 的采样值。

pm.traceplot(trace);

这里写图片描述

pm.traceplot(trace,[nu,sigma]);

这里写图片描述

pm.traceplot(trace[0:200],[nu,sigma]);

这里写图片描述

将波动路径的分布和采样的波动路径画在同一个图上，采样点的绘制采用半透明方式，这样叠加点越多的地方颜色就越深。

fig,ax=plt.subplots(figsize=(15,8));
returns.plot(ax=ax);
ax.plot(returns.index, 1/np.exp(trace['s',::5].T), 'r', alpha=.03);
ax.set(title='volatility_process', xlabel='time', ylabel='volatility');
ax.legend(['S&P500', 'stochastic volatility process']);

这里写图片描述

案例学习2：煤矿灾难

下面是1851至1962年期间的英国煤矿灾难记录，灾难的数量随着安全防护的进步在逐年减少。但是有一些年份的数据缺失了。

这里建立模型估计灾难数量发生变化的时刻，同时使用多重采样器从随机变量中采样来处理缺失的数据。

丢失数据在NumPy的MaskedArray数据结构中使用-999来标记。

disaster_data = np.ma.masked_values([4, 5, 4, 0, 1, 4, 3, 4, 0, 6, 3, 3, 4, 0, 2, 6,
                            3, 3, 5, 4, 5, 3, 1, 4, 4, 1, 5, 5, 3, 4, 2, 5,
                            2, 2, 3, 4, 2, 1, 3, -999, 2, 1, 1, 1, 1, 3, 0, 0,
                            1, 0, 1, 1, 0, 0, 3, 1, 0, 3, 2, 2, 0, 1, 1, 1,
                            0, 1, 0, 1, 0, 0, 0, 2, 1, 0, 0, 0, 1, 1, 0, 2,
                            3, 3, 1, -999, 2, 1, 1, 1, 1, 2, 4, 2, 0, 0, 1, 4,
                            0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1], value=-999)
year = np.arange(1851, 1962)

plt.plot(year, disaster_data, 'o', markersize=8);
plt.ylabel("Disaster count");
plt.xlabel("Year");

这里写图片描述

时间序列中灾难时间的发生符合泊松分布，分布参数和这个时间序列早起的分布参数有很大的概率是相似的。这个感兴趣的是时间序列改变的位置，这个位置可能和采矿安全规定的变化位置是相关的。

模型如下：

Dtsel∼Pois(rt),rt={l,e,if t<sif t≥s∼Unif(tl,th)∼exp(1)∼exp(1)
<script type="math/tex; mode=display" id="MathJax-Element-28"> \begin{aligned} D_t &\sim \text{Pois}(r_t), r_t= \begin{cases} l, & \text{if } t \lt s \\ e, & \text{if } t \ge s \end{cases} \\ s &\sim \text{Unif}(t_l, t_h)\\ e &\sim \text{exp}(1)\\ l &\sim \text{exp}(1) \end{aligned} </script>

其中
* Dt <script type="math/tex" id="MathJax-Element-29">D_t</script>: 第 t <script type="math/tex" id="MathJax-Element-30">t</script> 年的灾难数量。
* rt<script type="math/tex" id="MathJax-Element-31">r_t</script>: 第 t <script type="math/tex" id="MathJax-Element-32">t</script> 年的泊松分布参数。
* s<script type="math/tex" id="MathJax-Element-33">s</script>: 泊松分布参数改变的那一年。
* e <script type="math/tex" id="MathJax-Element-34">e</script>: 在变化点 s<script type="math/tex" id="MathJax-Element-35">s</script> 之前的泊松分布参数。
* l <script type="math/tex" id="MathJax-Element-36">l</script>: 在变化点 s<script type="math/tex" id="MathJax-Element-37">s</script> 之后的泊松分布参数。
* tl <script type="math/tex" id="MathJax-Element-38">t_l</script>, th <script type="math/tex" id="MathJax-Element-39">t_h</script>: 年份 t <script type="math/tex" id="MathJax-Element-40">t</script> 的上下界。

with pm.Model() as disaster_model:
    switchpoint = pm.DiscreteUniform('switchpoint', lower=year.min(), upper=year.max(), testval=1900)

    # 变化点前后的参数
    early_rate = pm.Exponential('early_rate', 1)
    late_rate = pm.Exponential('late_rate', 1)

    # 定义泊松分布的参数
    rate = pm.math.switch(switchpoint >= year, early_rate, late_rate)

    # 定义分布
    disasters = pm.Poisson('disasters', rate, observed=disaster_data)

当使用 MaskedArray 或 pandas.DataFrame 传递数据是，观测数据 disaster_data 中的缺失数据被标记成NaN，因此将被透明的处理。其背后的处理是创建一个新的随机变量 diasters.missing_values，并从这个随机变量中采样。由于这是个离散随机变量，需要使用Metropolis进行采样。

with disaster_model:
    trace=pm.sample(10000)
Assigned Metropolis to switchpoint
Assigned NUTS to early_rate_log_
Assigned NUTS to late_rate_log_
Assigned Metropolis to disasters_missing
100%|██████████| 10000/10000 [00:12<00:00, 828.95it/s]

从下图(switchpoint)可以看出，对可能的改变时间点的估计大约有10年的跨度。主要概率集中在5年内（1888~1892），switchpoint分布的台阶式跳跃是由于switchpoint和似然函数的关系而不是因为采样误差。

pm.traceplot(trace);

这里写图片描述

via

https://pymc-devs.github.io/pymc3/getting_started.html

阅读全文
AI总结
# python # 开源
Logo
GitCode 开源社区

旨在为数千万中国开发者提供一个无缝且高效的云端环境，以支持学习、使用和贡献开源项目。

更多推荐

cover

“盲水印”技术新标杆：blind_watermark加入GitCode

avatar GitCode 开源社区
cover

颠覆开发效率！国内首个开源微服务编排框架Juggle入驻GitCode

avatar GitCode 开源社区
cover

开源平台GitCode招人啦，大量好岗！

avatar GitCode 开源社区