Python实战：深入理解Gibbs采样及其代码实现16

作为一名专业的程序员，我们经常需要处理复杂的数据和模型，尤其是在机器学习、贝叶斯统计和概率图模型领域。在这些场景下，我们常常面临从高维、复杂或难以直接采样的概率分布中获取样本的问题。马尔可夫链蒙特卡洛（MCMC）方法正是解决这类问题的强大工具，而Gibbs采样作为MCMC家族中的重要一员，因其简洁性和高效性而备受青睐。

本文将深入探讨Gibbs采样的原理、算法，并通过Python代码实现一个经典的二元正态分布采样案例，帮助读者从理论到实践全面掌握Gibbs采样。

1. Gibbs采样：核心思想与背景

Gibbs采样是由统计物理学家Josiah Willard Gibbs命名的，但其在统计学领域的应用主要归功于Stuart Geman和Donald Geman在1984年提出的图像处理中的MCMC方法。它的核心思想是：当一个多维随机变量的联合概率分布难以直接采样时，如果它的所有“完全条件概率分布”（full conditional distributions）都相对容易采样，那么我们就可以通过依次从这些条件分布中采样来间接获得联合分布的样本。

想象一个具有N个随机变量的系统 X = (X1, X2, ..., XN)。我们想要从这个联合概率分布 p(X) 中采样。直接采样可能非常困难。Gibbs采样的策略是，在每个步骤中，我们只改变其中一个变量的值，而保持其他所有变量不变。具体来说，我们从 Xi 的条件分布 p(Xi | X-i) 中采样一个新的值 Xi'，其中 X-i 表示除 Xi 之外的所有其他变量。

通过这种迭代过程，Gibbs采样构造了一个马尔可夫链。在经过足够多的“预热”（burn-in）迭代后，这个马尔可夫链的平稳分布就是我们想要采样的目标联合分布 p(X)。

2. Gibbs采样的基本原理

要理解Gibbs采样的有效性，我们需要掌握两个关键概念：

2.1 马尔可夫链（Markov Chain）

Gibbs采样生成的序列 (X(0), X(1), X(2), ...) 构成了一个马尔可夫链。这意味着链中任何状态 X(t+1) 的概率只依赖于前一个状态 X(t)，而与更早的状态无关。Gibbs采样的巧妙之处在于，它设计的转移核（transition kernel）能够保证这个马尔可夫链是遍历的（irreducible）和非周期的（aperiodic），并且其平稳分布就是目标联合分布 p(X)。

2.2 完全条件概率分布（Full Conditional Distributions）

这是Gibbs采样的基石。对于变量 Xi，其完全条件概率分布定义为：

p(Xi | X1, ..., Xi-1, Xi+1, ..., XN)

简写为 p(Xi | X-i)。在实际应用中，计算这些条件分布通常比计算联合分布 p(X) 或其他边缘分布 p(Xi) 简单得多。例如，在很多贝叶斯模型中，如果选择了共轭先验，那么后验的完全条件分布通常会是标准分布（如正态、伽马、贝塔等），从而可以直接采样。

3. Gibbs采样的算法步骤

假设我们希望从 N 维联合分布 p(X1, X2, ..., XN) 中抽取 K 个样本，其算法步骤如下：

初始化： 为每个变量 X1, ..., XN 设置一个初始值 (X1(0), ..., XN(0))。这些初始值可以随机选取，但最好在分布的合理范围内。

迭代采样： 对于 t = 0 到 K-1 执行以下操作：

从 p(X1 | X2(t), X3(t), ..., XN(t)) 中采样一个新的 X1(t+1)。

从 p(X2 | X1(t+1), X3(t), ..., XN(t)) 中采样一个新的 X2(t+1)。

...

从 p(XN | X1(t+1), ..., XN-1(t+1)) 中采样一个新的 XN(t+1)。

注意：在采样 Xi 时，我们使用最新的已采样值来作为条件变量。也就是说，当采样 Xi 时，已经更新过的 X1 到 Xi-1 的值将用于条件，而 Xi+1 到 XN 仍然使用上一轮迭代 (t) 的值。

收集样本： 经过一定数量的“预热期”（burn-in period）迭代后，后续的样本 (X1(t), ..., XN(t)) 被视为近似来自目标联合分布的样本并被收集起来。

稀疏化（可选）： 为了减少样本之间的自相关性，可以每隔一定数量的迭代才收集一个样本（即进行“稀疏化”或“thinning”）。

4. 案例分析：二元正态分布的Gibbs采样

为了具体说明Gibbs采样，我们将使用一个经典的例子：从一个二元正态分布中采样。二元正态分布的联合概率密度函数形式比较复杂，但其条件分布是易于采样的。

4.1 数学推导：二元正态分布的条件分布

假设我们有一个二维随机变量 (X, Y) 服从二元正态分布 N(μ, Σ)，其中：

μ = [μx, μy]T 是均值向量

Σ = [[σx^2, ρσxσy], [ρσxσy, σy^2]] 是协方差矩阵，其中 ρ 是 X 和 Y 之间的相关系数。

则 X 在给定 Y=y 时的条件分布 p(X | Y=y) 仍然是正态分布，其均值和方差为：

E[X | Y=y] = μx + ρ * (σx / σy) * (y - μy)

Var[X | Y=y] = σx^2 * (1 - ρ^2)

同理，Y 在给定 X=x 时的条件分布 p(Y | X=x) 也是正态分布，其均值和方差为：

E[Y | X=x] = μy + ρ * (σy / σx) * (x - μx)

Var[Y | X=x] = σy^2 * (1 - ρ^2)

可以看到，给定另一个变量的值后，每个变量的条件分布都是一个简单的一元正态分布，这使得它们非常容易采样。

4.2 Python代码实现

我们将使用 `numpy` 进行数值计算，`matplotlib` 进行可视化。```python
import numpy as np
import as plt
import seaborn as sns
# 设置matplotlib中文显示
['-serif'] = ['SimHei'] # 指定默认字体
['axes.unicode_minus'] = False # 解决保存图像是负号'-'显示为方块的问题
def gibbs_sampler_bivariate_normal(
mean_vec, cov_matrix, num_iterations, initial_val=None, burn_in=1000, thinning=1
):
"""
实现二元正态分布的Gibbs采样。
参数:
mean_vec (): 均值向量 [mu_x, mu_y]。
cov_matrix (): 协方差矩阵 [[sigma_x^2, rho*sigma_x*sigma_y], [rho*sigma_x*sigma_y, sigma_y^2]]。
num_iterations (int): 总迭代次数。
initial_val (list/, optional): 初始值 [x0, y0]。如果为None，则随机初始化。
burn_in (int): 预热期迭代次数，在此期间的样本将被丢弃。
thinning (int): 稀疏化参数，每thinning个样本保留一个。
返回:
: 收集到的 (x, y) 样本数组。
"""
mu_x, mu_y = mean_vec[0], mean_vec[1]
sigma_x_sq = cov_matrix[0, 0]
sigma_y_sq = cov_matrix[1, 1]

# 计算协方差、标准差和相关系数
sigma_x = (sigma_x_sq)
sigma_y = (sigma_y_sq)
cov_xy = cov_matrix[0, 1] # 或 cov_matrix[1, 0]

# 防止除零错误，如果标准差为0则设为接近0的小值
if sigma_x == 0: sigma_x = 1e-9
if sigma_y == 0: sigma_y = 1e-9
rho = cov_xy / (sigma_x * sigma_y) if (sigma_x * sigma_y) != 0 else 0
# 初始化当前状态
if initial_val is None:
x_current = (mu_x, sigma_x)
y_current = (mu_y, sigma_y)
else:
x_current, y_current = initial_val[0], initial_val[1]
# 存储所有迭代的样本，用于可视化迹图
all_samples = []

# 存储经过burn_in和thinning处理的样本
collected_samples = []
for i in range(num_iterations):
# 1. 从 p(X | Y=y_current) 中采样 X
# 条件均值: mu_x_cond = mu_x + rho * (sigma_x / sigma_y) * (y_current - mu_y)
# 条件方差: sigma_x_cond_sq = sigma_x^2 * (1 - rho^2)
mu_x_cond = mu_x + rho * (sigma_x / sigma_y) * (y_current - mu_y)
sigma_x_cond = (sigma_x_sq * (1 - rho2))
x_current = (mu_x_cond, sigma_x_cond)
# 2. 从 p(Y | X=x_current) 中采样 Y
# 条件均值: mu_y_cond = mu_y + rho * (sigma_y / sigma_x) * (x_current - mu_x)
# 条件方差: sigma_y_cond_sq = sigma_y^2 * (1 - rho^2)
mu_y_cond = mu_y + rho * (sigma_y / sigma_x) * (x_current - mu_x)
sigma_y_cond = (sigma_y_sq * (1 - rho2))
y_current = (mu_y_cond, sigma_y_cond)

([x_current, y_current])
if i >= burn_in and (i - burn_in) % thinning == 0:
([x_current, y_current])
return (collected_samples), (all_samples)
# 定义二元正态分布的参数
mean = ([0, 0])
# 协方差矩阵，这里设置X和Y之间有较强的正相关性
covariance = ([[1.0, 0.8],
[0.8, 1.0]])
# 运行Gibbs采样
num_iterations = 20000 # 总迭代次数
burn_in_period = 2000 # 预热期
thinning_factor = 10 # 稀疏化
final_samples, all_iterations_samples = gibbs_sampler_bivariate_normal(
mean, covariance, num_iterations, burn_in=burn_in_period, thinning=thinning_factor
)
print(f"Gibbs采样完成。总共生成 {num_iterations} 个样本，经过预热和稀疏化后收集了 {len(final_samples)} 个有效样本。")
# ----------------------------------------------------
# 结果可视化
# ----------------------------------------------------
(figsize=(15, 6))
# 1. 迹图 (Trace Plot) - 显示马尔可夫链的收敛性
(1, 3, 1)
(all_iterations_samples[:burn_in_period, 0], color='red', alpha=0.5, label='X (Burn-in)')
(all_iterations_samples[burn_in_period:, 0], color='blue', alpha=0.5, label='X (Collected)')
(all_iterations_samples[:burn_in_period, 1], color='orange', alpha=0.5, label='Y (Burn-in)')
(all_iterations_samples[burn_in_period:, 1], color='green', alpha=0.5, label='Y (Collected)')
('Gibbs采样迹图 (Trace Plot)')
('迭代次数')
('样本值')
()
(True, linestyle='--', alpha=0.6)
# 2. 边缘分布直方图 (Marginal Histograms) - 检查是否符合目标分布的边缘
(1, 3, 2)
(final_samples[:, 0], kde=True, color='skyblue', label='X 边缘分布', stat='density')
(final_samples[:, 1], kde=True, color='lightcoral', label='Y 边缘分布', stat='density')
('Gibbs采样边缘分布')
('样本值')
('密度')
()
(True, linestyle='--', alpha=0.6)
# 3. 联合分布散点图 (Joint Distribution Scatter Plot) - 检查联合分布形态
(1, 3, 3)
(final_samples[:, 0], final_samples[:, 1], alpha=0.3, s=5, color='purple')
('Gibbs采样联合分布')
('X')
('Y')
(True, linestyle='--', alpha=0.6)
plt.tight_layout()
()
# 验证采样结果的统计特性
print("采样结果的均值:")
print(f"X的均值: {(final_samples[:, 0]):.4f}, 目标均值: {mean[0]:.4f}")
print(f"Y的均值: {(final_samples[:, 1]):.4f}, 目标均值: {mean[1]:.4f}")
print("采样结果的协方差矩阵:")
sampled_cov = (final_samples[:, 0], final_samples[:, 1])
print(sampled_cov)
print("目标协方差矩阵:")
print(covariance)
```

4.3 代码解释

1. 参数设置： 定义了二元正态分布的均值向量 `mean_vec` 和协方差矩阵 `cov_matrix`。我们特意设置了0.8的相关系数，以便观察采样结果的形状。

2. 条件分布参数计算： 在 `gibbs_sampler_bivariate_normal` 函数内部，首先从协方差矩阵中提取出方差 `sigma_x_sq`, `sigma_y_sq` 和协方差 `cov_xy`，然后计算出相关系数 `rho`。接着，利用上一节的数学推导公式，计算出 X 和 Y 在给定另一个变量时的条件均值和条件标准差。这是Gibbs采样最核心的部分。

3. 初始化： 随机初始化 `x_current` 和 `y_current`，或者使用用户提供的初始值。良好的初始化有助于链更快收敛，尽管理论上任何值都可以。

4. 迭代采样：

在每次迭代中，首先利用当前 `y_current` 的值计算 X 的条件均值 `mu_x_cond` 和条件标准差 `sigma_x_cond`。然后使用 `()` 从这个条件正态分布中采样一个新的 `x_current`。
接着，利用刚刚更新的 `x_current` 的值计算 Y 的条件均值 `mu_y_cond` 和条件标准差 `sigma_y_cond`。再次使用 `()` 从这个条件正态分布中采样一个新的 `y_current`。
这个“先更新X再更新Y”的顺序是Gibbs采样的精髓，确保了每次采样都基于最新的信息。

5. 预热和稀疏化：

`burn_in` 参数用于指定丢弃链开头样本的数量。在预热期内，马尔可夫链可能还没有收敛到目标分布，因此这些样本通常不具有代表性。
`thinning` 参数用于稀疏化采样结果，每隔 `thinning` 个样本才保留一个。这有助于减少样本之间的自相关性，使收集到的样本更接近独立同分布。

6. 结果可视化：

迹图 (Trace Plot)： 绘制了 X 和 Y 变量随迭代次数变化的轨迹。理想的迹图应该在预热期后呈现出“毛茸茸”的随机游走模式，且没有明显的趋势或周期性，表明链已收敛。
边缘分布直方图： 使用 `` 绘制 X 和 Y 变量的边缘分布直方图和核密度估计 (KDE)。这些图应该近似于目标二元正态分布的边缘正态分布。
联合分布散点图： 绘制 X 和 Y 的散点图，展示它们的联合分布。此图应该呈现出与目标二元正态分布相似的椭圆形状和相关性。

5. Gibbs采样的优点与局限性