14 贝叶斯回归之先验选择与后验分析

在上一篇中，我们探讨了贝叶斯回归的线性回归模型，了解了线性回归在贝叶斯框架下的理论基础与应用。在本篇中，我们将聚焦于贝叶斯回归中的先验选择与后验分析，帮助您更好地理解模型在不同先验下的行为，以及如何进行后验推断。

先验选择

在贝叶斯统计中，选择合适的先验分布是至关重要的一步。先验反映了我们在观察数据之前对模型参数的不确定性或假设。我们将通过两个基本的先验进行比较：非信息性先验和信息性先验。

非信息性先验

非信息性先验，或称为“扩散先验”，用于表示对参数的最大无偏见态度。例如，对于线性回归模型的参数，我们可以选择以下形式的先验：

$$
\beta \sim \mathcal{N}(0, \tau^2)
$$

其中，$\tau$表示较大的方差，即使得先验对后续数据的影响较小，常见于理论上对参数的“无信息”假设。

信息性先验

信息性先验用于表示对某些参数的先验知识。比如，如果我们认为某个变量的回归系数可能在某个特定值附近，我们可以选择较小的方差，例如：

$$
\beta \sim \mathcal{N}(\mu, \sigma^2)
$$

在这里，$\mu$是我们的信念，而$\sigma$是对此信念的确定性程度。

案例分析

假设我们有一组经济数据，试图通过线性回归来估计某个经济指标（如消费支出）与多个自变量（如收入、教育程度等）之间的关系。使用非信息性先验对参数进行贝叶斯回归，可能会产生宽广的后验分布，鼎盛的展示对数据的依赖性。

而采用信息性先验，假如我们对收入影响消费支出的系数形成了先前的看法（例如，相信这一系数应该超过0.5），后验将更具方向性，快速收敛至我们认为合适的范围。

下面是一个简单的Python示例，展示如何使用PyMC3进行贝叶斯回归并选择不同的先验：

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

# 生成示例数据
np.random.seed(42)
X = np.random.rand(100, 1) * 10
y = 2.5 * X.squeeze() + np.random.randn(100) * 2

# 使用非信息性先验进行贝叶斯回归
with pm.Model() as model_noninformative:
    # 先验设置
    alpha = pm.Normal('alpha', mu=0, sigma=10)
    beta = pm.Normal('beta', mu=0, sigma=10)
    sigma = pm.HalfNormal('sigma', sigma=1)

    # 线性模型
    y_obs = pm.Normal('y_obs', mu=alpha + beta * X.flatten(), sigma=sigma, observed=y)

    # 后验抽样
    trace_noninformative = pm.sample(2000, tune=1000)

# 使用信息性先验进行贝叶斯回归
with pm.Model() as model_informative:
    # 先验设置
    alpha = pm.Normal('alpha', mu=0, sigma=1)
    beta = pm.Normal('beta', mu=0.5, sigma=0.1) # 信息性先验
    sigma = pm.HalfNormal('sigma', sigma=1)

    # 线性模型
    y_obs = pm.Normal('y_obs', mu=alpha + beta * X.flatten(), sigma=sigma, observed=y)

    # 后验抽样
    trace_informative = pm.sample(2000, tune=1000)

# 绘制后验分布比较
pm.plot_posterior(trace_noninformative, fig=plt.figure(figsize=(12, 6)), color='blue')
pm.plot_posterior(trace_informative, fig=myfig, color='orange')
plt.title('Posterior Distributions Comparison')
plt.show()

在上面的代码中，我们分别使用非信息性和信息性先验对相同数据进行贝叶斯回归。通过结果的后验分布的可视化比较，我们可以观察到两种先验对参数估计的影响。

后验分析

在获取了先验和数据的结合后，我们通过贝叶斯定理得到后验分布：

$$
p(\theta | D) = \frac{p(D | \theta) p(\theta)}{p(D)}
$$

其中，$p(D | \theta)$是似然函数，$p(\theta)$是先验分布，而$p(D)$是证据（常数）。

后验推断

一旦得到了后验分布，我们可以通过后验分析来进行推断。重要的推断包括：

1. 点估计：常见的是后验均值、众数或中位数。
2. 可信区间：例如计算95%可信区间以量化不确定性。
3. 模型评估：通过计算后验预测分布与真实数据的拟合程度进行模型比较。

示例分析

假设我们想了解回归系数$\beta$的可信区间。我们可以用以下代码计算并绘制95%可信区间：

import arviz as az

# 提取回归系数的后验样本
beta_samples = trace_noninformative['beta']
pm.plot_posterior(trace_noninformative, var_names=['beta'])

# 计算95%可信区间
ci = np.percentile(beta_samples, [2.5, 97.5])
print(f'95% 可信区间: {ci}')

上述代码使用ArviZ进行后验分析，并绘制了回归系数的后验分布和95%可信区间。

总结

在这篇文章中，我们探讨了贝叶斯回归中先验的选择及其对后验分析的影响。正确地选择先验和进行后验分析是贝叶斯学习中的关键。接下来，我们将在下一篇中进一步讨论贝叶斯回归的预测机制与不确定性量化。这将是您在应用贝叶斯回归时必需掌握的知识。