`Python`的`概率`编程：`PyMC3`和`Stan`在`贝叶斯`建模中的`应用`。

Python 概率编程：PyMC3 和 Stan 在贝叶斯建模中的应用

大家好，今天我们来探讨 Python 中的概率编程，特别是聚焦于两个强大的库：PyMC3 和 Stan，以及它们在贝叶斯建模中的应用。贝叶斯建模提供了一种量化不确定性、整合先验知识并进行预测的强大框架。PyMC3 和 Stan 提供了灵活的工具，使我们能够构建、拟合和分析复杂的贝叶斯模型。

1. 贝叶斯建模基础

在深入讨论 PyMC3 和 Stan 之前，我们先快速回顾一下贝叶斯建模的核心概念。贝叶斯定理是贝叶斯统计的基石：

P(θ|D) = [P(D|θ) * P(θ)] / P(D)

其中：

• P(θ|D)：后验概率（Posterior probability），给定数据 D 后，参数 θ 的概率。这是我们最感兴趣的部分，它反映了在观察到数据后我们对参数的信念。
• P(D|θ)：似然函数（Likelihood function），给定参数 θ 时，观察到数据 D 的概率。它衡量了模型对数据的拟合程度。
• P(θ)：先验概率（Prior probability），在观察到任何数据之前，参数 θ 的概率。它代表了我们对参数的初始信念。
• P(D)：证据（Evidence）或归一化常数，观察到数据 D 的概率。在实际计算中，通常不需要显式计算它，因为后验概率可以通过比例关系确定：P(θ|D) ∝ P(D|θ) * P(θ)。

贝叶斯建模的核心思想是利用先验知识和观测数据来更新我们对参数的信念，最终得到后验分布。

2. PyMC3 简介

PyMC3 是一个基于 Theano 的 Python 概率编程库。它允许我们使用 Python 语法定义贝叶斯模型，并使用马尔可夫链蒙特卡洛 (MCMC) 等方法进行后验推断。

2.1 PyMC3 的主要特点：

• 模型定义简易性： 使用 Python 语法定义模型，易于理解和修改。
• 多种推断算法： 支持多种 MCMC 采样器，如 Metropolis-Hastings、NUTS (No-U-Turn Sampler) 等。
• 扩展性： 基于 Theano，可以利用 GPU 加速计算。
• 可视化工具： 集成了多种可视化工具，方便分析后验分布。

2.2 PyMC3 模型构建示例：线性回归

让我们通过一个简单的线性回归模型来演示 PyMC3 的使用。假设我们有以下数据：

import numpy as np
import pymc3 as pm
import matplotlib.pyplot as plt
import seaborn as sns

# 生成模拟数据
np.random.seed(123)
size = 100
true_intercept = 1
true_slope = 2
x = np.linspace(0, 1, size)
y = true_intercept + true_slope * x + np.random.normal(loc=0, scale=0.5, size=size)

data = dict(x=x, y=y)

# 可视化数据
plt.figure(figsize=(7, 5))
plt.scatter(x, y, label="Data")
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.show()

现在，我们使用 PyMC3 构建线性回归模型：

with pm.Model() as model:
    # 定义先验分布
    alpha = pm.Normal("alpha", mu=0, sigma=10) #截距
    beta = pm.Normal("beta", mu=0, sigma=10) #斜率
    sigma = pm.HalfNormal("sigma", sigma=1) #标准差

    # 定义似然函数
    mu = alpha + beta * data["x"]
    likelihood = pm.Normal("y", mu=mu, sigma=sigma, observed=data["y"])

    # 进行后验推断
    trace = pm.sample(1000, tune=1000, cores=4)

代码解释：

• pm.Model()：创建一个 PyMC3 模型上下文。
• pm.Normal("alpha", mu=0, sigma=10)：定义 alpha (截距) 的先验分布为均值为 0，标准差为 10 的正态分布。
• pm.Normal("beta", mu=0, sigma=10)：定义 beta (斜率) 的先验分布为均值为 0，标准差为 10 的正态分布。
• pm.HalfNormal("sigma", sigma=1)：定义 sigma (标准差) 的先验分布为标准差为 1 的半正态分布（只取正值）。
• mu = alpha + beta * data["x"]：定义线性回归模型的均值。
• pm.Normal("y", mu=mu, sigma=sigma, observed=data["y"])：定义似然函数，观测数据 y 服从均值为 mu，标准差为 sigma 的正态分布。
• pm.sample(1000, tune=1000, cores=4)：使用 MCMC 采样器进行后验推断，采样 1000 个样本，tune 1000 个样本用于调整采样器参数，并使用 4 个 CPU 核心并行计算。

2.3 PyMC3 模型分析：

采样完成后，我们可以使用 PyMC3 的可视化工具分析后验分布：

pm.traceplot(trace, var_names=["alpha", "beta", "sigma"])
plt.show()
pm.summary(trace, var_names=["alpha", "beta", "sigma"])

pm.traceplot() 可以绘制每个参数的采样轨迹和后验分布直方图，用于检查 MCMC 采样的收敛性。pm.summary() 可以提供后验分布的统计摘要，如均值、标准差、HPD (Highest Posterior Density) 区间等。

2.4 PyMC3 预测：

获得后验分布后，我们可以使用它进行预测：

# 生成新的 x 值用于预测
x_new = np.linspace(0, 1, 100)

# 使用后验样本进行预测
with model:
    ppc = pm.sample_posterior_predictive(trace, var_names=["y"], samples=500,  progressbar=False)

# 计算预测均值和置信区间
y_pred = ppc["y"].mean(axis=0)
y_lower = np.percentile(ppc["y"], 2.5, axis=0)
y_upper = np.percentile(ppc["y"], 97.5, axis=0)

# 可视化预测结果
plt.figure(figsize=(7, 5))
plt.scatter(x, y, label="Data")
plt.plot(x_new, y_pred, "r-", label="Prediction")
plt.fill_between(x_new, y_lower, y_upper, color="r", alpha=0.2, label="95% CI")
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.show()

pm.sample_posterior_predictive() 使用后验样本生成新的预测数据，然后我们可以计算预测均值和置信区间，并将其可视化。

2.5 复杂模型示例：混合模型

PyMC3 不仅可以处理简单的线性回归，还可以构建更复杂的模型。例如，我们可以构建一个混合模型，假设数据来自多个不同的分布：

# 生成模拟数据
np.random.seed(123)
n_samples = 500
mixture_probs = np.array([0.3, 0.7])
means = np.array([0, 5])
stds = np.array([1, 1])

components = np.random.choice(2, size=n_samples, p=mixture_probs)
data = np.random.normal(means[components], stds[components], size=n_samples)

# 可视化数据
plt.figure(figsize=(7, 5))
sns.histplot(data, kde=True)
plt.xlabel("Data")
plt.ylabel("Frequency")
plt.show()

# 构建混合模型
with pm.Model() as model_mixture:
    # 定义先验分布
    p = pm.Dirichlet("p", a=np.array([1, 1])) #混合比例
    mu = pm.Normal("mu", mu=np.array([0, 5]), sigma=10, shape=2) #每个分布的均值
    sigma = pm.HalfNormal("sigma", sigma=1, shape=2) #每个分布的标准差

    # 定义似然函数
    obs = pm.Mixture("obs", w=p, mu=mu, sigma=sigma, observed=data)

    # 进行后验推断
    trace_mixture = pm.sample(1000, tune=1000, cores=4)

在这个例子中，我们使用 pm.Dirichlet 定义混合比例的先验分布，使用 pm.Normal 和 pm.HalfNormal 定义每个分布的均值和标准差的先验分布，然后使用 pm.Mixture 定义似然函数，将数据建模为来自多个正态分布的混合。

3. Stan 简介

Stan 是一个用于统计建模和高性能计算的概率编程语言。它提供了一种灵活的语法来定义贝叶斯模型，并使用高效的 MCMC 算法（如 NUTS）进行后验推断。与 PyMC3 不同，Stan 本身不是一个 Python 库，而是一种独立的语言和编译器。我们需要使用 Python 接口（如 PyStan 或 CmdStanPy）来与 Stan 交互。

3.1 Stan 的主要特点：

• 高性能： 使用 C++ 编写，编译成高效的机器码，运行速度快。
• 强大的 MCMC 算法： 内置了 NUTS 等先进的 MCMC 算法，能够高效地探索高维后验分布。
• 灵活的建模语言： Stan 的建模语言非常灵活，可以表达复杂的统计模型。
• 自动微分： Stan 能够自动计算模型的梯度，方便进行优化和推断。

3.2 Stan 模型构建示例：线性回归

我们使用与 PyMC3 相同的线性回归模型来演示 Stan 的使用。首先，我们需要编写 Stan 代码来定义模型：

data {
  int<lower=0> N;
  vector[N] x;
  vector[N] y;
}
parameters {
  real alpha;
  real beta;
  real<lower=0> sigma;
}
model {
  y ~ normal(alpha + beta * x, sigma);
  alpha ~ normal(0, 10);
  beta ~ normal(0, 10);
  sigma ~ exponential(1);
}

代码解释：

• data 块：声明输入数据，包括数据点的数量 N，自变量 x 和因变量 y。
• parameters 块：声明需要推断的参数，包括截距 alpha，斜率 beta 和标准差 sigma。
• model 块：定义模型，包括似然函数和先验分布。
  • y ~ normal(alpha + beta * x, sigma)：定义似然函数，观测数据 y 服从均值为 alpha + beta * x，标准差为 sigma 的正态分布。
  • alpha ~ normal(0, 10)：定义 alpha 的先验分布为均值为 0，标准差为 10 的正态分布。
  • beta ~ normal(0, 10)：定义 beta 的先验分布为均值为 0，标准差为 10 的正态分布。
  • sigma ~ exponential(1)：定义 sigma 的先验分布为均值为 1 的指数分布。

然后，我们可以使用 PyStan 或 CmdStanPy 从 Python 中调用 Stan：

import stan
import numpy as np

# Stan 代码
stan_code = """
data {
  int<lower=0> N;
  vector[N] x;
  vector[N] y;
}
parameters {
  real alpha;
  real beta;
  real<lower=0> sigma;
}
model {
  y ~ normal(alpha + beta * x, sigma);
  alpha ~ normal(0, 10);
  beta ~ normal(0, 10);
  sigma ~ exponential(1);
}
"""

# 准备数据
data = {
    "N": len(x),
    "x": x,
    "y": y
}

# 编译 Stan 模型
model = stan.build(stan_code, data=data)

# 进行后验推断
fit = model.sample(num_chains=4, num_samples=1000)

代码解释：

• stan.build(stan_code, data=data)：编译 Stan 模型，将 Stan 代码编译成 C++ 代码，并生成可执行文件。
• fit = model.sample(num_chains=4, num_samples=1000)：使用 MCMC 采样器进行后验推断，采样 1000 个样本，并使用 4 个链并行计算。

3.3 Stan 模型分析：

采样完成后，我们可以使用 PyStan 提供的函数分析后验分布：

# 获取后验样本
posterior = fit.to_frame()

# 打印后验摘要
print(posterior.describe())

# 可视化后验分布
import matplotlib.pyplot as plt
import seaborn as sns

plt.figure(figsize=(12, 4))
sns.kdeplot(posterior["alpha"], label="alpha")
sns.kdeplot(posterior["beta"], label="beta")
sns.kdeplot(posterior["sigma"], label="sigma")
plt.legend()
plt.show()

3.4 Stan 预测：

与 PyMC3 类似，我们可以使用后验分布进行预测。我们需要修改 Stan 代码，添加 generated quantities 块：

data {
  int<lower=0> N;
  vector[N] x;
  vector[N] y;
  int<lower=0> N_new;
  vector[N_new] x_new;
}
parameters {
  real alpha;
  real beta;
  real<lower=0> sigma;
}
model {
  y ~ normal(alpha + beta * x, sigma);
  alpha ~ normal(0, 10);
  beta ~ normal(0, 10);
  sigma ~ exponential(1);
}
generated quantities {
  vector[N_new] y_pred;
  for (i in 1:N_new) {
    y_pred[i] = normal_rng(alpha + beta * x_new[i], sigma);
  }
}

代码解释：

• data 块：添加新的自变量 x_new 和数量 N_new。
• generated quantities 块：生成新的预测数据 y_pred。
  • y_pred[i] = normal_rng(alpha + beta * x_new[i], sigma)：使用后验样本生成新的预测数据，服从均值为 alpha + beta * x_new[i]，标准差为 sigma 的正态分布。

然后，我们可以使用 PyStan 或 CmdStanPy 再次调用 Stan，并获取预测结果：

# 准备新的 x 值用于预测
x_new = np.linspace(0, 1, 100)

# 准备数据
data = {
    "N": len(x),
    "x": x,
    "y": y,
    "N_new": len(x_new),
    "x_new": x_new
}

# 编译 Stan 模型
model = stan.build(stan_code, data=data)

# 进行后验推断
fit = model.sample(num_chains=4, num_samples=1000)

# 获取后验样本
posterior = fit.to_frame()

# 获取预测结果
y_pred = posterior[[col for col in posterior.columns if col.startswith('y_pred')]].values

# 计算预测均值和置信区间
y_pred_mean = y_pred.mean(axis=0)
y_pred_lower = np.percentile(y_pred, 2.5, axis=0)
y_pred_upper = np.percentile(y_pred, 97.5, axis=0)

# 可视化预测结果
plt.figure(figsize=(7, 5))
plt.scatter(x, y, label="Data")
plt.plot(x_new, y_pred_mean, "r-", label="Prediction")
plt.fill_between(x_new, y_pred_lower, y_pred_upper, color="r", alpha=0.2, label="95% CI")
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.show()

4. PyMC3 vs. Stan：如何选择？

PyMC3 和 Stan 都是强大的贝叶斯建模工具，但它们在某些方面有所不同。选择哪个工具取决于具体的需求和偏好。

特性	PyMC3	Stan
语言	Python	Stan (独立的建模语言)
易用性	更容易上手，Python 语法更易理解	学习曲线较陡峭，需要学习 Stan 建模语言
性能	较慢，依赖 Theano 或 Aesara 进行计算	更快，使用 C++ 编写，编译成机器码
MCMC 算法	支持多种 MCMC 算法，但 NUTS 的实现较慢	内置高效的 NUTS 算法，能够处理高维问题
扩展性	基于 Theano/Aesara，可以利用 GPU 加速	扩展性强，可以自定义函数和分布
社区支持	活跃的 Python 社区	活跃的 Stan 社区
模型复杂度	适合中等复杂度的模型	适合复杂、高性能要求的模型

建议：

• PyMC3： 如果你熟悉 Python，需要快速构建和迭代模型，或者模型复杂度不高，PyMC3 是一个不错的选择。
• Stan： 如果你需要处理复杂的模型，对性能有较高要求，或者需要使用高级的 MCMC 算法，Stan 是一个更好的选择。

5. 贝叶斯建模中的常见问题

在贝叶斯建模中，有一些常见的问题需要注意：

• 先验选择： 先验分布的选择会影响后验分布。需要根据领域知识和建模目标选择合适的先验分布。可以使用信息先验（informative prior）来加入先验知识，也可以使用无信息先验（uninformative prior）来尽量减少先验对后验的影响。
• 模型诊断： MCMC 采样的收敛性需要进行诊断，可以使用迹图 (traceplot)、自相关图 (autocorrelation plot) 和 R-hat 统计量等工具进行检查。
• 模型验证： 需要验证模型的预测能力，可以使用后验预测检验 (posterior predictive check) 等方法。
• 计算成本： 复杂的贝叶斯模型可能需要大量的计算资源和时间。需要优化模型和采样算法，或者使用 GPU 加速等方法。

6. 总结与展望

我们探讨了 Python 中的概率编程，特别是 PyMC3 和 Stan 在贝叶斯建模中的应用。PyMC3 提供了简单易用的 Python 接口，适合快速构建和迭代模型。Stan 提供了高性能的计算能力和灵活的建模语言，适合处理复杂的模型。在实际应用中，需要根据具体的需求和偏好选择合适的工具。

随着概率编程的发展，我们可以期待更多强大的工具和算法的出现，这将使得贝叶斯建模更加便捷和高效。未来，我们可以探索更复杂的模型，例如深度贝叶斯模型、非参数贝叶斯模型等，并将其应用于更广泛的领域，例如自然语言处理、计算机视觉和强化学习等。

7. 更多资源

• PyMC3 官方文档：https://www.pymc.io/projects/docs/en/v3/
• Stan 官方网站：https://mc-stan.org/
• PyStan 文档：https://pystan.readthedocs.io/en/latest/
• CmdStanPy 文档：https://cmdstanpy.readthedocs.io/en/latest/

希望今天的讲座对您有所帮助！谢谢大家！

8. 理解建模工具，持续探索

贝叶斯建模是一个强大的工具，PyMC3和Stan各有优势。选择合适的工具，不断学习和实践，将有助于我们更好地理解数据和做出预测。