贝叶斯线性模型实战：从联合高斯到条件高斯的完整推导与应用

贝叶斯线性模型实战：从联合高斯到条件高斯的完整推导与应用

如果你接触过机器学习，大概率听说过贝叶斯方法。它不像频率学派那样把模型参数当作固定值，而是将其视为随机变量，用概率分布来描述我们的“信念”。这种视角的转换，带来的不仅是理论上的优雅，更是在面对数据稀缺、噪声复杂或需要量化预测不确定性时，提供了强大的工具。今天，我们不打算只停留在理论层面，而是想和你一起，亲手推导一遍贝叶斯线性回归的核心——联合高斯与条件高斯分布的关系，并用Python代码将其落地，解决一个真实的回归预测问题。你会发现，那些看似复杂的矩阵运算背后，隐藏着极其直观的统计直觉。

这篇文章面向的是希望深入理解贝叶斯模型底层逻辑，并渴望将其应用于实际项目的工程师和统计学习者。我们将从零开始，构建一个完整的贝叶斯线性回归工作流，涵盖数据生成、先验设定、后验推导与计算，并深入探讨一个常被忽视但对结果影响巨大的细节：协方差矩阵的选择。准备好了吗？让我们开始这场从理论到实战的旅程。

1. 重温高斯分布：从一元到多元的跃迁

在深入贝叶斯线性模型之前，我们必须对高斯（正态）分布有坚实的理解。它不仅是自然界中最常见的分布，更是贝叶斯推断中“共轭先验”这一美妙性质的核心载体。

1.1 多元高斯分布的本质

一元高斯分布由均值和方差决定。到了多元世界，我们面对的是随机向量 $\boldsymbol{x} \in \mathbb{R}^p$。其概率密度函数为：

$$ p(\boldsymbol{x}) = \frac{1}{(2\pi)^{p/2} |\boldsymbol{\Sigma}|^{1/2}} \exp\left( -\frac{1}{2} (\boldsymbol{x} - \boldsymbol{\mu})^T \boldsymbol{\Sigma}^{-1} (\boldsymbol{x} - \boldsymbol{\mu}) \right) $$

这里，$\boldsymbol{\mu}$ 是均值向量，$\boldsymbol{\Sigma}$ 是协方差矩阵。协方差矩阵不仅包含了各个维度的方差（对角线元素），还包含了任意两个维度之间的协方差（非对角线元素），它完整刻画了随机向量各分量之间的线性相关关系。

注意：协方差矩阵 $\boldsymbol{\Sigma}$ 必须是对称且正定的。正定性保证了其逆矩阵存在，并且指数项内的二次型始终为正，这在计算中至关重要。

一个关键且强大的性质是：高斯随机向量的任意线性变换，其结果仍然是高斯分布。假设 $\boldsymbol{x} \sim \mathcal{N}(\boldsymbol{\mu}, \boldsymbol{\Sigma})$，对于矩阵 $\boldsymbol{A} \in \mathbb{R}^{m \times p}$ 和向量 $\boldsymbol{b} \in \mathbb{R}^{m}$，有：

$$ \boldsymbol{y} = \boldsymbol{A}\boldsymbol{x} + \boldsymbol{b} \quad \Rightarrow \quad \boldsymbol{y} \sim \mathcal{N}(\boldsymbol{A}\boldsymbol{\mu} + \boldsymbol{b}, \boldsymbol{A}\boldsymbol{\Sigma}\boldsymbol{A}^T) $$

这个性质是后续所有推导的基石。我们可以通过构造不同的 $\boldsymbol{A}$ 来方便地求取边际分布。例如，想得到 $\boldsymbol{x}$ 的前 $k$ 个分量的边际分布，只需令 $\boldsymbol{A} = [\boldsymbol{I}_k, \boldsymbol{0}]$。

1.2 联合、边际与条件分布：一个三位一体的框架

当我们有两个联合高斯分布的随机向量 $\boldsymbol{x}$ 和 $\boldsymbol{y}$ 时，它们的关系可以通过联合分布完整描述：

$$ \begin{bmatrix} \boldsymbol{x} \ \boldsymbol{y} \end{bmatrix} \sim \mathcal{N}\left( \begin{bmatrix} \boldsymbol{\mu}x \ \boldsymbol{\mu}y \end{bmatrix}, \begin{bmatrix} \boldsymbol{\Sigma}{xx} & \boldsymbol{\Sigma}{xy} \ \boldsymbol{\Sigma}{yx} & \boldsymbol{\Sigma}{yy} \end{bmatrix} \right) $$

• 边际分布：直接对联合分布“积分”掉不关心的变量。利用线性变换性质，$\boldsymbol{x}$ 的边际分布就是 $\mathcal{N}(\boldsymbol{\mu}x, \boldsymbol{\Sigma}{xx})$。这非常直观。
• 条件分布：这是在已知 $\boldsymbol{y}$ 取某个值时，$\boldsymbol{x}$ 的分布。其公式是贝叶斯推断的核心：

$$ \boldsymbol{x} | \boldsymbol{y} \sim \mathcal{N}\left( \boldsymbol{\mu}{x|y},\ \boldsymbol{\Sigma}{x|y} \right) $$

其中： $$ \begin{aligned} \boldsymbol{\mu}{x|y} &= \boldsymbol{\mu}x + \boldsymbol{\Sigma}{xy} \boldsymbol{\Sigma}{yy}^{-1} (\boldsymbol{y} - \boldsymbol{\mu}y) \ \boldsymbol{\Sigma}{x|y} &= \boldsymbol{\Sigma}{xx} - \boldsymbol{\Sigma}{xy} \boldsymbol{\Sigma}{yy}^{-1} \boldsymbol{\Sigma}{yx} \end{aligned} $$

这个公式的优美之处在于，条件均值 $\boldsymbol{\mu}{x|y}$ 是 $\boldsymbol{y}$ 的线性函数，而条件协方差 $\boldsymbol{\Sigma}{x|y}$ 与 $\boldsymbol{y}$ 的具体取值无关。这意味着，对于联合高斯分布，条件分布的计算是解析的、高效的。

为了更直观地理解这三个概念的关系，我们用一个简单的二维例子来可视化：

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal

# 设置联合高斯分布的参数
mu = np.array([0, 0])  # 均值
Sigma = np.array([[1, 0.8],  # 协方差矩阵，存在正相关
                  [0.8, 1]])

# 生成网格点用于绘制概率密度
x = np.linspace(-3, 3, 100)
y = np.linspace(-3, 3, 100)
X, Y = np.meshgrid(x, y)
pos = np.dstack((X, Y))

# 计算联合概率密度
rv = multivariate_normal(mu, Sigma)
Z_joint = rv.pdf(pos)

# 计算边际分布 (x的边际分布，即对y积分)
Z_marginal_x = multivariate_normal(mu[0], Sigma[0, 0]).pdf(x)

# 计算条件分布 (给定y=0.5时，x的条件分布)
y_cond = 0.5
mu_x_given_y = mu[0] + Sigma[0,1]/Sigma[1,1] * (y_cond - mu[1])
sigma2_x_given_y = Sigma[0,0] - Sigma[0,1]**2 / Sigma[1,1]
Z_conditional = multivariate_normal(mu_x_given_y, sigma2_x_given_y).pdf(x)

# 绘图
fig, axes = plt.subplots(1, 3, figsize=(15, 4))

# 1. 联合分布
contour = axes[0].contourf(X, Y, Z_joint, levels=20, cmap='viridis')
axes[0].set_xlabel('x')
axes[0].set_ylabel('y')
axes[0].set_title('联合分布 p(x, y)')
plt.colorbar(contour, ax=axes[0])

# 2. 边际分布 p(x)
axes[1].plot(x, Z_marginal_x, 'r-', linewidth=2)
axes[1].fill_between(x, Z_marginal_x, alpha=0.3, color='red')
axes[1].set_xlabel('x')
axes[1].set_ylabel('概率密度')
axes[1].set_title('边际分布 p(x)')
axes[1].grid(True, alpha=0.3)

# 3. 条件分布 p(x|y=0.5)
axes[2].plot(x, Z_conditional, 'b-', linewidth=2, label=f'p(x|y={y_cond})')
axes[2].axvline(mu_x_given_y, color='b', linestyle='--', label=f'条件均值={mu_x_given_y:.2f}')