PyMC3 Models 入门指南：从线性回归开始

2025-07-07 08:09:44作者：齐冠琰

前言

PyMC3 Models 是一个基于 PyMC3 的高级建模库，它提供了类似于 scikit-learn 的简洁 API，让贝叶斯建模变得更加容易。本文将带你从零开始，通过一个简单的线性回归示例，了解如何使用 PyMC3 Models 进行贝叶斯建模。

准备工作

在开始之前，我们需要准备一些模拟数据来演示建模过程：

import numpy as np

# 生成1000个样本的模拟数据
X = np.random.randn(1000, 1)  # 特征数据
noise = 2 * np.random.randn(1000, 1)  # 噪声
Y = 4 * X + 3 + noise  # 目标变量，真实斜率为4，截距为3

这段代码生成了一个简单的线性关系数据集，其中包含高斯噪声。我们可以通过可视化来观察数据的分布情况。

传统PyMC3建模方式

为了更好地理解 PyMC3 Models 的价值，我们先看看如何使用原生 PyMC3 进行建模。

1. 模型定义

import pymc3 as pm
import theano

lin_reg_model = pm.Model()

# 使用theano共享变量以便后续使用minibatch
model_input = theano.shared(X)
model_output = theano.shared(Y)

with lin_reg_model:
    # 定义先验分布
    alpha = pm.Normal('alpha', mu=0, sd=100, shape=(1))  # 截距
    beta = pm.Normal('beta', mu=0, sd=100, shape=(1))  # 斜率
    
    # 定义噪声的标准差
    s = pm.HalfNormal('s', tau=1)
    
    # 定义线性模型
    mean = alpha + beta * model_input
    
    # 定义似然函数
    y = pm.Normal('y', mu=mean, sd=s, observed=model_output)

2. 参数推断

PyMC3 提供了多种推断方法，这里我们使用变分推断(ADVI)：

with lin_reg_model:
    inference = pm.ADVI()
    approx = pm.fit(
        n=20000, 
        method=inference, 
        more_replacements={
            model_input: pm.Minibatch(X), 
            model_output: pm.Minibatch(Y)
        }
    )

我们可以通过绘制ELBO(证据下界)来检查收敛情况：

import matplotlib.pyplot as plt

plt.plot(-inference.hist)
plt.ylabel('ELBO')
plt.xlabel('iteration')
plt.show()

3. 结果解释

从近似后验分布中采样并查看结果：

trace = approx.sample(draws=5000)
pm.summary(trace)
pm.traceplot(trace)

4. 预测新数据

使用后验预测检查(PPC)进行预测：

from sklearn.metrics import r2_score

ppc = pm.sample_ppc(trace[1000:], model=lin_reg_model, samples=2000)
pred = ppc['y'].mean(axis=0)
print(r2_score(Y, pred))  # 输出R²分数

使用PyMC3 Models简化建模

现在，让我们看看如何使用 PyMC3 Models 来完成同样的任务，代码将更加简洁。

1. 模型拟合

from pymc3_models import LinearRegression

LR = LinearRegression()
LR.fit(X, Y, minibatch_size=100)
LR.plot_elbo()  # 可视化训练过程

2. 结果分析

PyMC3 Models 自动保存了trace，我们可以直接使用PyMC3的函数进行分析：

pm.traceplot(LR.trace)

3. 预测与评估

Y_predict = LR.predict(X)
print(LR.score(X, Y))  # 计算R²分数

如果需要获取点估计值：

beta = LR.summary['mean']['betas__0_0']
alpha = LR.summary['mean']['alpha__0']

高级功能

模型保存与加载

PyMC3 Models 提供了方便的模型保存和加载功能：

LR.save('model_save_path/')  # 保存模型

# 加载模型
new_LR = LinearRegression()
new_LR.load('model_save_path/')

使用NUTS采样器

虽然默认使用ADVI变分推断，但也可以选择NUTS采样器：

LR_nuts = LinearRegression()
LR_nuts.fit(X, Y, inference_type='nuts', inference_args={'draws': 2000})

自定义推断参数

可以通过inference_args参数自定义推断过程的参数：

LR_custom = LinearRegression()
LR_custom.fit(X, Y, inference_args={'n': 30000, 'method': pm.ADVI()})

构建自定义模型

PyMC3 Models 提供了BayesianModel基类，可以基于它构建自己的模型：

from pymc3_models import BayesianModel

class CustomModel(BayesianModel):
    def __init__(self):
        super(CustomModel, self).__init__()
        
    def create_model(self):
        # 在这里定义你的PyMC3模型
        pass
        
    def fit(self, X, y, **kwargs):
        # 自定义拟合过程
        pass
        
    def predict(self, X):
        # 自定义预测方法
        pass