PyMC3 - 简介和入门马尔可夫蒙特卡洛

PyMC3 - 简介和入门

https://blog.csdn.net/jackxu8/article/details/71080865

线性回归

考虑一个简单的贝叶斯线性回归模型，其参数具有正态分布（Normal）先验。我们预测的具有正态分布的观测值 ，其期望$\mu$ 是两个预测变量$(X_1,X_2)$的线性组合
$$\begin{gathered} YN(\mu ,{\sigma }^2) \\ \mu =\alpha +{\beta }_1{X}_1+{\beta }_2{X}_2 \end{gathered}$$

模拟观察数据

使用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');

<Figure size 640x480 with 0 Axes>

定义模型

PyMC3中的模型定义语言和统计中的公式很接近，一般是一条语句对应一个公式

import pymc3 as pm
from scipy import optimize

basic_model = pm.Model()
with basic_model:
    # 定义三个具有正态分布先验的随机性随机变量（stochastic random veriables）
    alpha=pm.Normal('alpha',mu=0,sd=10)
    beta=pm.Normal('beta',mu=0,sd=10,shape=2)
    sigma=pm.HalfNormal('sigma',sd=1)
    
    # 定义了一个确定性随机变量（deterministic random veriable）
    mu=alpha+beta[0]*X1+beta[1]*X2

    Y_obs=pm.Normal('Y_obs',mu=mu,sd=sigma,observed=Y)
    
    map_estimate = pm.find_MAP(model=basic_model)    
    print(map_estimate)

100.00% [18/18 00:00<00:00 logp = -155.59, ||grad|| = 51.974]

{'alpha': array(0.90661753), 'beta': array([0.94850377, 2.52246124]), 'sigma_log__': array(-0.03771521), 'sigma': array(0.96298715)}

最后定义了 的观测值，这是一个特殊的观测随机变量（observed stochastic），表示模型数据的可能性。通过 observed 参数来告诉这个变量其值是已经被观测到了的（就是 Y 变量），不会被拟合算法改变。可以使用 numpy.ndarry 或 pandas.DataFrame 对象作为数据参数传入。

模型拟合

获得模型中未知变量的后验估计。考虑两个方法:
（1）使用优化方法找到参数的极大后验估计点(maximum a posteriori(MAP))
（2）使用MCMC采样方法获得后验分布来计算。

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

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,return_inferencedata=False)

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'])

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [sigma, beta, alpha]

100.00% [12000/12000 00:01<00:00 Sampling 4 chains, 0 divergences]

Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 2 seconds.


alpha (8000,)
[0.79548864 0.81073319 0.99653017 0.90949895 0.92560485] 

beta (8000, 2)
[[1.02158413 2.29588267]
 [1.02758957 2.28167272]
 [0.96882434 2.83810961]
 ...
 [0.95265641 2.5653976 ]
 [1.09520488 2.51481694]
 [0.7712951  2.54223799]] 

sigma (8000,)
[0.9931085  0.98117391 1.05126755 ... 0.90858098 0.94857458 1.02847602]

可以通过 step 参数指定特定的采样器（迭代器）来替换默认的迭代器NUTS。这里给 sigma 变量指定 Slice 采样器（step），那么其他两个变量（ alpha 和 beta）的采样器会自动指定（NUTS）。

with basic_model:

    # 实例化采样器
    step = pm.Slice(vars=[sigma])

    # 对后验分布进行5000次采样
    trace = pm.sample(5000, step=step,return_inferencedata=False)

Multiprocess sampling (4 chains in 4 jobs)
CompoundStep
>Slice: [sigma]
>NUTS: [beta, alpha]

100.00% [24000/24000 00:04<00:00 Sampling 4 chains, 0 divergences]

Sampling 4 chains for 1_000 tune and 5_000 draw iterations (4_000 + 20_000 draws total) took 5 seconds.

后验分析

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

with basic_model:
    pm.plot_trace(trace);