深度学习学习教程，从入门到精通，蒙特卡罗方法 — 完整知识点与代码案例(17）
·
蒙特卡罗方法 — 完整知识点与代码案例
一、采样和蒙特卡罗方法
1.1 基本概念
蒙特卡罗方法（Monte Carlo Methods） 是一类基于随机采样来近似计算数学问题数值解的方法。其核心思想是：通过从目标概率分布中抽取大量样本，利用样本的统计特性来估计期望、积分、概率等量。
关键公式：
Ex∼P[f(x)]=∫f(x)P(x)dx≈1N∑i=1Nf(x(i)),x(i)∼P(x)E_{x \sim P}[f(x)] = \int f(x) P(x) dx \approx \frac{1}{N} \sum_{i=1}^{N} f(x^{(i)}), \quad x^{(i)} \sim P(x)Ex∼P[f(x)]=∫f(x)P(x)dx≈N1i=1∑Nf(x(i)),x(i)∼P(x)
当样本数 N→∞N \to \inftyN→∞ 时，根据大数定律，样本均值收敛于真实期望。
1.2 基本采样方法
（1）逆变换采样（Inverse Transform Sampling）
原理：若已知累积分布函数（CDF）F(x)F(x)F(x) 及其逆函数 F−1(u)F^{-1}(u)F−1(u)，则 X=F−1(U)X = F^{-1}(U)X=F−1(U)，其中 U∼Uniform(0,1)U \sim \text{Uniform}(0,1)U∼Uniform(0,1)。
步骤：
- 生成均匀随机数 u∼U(0,1)u \sim U(0,1)u∼U(0,1)
- 计算 x=F−1(u)x = F^{-1}(u)x=F−1(u)
- xxx 即为目标分布的一个样本
（2）拒绝采样（Rejection Sampling）
原理：当目标分布 p(x)p(x)p(x) 难以直接采样时，选择一个容易采样的提议分布 q(x)q(x)q(x)，使得存在常数 MMM 满足 Mq(x)≥p(x)Mq(x) \geq p(x)Mq(x)≥p(x) 对所有 xxx 成立。
步骤：
- 从 q(x)q(x)q(x) 中采样 x∗x^*x∗
- 从 U(0,1)U(0,1)U(0,1) 中采样 uuu
- 若 u≤p(x∗)Mq(x∗)u \leq \frac{p(x^*)}{Mq(x^*)}u≤Mq(x∗)p(x∗)，接受 x∗x^*x∗；否则拒绝，回到步骤1
1.3 代码案例
import numpy as np
import matplotlib.pyplot as plt
# ============================================================
# 案例1：逆变换采样 —— 从指数分布中采样
# ============================================================
# 目标分布：指数分布 Exponential(lambda=1.5)
# PDF: p(x) = lambda * exp(-lambda * x), x >= 0
# CDF: F(x) = 1 - exp(-lambda * x)
# 逆CDF: F^{-1}(u) = -ln(1 - u) / lambda
lambda_param = 1.5 # 指数分布的参数 lambda
n_samples = 10000 # 需要生成的样本数量
# 步骤1：从均匀分布 U(0,1) 生成随机数
u = np.random.uniform(0, 1, n_samples) # 生成 n_samples 个 [0,1) 之间的均匀随机数
# 步骤2：通过逆CDF变换，将均匀随机数转换为指数分布样本
samples_exp = -np.log(1 - u) / lambda_param # 应用逆变换公式 x = -ln(1-u)/lambda
# 可视化：对比采样结果与真实PDF
x_range = np.linspace(0, 6, 500) # 定义x轴的绘图范围
true_pdf = lambda_param * np.exp(-lambda_param * x_range) # 计算真实概率密度函数值
plt.figure(figsize=(12, 5))
# 子图1：直方图 vs 真实PDF
plt.subplot(1, 2, 1)
plt.hist(samples_exp, bins=80, density=True, alpha=0.7, color='steelblue',
label='采样直方图') # 绘制采样结果的归一化直方图
plt.plot(x_range, true_pdf, 'r-', linewidth=2,
label=f'真实PDF (λ={lambda_param})') # 绘制真实概率密度曲线
plt.xlabel('x') # 设置x轴标签
plt.ylabel('概率密度') # 设置y轴标签
plt.title('逆变换采样：指数分布') # 设置图标题
plt.legend() # 显示图例
# 子图2：均匀分布输入 vs 指数分布输出
plt.subplot(1, 2, 2)
plt.scatter(u[:1000], samples_exp[:1000], s=1, alpha=0.3, color='coral')
plt.xlabel('u ~ Uniform(0,1)') # x轴为原始均匀分布
plt.ylabel('x = F⁻¹(u)') # y轴为变换后的指数分布样本
plt.title('逆变换函数映射关系') # 展示映射过程
plt.tight_layout() # 自动调整子图间距
plt.show() # 显示图形
# ============================================================
# 案例2：拒绝采样 —— 从复杂分布中采样
# ============================================================
# 目标分布 p(x)：混合高斯分布（双峰）
# p(x) = 0.3 * N(-2, 0.5^2) + 0.7 * N(3, 1.0^2)
# 这个分布不容易直接采样，需要用拒绝采样
def target_pdf(x):
"""
计算目标分布（混合高斯）的概率密度值
参数:
x: 输入点（标量或数组）
返回:
对应的概率密度值
"""
# 第一个高斯分量：均值=-2, 标准差=0.5, 权重=0.3
comp1 = 0.3 * np.exp(-0.5 * ((x + 2) / 0.5) ** 2) / (0.5 * np.sqrt(2 * np.pi))
# 第二个高斯分量：均值=3, 标准差=1.0, 权重=0.7
comp2 = 0.7 * np.exp(-0.5 * ((x - 3) / 1.0) ** 2) / (1.0 * np.sqrt(2 * np.pi))
return comp1 + comp2 # 返回两个分量之和
# 提议分布 q(x)：均匀分布 U(-6, 8)
# q(x) = 1 / (8 - (-6)) = 1/14, 当 x ∈ [-6, 8]
a, b = -6, 8 # 均匀分布的范围
M = 15.0 # 包络常数 M，使得 M*q(x) >= p(x) 对所有x成立
n_samples_rejection = 50000 # 尝试采样的最大次数
accepted_samples = [] # 存储接受的样本
n_trials = 0 # 记录总尝试次数
while len(accepted_samples) < 10000: # 直到接受10000个样本
# 步骤1：从提议分布 q(x) = Uniform(-6, 8) 中采样候选点
x_star = np.random.uniform(a, b) # 生成候选样本
# 步骤2：从 Uniform(0, 1) 中采样用于接受/拒绝判断的随机数
u_accept = np.random.uniform(0, 1) # 生成 [0,1) 均匀随机数
# 步骤3：计算接受概率 = p(x*) / (M * q(x*))
# 其中 q(x*) = 1/(b-a) = 1/14 为常数
acceptance_ratio = target_pdf(x_star) / (M * (1.0 / (b - a))) # 计算接受比率
# 步骤4：接受/拒绝判断
if u_accept <= acceptance_ratio: # 若 u <= 接受概率，则接受
accepted_samples.append(x_star) # 将候选样本加入接受列表
n_trials += 1 # 累加尝试次数
accepted_samples = np.array(accepted_samples) # 转换为numpy数组
# 计算接受率（有效采样效率）
acceptance_rate = len(accepted_samples) / n_trials # 接受率 = 接受样本数 / 总尝试次数
print(f"拒绝采样接受率: {acceptance_rate:.4f}") # 输出接受率
print(f"总尝试次数: {n_trials}") # 输出总尝试次数
# 可视化拒绝采样结果
x_plot = np.linspace(-6, 8, 1000) # 定义绘图的x范围
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.hist(accepted_samples, bins=100, density=True, alpha=0.7, color='mediumseagreen',
label='拒绝采样结果') # 绘制接受样本的直方图
plt.plot(x_plot, target_pdf(x_plot), 'r-', linewidth=2,
label='目标分布 p(x)') # 绘制真实目标分布曲线
plt.plot(x_plot, M * np.ones_like(x_plot) / (b - a), 'b--', linewidth=1.5,
label=f'包络 M·q(x) = {M/(b-a):.2f}') # 绘制包络线
plt.xlabel('x')
plt.ylabel('概率密度')
plt.title('拒绝采样：混合高斯分布')
plt.legend()
# 子图2：展示接受与拒绝的过程（仅前500次尝试）
plt.subplot(1, 2, 2)
# 重新进行少量采样用于可视化过程
np.random.seed(42) # 固定随机种子以保证可复现
x_candidates = np.random.uniform(a, b, 500) # 生成500个候选样本
u_vals = np.random.uniform(0, 1, 500) # 生成500个判断用随机数
ratios = target_pdf(x_candidates) / (M / (b - a)) # 计算每个候选的接受比率
accepted_mask = u_vals <= ratios # 布尔掩码：True表示接受
plt.scatter(x_candidates[accepted_mask], u_vals[accepted_mask],
s=5, c='green', alpha=0.5, label='接受') # 绘制接受的点（绿色）
plt.scatter(x_candidates[~accepted_mask], u_vals[~accepted_mask],
s=5, c='red', alpha=0.3, label='拒绝') # 绘制拒绝的点（红色）
plt.plot(x_plot, ratios[np.searchsorted(x_candidates, x_plot).clip(0, 499)],
'b-', alpha=0.3) # 绘制接受概率曲线的近似
plt.xlabel('候选样本 x*')
plt.ylabel('u (用于判断的随机数)')
plt.title('拒绝采样过程可视化')
plt.legend()
plt.tight_layout()
plt.show()
二、重要采样（Importance Sampling）
2.1 核心概念
重要采样是蒙特卡罗方法中最常用的技术之一。其核心思想是：当无法从目标分布 p(x)p(x)p(x) 直接采样时，可以从一个提议分布 q(x)q(x)q(x) 中采样，并通过重要性权重来修正偏差。
公式推导：
Ex∼p[f(x)]=∫f(x)p(x)dx=∫f(x)p(x)q(x)q(x)dx=Ex∼q[f(x)p(x)q(x)]E_{x \sim p}[f(x)] = \int f(x) p(x) dx = \int f(x) \frac{p(x)}{q(x)} q(x) dx = E_{x \sim q}\left[f(x) \frac{p(x)}{q(x)}\right]Ex∼p[f(x)]=∫f(x)p(x)dx=∫f(x)q(x)p(x)q(x)dx=Ex∼q[f(x)q(x)p(x)]
估计量：
μ^=1N∑i=1Nf(x(i))p(x(i))q(x(i)),x(i)∼q(x)\hat{\mu} = \frac{1}{N} \sum_{i=1}^{N} f(x^{(i)}) \frac{p(x^{(i)})}{q(x^{(i)})}, \quad x^{(i)} \sim q(x)μ^=N1i=1∑Nf(x(i))q(x(i))p(x(i)),x(i)∼q(x)
其中 w(x(i))=p(x(i))q(x(i))w(x^{(i)}) = \frac{p(x^{(i)})}{q(x^{(i)})}w(x(i))=q(x(i))p(x(i)) 称为重要性权重（importance weight）。
2.2 提议分布的选择
最优提议分布（最小化方差）：
q∗(x)=∣f(x)∣p(x)∫∣f(x)∣p(x)dxq^*(x) = \frac{|f(x)| p(x)}{\int |f(x)| p(x) dx}q∗(x)=∫∣f(x)∣p(x)dx∣f(x)∣p(x)
实践中，q(x)q(x)q(x) 应满足：
- q(x)>0q(x) > 0q(x)>0 当 f(x)p(x)≠0f(x)p(x) \neq 0f(x)p(x)=0
- q(x)q(x)q(x) 的尾部应比 p(x)p(x)p(x) 更重（heavy-tailed），避免权重爆炸
- q(x)q(x)q(x) 应易于采样
2.3 标准化重要采样（Self-Normalized Importance Sampling）
当 p(x)p(x)p(x) 仅已知到一个归一化常数时（即只知道 p~(x)\tilde{p}(x)p~(x)），使用标准化重要采样：
μ^=∑i=1Nf(x(i))w~(x(i))∑i=1Nw~(x(i)),w~(x(i))=p~(x(i))q(x(i))\hat{\mu} = \frac{\sum_{i=1}^{N} f(x^{(i)}) \tilde{w}(x^{(i)})}{\sum_{i=1}^{N} \tilde{w}(x^{(i)})}, \quad \tilde{w}(x^{(i)}) = \frac{\tilde{p}(x^{(i)})}{q(x^{(i)})}μ^=∑i=1Nw~(x(i))∑i=1Nf(x(i))w~(x(i)),w~(x(i))=q(x(i))p~(x(i))
2.4 代码案例
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
# ============================================================
# 案例1：基础重要采样 —— 估计高斯分布尾部概率
# ============================================================
# 问题：计算 P(X > 4), X ~ N(0, 1)
# 这是一个极小概率事件，直接蒙特卡罗估计效率很低
np.random.seed(42) # 固定随机种子
n_samples = 10000 # 样本数量
# --- 方法1：直接蒙特卡罗采样 ---
samples_direct = np.random.normal(0, 1, n_samples) # 从标准正态分布直接采样
# 计算有多少样本落在 x > 4 区域
estimate_direct = np.mean(samples_direct > 4) # 直接估计 P(X > 4)
print(f"直接蒙特卡罗估计 P(X>4): {estimate_direct}") # 大部分时候为0，估计极不准确
# --- 方法2：重要采样 ---
# 提议分布 q(x) = N(4, 1)，将采样重心移到感兴趣区域
mu_q = 4.0 # 提议分布的均值
sigma_q = 1.0 # 提议分布的标准差
# 步骤1：从提议分布 q(x) 中采样
samples_is = np.random.normal(mu_q, sigma_q, n_samples) # 从 N(4,1) 采样
# 步骤2：计算重要性权重 w(x) = p(x) / q(x)
# p(x) = N(0,1) 的PDF值
p_values = stats.norm.pdf(samples_is, loc=0, scale=1) # 目标分布 p(x) = N(0,1)
# q(x) = N(4,1) 的PDF值
q_values = stats.norm.pdf(samples_is, loc=mu_q, scale=sigma_q) # 提议分布 q(x) = N(4,1)
# 重要性权重
weights = p_values / q_values # w_i = p(x_i) / q(x_i)
# 步骤3：计算指示函数 f(x) = I(x > 4)
indicator = (samples_is > 4).astype(float) # f(x) = 1 if x > 4, else 0
# 步骤4：加权估计
estimate_is = np.mean(indicator * weights) # E_p[I(x>4)] ≈ mean(f(x)*w(x))
print(f"重要采样估计 P(X>4): {estimate_is:.6f}") # 输出重要采样估计值
# 真实值（通过scipy精确计算）
true_value = 1 - stats.norm.cdf(4) # 真实的 P(X > 4)
print(f"真实值 P(X>4): {true_value:.6f}") # 输出精确值作为对比
# 方差分析
# 有效样本数（Effective Sample Size, ESS）
ESS = (np.sum(weights) ** 2) / np.sum(weights ** 2) # ESS = (sum w)^2 / sum(w^2)
print(f"有效样本数 ESS: {ESS:.1f} / {n_samples}") # ESS越接近N越好
# ============================================================
# 案例2：标准化重要采样 —— 估计未归一化分布的期望
# ============================================================
# 场景：目标分布 p(x) 只知道未归一化的形式 tilde_p(x)
# 例如：tilde_p(x) = exp(-x^2/2) * (1 + sin(2x))^2
# 我们想计算 E_p[x] 和 E_p[x^2]
def unnormalized_log_pdf(x):
"""
计算未归一化对数概率密度
参数:
x: 输入值
返回:
log(tilde_p(x))
"""
return -x**2 / 2 + 2 * np.log(np.abs(1 + np.sin(2 * x)) + 1e-10) # 加小常数避免log(0)
def unnormalized_pdf(x):
"""
计算未归一化概率密度（直接形式）
参数:
x: 输入值
返回:
tilde_p(x)
"""
return np.exp(-x**2 / 2) * (1 + np.sin(2 * x))**2 # 未归一化的密度函数
np.random.seed(42) # 固定随机种子
n_samples_snis = 50000 # 采样数量
# 提议分布 q(x) = N(0, 2)，选择较宽的高斯以覆盖目标分布
mu_proposal = 0.0 # 提议分布均值
sigma_proposal = 2.0 # 提议分布标准差（较宽）
# 步骤1：从提议分布中采样
x_samples = np.random.normal(mu_proposal, sigma_proposal, n_samples_snis) # 从q(x)采样
# 步骤2：计算未归一化重要性权重
p_tilde = unnormalized_pdf(x_samples) # 计算 tilde_p(x_i)
q_vals = stats.norm.pdf(x_samples, mu_proposal, sigma_proposal) # 计算 q(x_i)
unnormalized_weights = p_tilde / q_vals # w_tilde_i = tilde_p(x_i) / q(x_i)
# 步骤3：标准化重要采样估计
# E_p[x] 的估计
x_expectation = np.sum(x_samples * unnormalized_weights) / np.sum(unnormalized_weights)
print(f"标准化重要采样估计 E_p[x]: {x_expectation:.6f}") # 期望的估计值
# E_p[x^2] 的估计
x2_expectation = np.sum(x_samples**2 * unnormalized_weights) / np.sum(unnormalized_weights)
print(f"标准化重要采样估计 E_p[x^2]: {x2_expectation:.6f}") # 二阶矩的估计值
# 可视化：目标分布 vs 提议分布 vs 采样权重
x_plot = np.linspace(-6, 6, 1000) # 绘图范围
p_tilde_plot = unnormalized_pdf(x_plot) # 未归一化目标分布
q_plot = stats.norm.pdf(x_plot, mu_proposal, sigma_proposal) # 提议分布
# 归一化目标分布用于对比（数值积分）
from scipy.integrate import quad
normalization_const, _ = quad(unnormalized_pdf, -20, 20) # 数值计算归一化常数
p_normalized = p_tilde_plot / normalization_const # 归一化后的目标分布
plt.figure(figsize=(14, 5))
plt.subplot(1, 2, 1)
plt.plot(x_plot, p_normalized, 'r-', linewidth=2, label='目标分布 p(x)（归一化后）')
plt.plot(x_plot, q_plot, 'b--', linewidth=2, label='提议分布 q(x) = N(0, 2)')
plt.fill_between(x_plot, 0, p_normalized, alpha=0.1, color='red')
plt.fill_between(x_plot, 0, q_plot, alpha=0.1, color='blue')
plt.xlabel('x')
plt.ylabel('概率密度')
plt.title('目标分布 vs 提议分布')
plt.legend()
plt.grid(True, alpha=0.3)
# 子图2：样本点的权重分布
plt.subplot(1, 2, 2)
# 按权重大小着色
norm_weights = unnormalized_weights / np.max(unnormalized_weights) # 归一化权重用于颜色映射
plt.scatter(x_samples[:2000], unnormalized_weights[:2000],
c=norm_weights[:2000], cmap='hot', s=3, alpha=0.5) # 权重散点图
plt.xlabel('样本值 x')
plt.ylabel('未归一化权重 w̃')
plt.title('重要性权重分布')
plt.colorbar(label='归一化权重')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# ============================================================
# 案例3：重要采样在强化学习中的应用 —— 策略评估
# ============================================================
# 用重要采样从行为策略(behavior policy)的轨迹评估目标策略(target policy)
# 这是off-policy方法的基础
def behavior_policy(state):
"""
行为策略：给定状态，返回各动作的概率（较随机的策略）
参数:
state: 状态索引
返回:
各动作的概率分布
"""
n_actions = 4 # 动作数量
probs = np.ones(n_actions) / n_actions # 均匀策略：各动作概率相等
return probs
def target_policy(state):
"""
目标策略：给定状态，返回各动作的概率（较确定性的策略）
参数:
state: 状态索引
返回:
各动作的概率分布
"""
n_actions = 4 # 动作数量
probs = np.zeros(n_actions) # 初始化概率为0
probs[state % n_actions] = 0.7 # 偏好某个动作，概率0.7
probs += 0.3 / n_actions # 其余动作分配剩余概率
return probs
# 模拟收集一条轨迹
np.random.seed(42)
n_steps = 20 # 轨迹长度
states = np.random.randint(0, 4, n_steps) # 随机生成状态序列
rewards = np.random.randn(n_steps) * 0.5 + 1 # 随机生成奖励（均值为1）
# 使用行为策略收集的动作
actions = [] # 存储每步的动作
for t in range(n_steps): # 遍历每个时间步
action_probs = behavior_policy(states[t]) # 获取行为策略的动作概率
action = np.random.choice(4, p=action_probs) # 按行为策略采样动作
actions.append(action) # 记录动作
actions = np.array(actions) # 转为numpy数组
# 计算每个时间步的重要性权重
# w_t = pi_target(a_t | s_t) / pi_behavior(a_t | s_t)
importance_weights = np.ones(n_steps) # 初始化权重为1
for t in range(n_steps): # 遍历每个时间步
pi_target = target_policy(states[t])[actions[t]] # 目标策略选择该动作的概率
pi_behavior = behavior_policy(states[t])[actions[t]] # 行为策略选择该动作的概率
importance_weights[t] = pi_target / pi_behavior # 计算重要性权重
# 普通重要采样估计每个状态的期望回报
# 在这个简单例子中，我们估计 E[R] under target policy
# 乘积权重（从轨迹开始到当前步）
cumulative_weights = np.cumprod(importance_weights) # 累积乘积权重
# 普通重要采样的回报估计
V_estimate = np.zeros(n_steps) # 存储每个时间步的价值估计
for t in range(n_steps): # 遍历每个时间步
# 从时间步t开始的折扣回报
gamma = 0.9 # 折扣因子
returns_t = sum(gamma**(k-t) * rewards[k] for k in range(t, n_steps)) # 计算折扣回报
V_estimate[t] = cumulative_weights[t] * returns_t # 加权回报
print("各时间步的重要性权重:", np.round(importance_weights, 4))
print("各时间步的累积权重:", np.round(cumulative_weights, 4))
print(f"普通IS估计的平均价值: {np.mean(V_estimate):.4f}")
# 对比：如果使用截断权重（减少方差）
truncation_threshold = 5.0 # 权重截断阈值
cumulative_weights_clipped = np.minimum(cumulative_weights, truncation_threshold) # 截断权重
V_estimate_clipped = np.zeros(n_steps) # 存储截断版本的价值估计
for t in range(n_steps):
gamma = 0.9
returns_t = sum(gamma**(k-t) * rewards[k] for k in range(t, n_steps))
V_estimate_clipped[t] = cumulative_weights_clipped[t] * returns_t # 使用截断权重
print(f"截断IS估计的平均价值: {np.mean(V_estimate_clipped):.4f}")
三、马尔可夫链蒙特卡罗方法（MCMC）
3.1 核心概念
MCMC（Markov Chain Monte Carlo） 是一类构造马尔可夫链、使其平稳分布等于目标分布 p(x)p(x)p(x) 的采样方法。
关键定理——细致平衡条件（Detailed Balance）：
p(x)T(x′∣x)=p(x′)T(x∣x′)p(x) T(x'|x) = p(x') T(x|x')p(x)T(x′∣x)=p(x′)T(x∣x′)
其中 T(x′∣x)T(x'|x)T(x′∣x) 是从状态 xxx 转移到 x′x'x′ 的转移概率。满足细致平衡条件的转移核 TTT 的平稳分布即为 p(x)p(x)p(x)。
3.2 Metropolis-Hastings 算法
MH算法是最经典的MCMC方法。
算法流程：
- 初始化 x(0)x^{(0)}x(0)
- 对于 t=0,1,2,…,T−1t = 0, 1, 2, \ldots, T-1t=0,1,2,…,T−1：
- 从提议分布 q(x′∣x(t))q(x'|x^{(t)})q(x′∣x(t)) 中采样候选点 x′x'x′
- 计算接受概率：
α=min(1,p(x′)q(x(t)∣x′)p(x(t))q(x′∣x(t)))\alpha = \min\left(1, \frac{p(x') q(x^{(t)}|x')}{p(x^{(t)}) q(x'|x^{(t)})}\right)α=min(1,p(x(t))q(x′∣x(t))p(x′)q(x(t)∣x′)) - 以概率 α\alphaα 接受：x(t+1)=x′x^{(t+1)} = x'x(t+1)=x′，否则 x(t+1)=x(t)x^{(t+1)} = x^{(t)}x(t+1)=x(t)
3.3 代码案例
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
# ============================================================
# 案例1：Metropolis-Hastings 采样 —— 从二维混合高斯中采样
# ============================================================
def target_distribution_2d(x, y):
"""
计算目标分布（二维混合高斯）的概率密度值
p(x,y) = 0.4 * N(mu1, Sigma1) + 0.6 * N(mu2, Sigma2)
参数:
x, y: 二维坐标
返回:
概率密度值
"""
# 第一个高斯分量：均值=[-2, -2], 标准差=[1, 1]
comp1 = 0.4 * np.exp(-0.5 * ((x+2)**2 + (y+2)**2)) # 各向同性高斯
# 第二个高斯分量：均值=[3, 3], 标准差=[1.5, 1.5]
comp2 = 0.6 * np.exp(-0.5 * ((x-3)**2/1.5**2 + (y-3)**2/1.5**2))
return comp1 + comp2 # 返回混合密度
def metropolis_hastings_2d(n_samples, proposal_std=1.0, initial_point=None):
"""
二维Metropolis-Hastings采样器
参数:
n_samples: 需要生成的样本数量
proposal_std: 提议分布（各向同性高斯）的标准差
initial_point: 初始点，默认为原点
返回:
samples: 采样结果数组, shape=(n_samples, 2)
accept_rate: 接受率
"""
if initial_point is None: # 若未指定初始点
initial_point = np.array([0.0, 0.0]) # 默认从原点开始
samples = np.zeros((n_samples, 2)) # 初始化样本数组
samples[0] = initial_point # 设置初始样本
n_accepted = 0 # 接受次数计数器
for i in range(1, n_samples): # 从第2个样本开始迭代
current = samples[i-1] # 当前状态 x_t
# 步骤1：从对称的高斯提议分布 q(x'|x) = N(x, sigma^2 * I) 中采样候选点
proposal = current + np.random.normal(0, proposal_std, size=2) # 随机游走提议
# 步骤2：计算接受概率 alpha
# 由于提议分布是对称的（高斯），q(x'|x) = q(x|x')
# 因此 MH 比率简化为 p(x') / p(x)
p_current = target_distribution_2d(current[0], current[1]) # p(x_t)
p_proposal = target_distribution_2d(proposal[0], proposal[1]) # p(x')
if p_current > 0: # 避免除以零
acceptance_ratio = p_proposal / p_current # MH比率（对称提议简化）
else:
acceptance_ratio = 1.0 # 若当前密度为0，直接接受
# 步骤3：接受/拒绝
alpha = min(1.0, acceptance_ratio) # 接受概率，截断在1
u = np.random.uniform(0, 1) # 生成均匀随机数
if u < alpha: # 若 u < alpha，接受候选
samples[i] = proposal # 新状态为候选点
n_accepted += 1 # 接受计数+1
else: # 否则拒绝
samples[i] = current # 新状态保持不变（拒绝移动）
accept_rate = n_accepted / (n_samples - 1) # 计算接受率
return samples, accept_rate # 返回样本和接受率
np.random.seed(42) # 固定随机种子
# 运行MH采样器
n_mcmc_samples = 30000 # 总采样数
samples_mh, acc_rate = metropolis_hastings_2d(
n_samples=n_mcmc_samples, # 样本数量
proposal_std=1.5, # 提议分布标准差
initial_point=np.array([-2.0, -2.0]) # 从第一个模态附近开始
)
print(f"MH采样接受率: {acc_rate:.4f}") # 输出接受率
# 去除burn-in期的样本（前20%的样本作为预热）
burn_in = int(0.2 * n_mcmc_samples) # burn-in长度
samples_after_burnin = samples_mh[burn_in:] # 丢弃burn-in样本
# 可视化
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
# 子图1：采样轨迹
axes[0].plot(samples_mh[:5000, 0], samples_mh[:5000, 1],
linewidth=0.3, alpha=0.5, color='blue') # 绘制前5000步的轨迹
axes[0].scatter(samples_mh[0, 0], samples_mh[0, 1],
c='green', s=100, zorder=5, label='起始点') # 标记起始点
axes[0].set_xlabel('x')
axes[0].set_ylabel('y')
axes[0].set_title('MH采样轨迹（前5000步）')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
# 子图2：采样散点图（去burn-in后）
axes[1].scatter(samples_after_burnin[:, 0], samples_after_burnin[:, 1],
s=1, alpha=0.1, c='coral') # 绘制所有样本点
axes[1].set_xlabel('x')
axes[1].set_ylabel('y')
axes[1].set_title(f'MH采样结果（去burn-in后, n={len(samples_after_burnin)}）')
axes[1].set_aspect('equal')
axes[1].grid(True, alpha=0.3)
# 子图3：核密度估计的等高线图
from scipy.stats import gaussian_kde
kde = gaussian_kde(samples_after_burnin.T) # 对采样结果做KDE
x_grid = np.linspace(-6, 7, 100) # x轴网格
y_grid = np.linspace(-6, 7, 100) # y轴网格
X, Y = np.meshgrid(x_grid, y_grid) # 生成二维网格
positions = np.vstack([X.ravel(), Y.ravel()]) # 展平为2行矩阵
Z = kde(positions).reshape(X.shape) # 计算KDE并reshape
axes[2].contourf(X, Y, Z, levels=20, cmap='YlOrRd') # 填充等高线图
axes[2].set_xlabel('x')
axes[2].set_ylabel('y')
axes[2].set_title('核密度估计等高线')
axes[2].set_aspect('equal')
plt.tight_layout()
plt.show()
# ============================================================
# 案例2：MH算法的诊断 —— 迹图和自相关函数
# ============================================================
def compute_autocorrelation(chain, max_lag=100):
"""
计算马尔可夫链的自相关函数（ACF）
参数:
chain: 一维采样链
max_lag: 最大滞后值
返回:
acf: 自相关系数数组
"""
n = len(chain) # 链的长度
mean = np.mean(chain) # 链的均值
var = np.var(chain) # 链的方差
acf = np.zeros(max_lag) # 初始化ACF数组
for lag in range(max_lag): # 遍历每个滞后值
if var > 0: # 避免除以零
# 计算滞后lag的自相关系数
cov = np.mean((chain[:n-lag] - mean) * (chain[lag:] - mean)) # 自协方差
acf[lag] = cov / var # 归一化为自相关系数
else:
acf[lag] = 0 # 方差为零时ACF为0
return acf # 返回ACF数组
# 使用不同proposal标准差运行多次MH
proposal_stds = [0.1, 1.5, 10.0] # 三种不同的提议标准差
results = {} # 存储各次结果
for pstd in proposal_stds: # 遍历每种标准差
samples_diag, ar = metropolis_hastings_2d(
n_samples=20000, # 样本数
proposal_std=pstd, # 提议标准差
initial_point=np.array([0.0, 0.0]) # 初始点
)
results[pstd] = { # 存储结果
'samples': samples_diag, # 采样链
'accept_rate': ar, # 接受率
'acf_x': compute_autocorrelation(samples_diag[:, 0], 100) # x分量的ACF
}
print(f"proposal_std={pstd:.1f}, 接受率={ar:.4f}") # 输出诊断信息
fig, axes = plt.subplots(len(proposal_stds), 3, figsize=(18, 4*len(proposal_stds)))
for idx, pstd in enumerate(proposal_stds): # 遍历每种提议标准差
s = results[pstd]['samples'] # 获取采样结果
# 迹图（trace plot）
axes[idx, 0].plot(s[:5000, 0], linewidth=0.5) # 绘制x分量的迹图
axes[idx, 0].set_title(f'σ_proposal={pstd}, AR={results[pstd]["accept_rate"]:.2f}')
axes[idx, 0].set_ylabel('x值') # y轴标签
# 自相关函数
axes[idx, 1].bar(range(50), results[pstd]['acf_x'][:50]) # 绘制ACF柱状图
axes[idx, 1].set_title(f'自相关函数 (lag=0~49)')
axes[idx, 1].axhline(y=0, color='red', linestyle='--') # 零线
# 采样直方图
axes[idx, 2].hist(s[2000:, 0], bins=80, density=True, alpha=0.7, color='steelblue')
axes[idx, 2].set_title('x分量的边际分布')
axes[idx, 2].set_xlabel('x')
plt.tight_layout()
plt.show()
# 理解：
# - proposal_std太小（0.1）：接受率高但移动慢，自相关高，混合差
# - proposal_std适中（1.5）：接受率适中，自相关衰减快，混合好
# - proposal_std太大（10.0）：大部分被拒绝，接受率极低，混合差
四、Gibbs采样（Gibbs Sampling）
4.1 核心概念
Gibbs采样是MCMC的一种特殊情况，它通过从条件分布中轮流采样每个变量来生成样本，不需要手动设定提议分布，且接受率为100%。
算法流程：
对于 kkk 维变量 x=(x1,x2,…,xk)\mathbf{x} = (x_1, x_2, \ldots, x_k)x=(x1,x2,…,xk)：
- 初始化 x(0)=(x1(0),x2(0),…,xk(0))\mathbf{x}^{(0)} = (x_1^{(0)}, x_2^{(0)}, \ldots, x_k^{(0)})x(0)=(x1(0),x2(0),…,xk(0))
- 对于 t=0,1,2,…t = 0, 1, 2, \ldotst=0,1,2,…：
- x1(t+1)∼p(x1∣x2(t),x3(t),…,xk(t))x_1^{(t+1)} \sim p(x_1 | x_2^{(t)}, x_3^{(t)}, \ldots, x_k^{(t)})x1(t+1)∼p(x1∣x2(t),x3(t),…,xk(t))
- x2(t+1)∼p(x2∣x1(t+1),x3(t),…,xk(t))x_2^{(t+1)} \sim p(x_2 | x_1^{(t+1)}, x_3^{(t)}, \ldots, x_k^{(t)})x2(t+1)∼p(x2∣x1(t+1),x3(t),…,xk(t))
- ⋮\vdots⋮
- xk(t+1)∼p(xk∣x1(t+1),x2(t+1),…,xk−1(t+1))x_k^{(t+1)} \sim p(x_k | x_1^{(t+1)}, x_2^{(t+1)}, \ldots, x_{k-1}^{(t+1)})xk(t+1)∼p(xk∣x1(t+1),x2(t+1),…,xk−1(t+1))
关键特点：
- 每一步只改变一个变量，其他变量固定
- 提议分布就是条件分布本身，因此接受率恒为1
- 适合变量之间存在强依赖关系的情况
4.2 代码案例
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
# ============================================================
# 案例1：二元高斯分布的Gibbs采样
# ============================================================
# 目标分布：二元正态分布 N(mu, Sigma)
# mu = [mu_x, mu_y], Sigma = [[sigma_x^2, rho*sigma_x*sigma_y],
# [rho*sigma_x*sigma_y, sigma_y^2]]
#
# 条件分布（已知推导结果）：
# p(x | y) = N(mu_x + rho*(sigma_x/sigma_y)*(y - mu_y), sigma_x^2*(1-rho^2))
# p(y | x) = N(mu_y + rho*(sigma_y/sigma_x)*(x - mu_x), sigma_y^2*(1-rho^2))
# 设置参数
mu_x, mu_y = 0.0, 0.0 # 两个变量的均值
sigma_x, sigma_y = 1.0, 1.5 # 两个变量的标准差
rho = 0.8 # 相关系数（强正相关）
# 计算条件分布参数（预先计算，避免重复计算）
cond_var_x = sigma_x**2 * (1 - rho**2) # p(x|y) 的条件方差
cond_var_y = sigma_y**2 * (1 - rho**2) # p(y|x) 的条件方差
cond_std_x = np.sqrt(cond_var_x) # p(x|y) 的条件标准差
cond_std_y = np.sqrt(cond_var_y) # p(y|x) 的条件标准差
def conditional_x_given_y(y):
"""
从条件分布 p(x | y) 采样
x | y ~ N(mu_x + rho*(sigma_x/sigma_y)*(y - mu_y), sigma_x^2*(1-rho^2))
参数:
y: 固定的y值
返回:
从条件分布中采样的x值
"""
cond_mean = mu_x + rho * (sigma_x / sigma_y) * (y - mu_y) # 条件均值
return np.random.normal(cond_mean, cond_std_x) # 从条件高斯中采样
def conditional_y_given_x(x):
"""
从条件分布 p(y | x) 采样
y | x ~ N(mu_y + rho*(sigma_y/sigma_x)*(x - mu_x), sigma_y^2*(1-rho^2))
参数:
x: 固定的x值
返回:
从条件分布中采样的y值
"""
cond_mean = mu_y + rho * (sigma_y / sigma_x) * (x - mu_x) # 条件均值
return np.random.normal(cond_mean, cond_std_y) # 从条件高斯中采样
def gibbs_sampler_bivariate(n_samples, burn_in=500, init_x=3.0, init_y=-2.0):
"""
二元高斯分布的Gibbs采样器
参数:
n_samples: 需要保留的样本数
burn_in: 预热期迭代数（丢弃前面的样本）
init_x: x的初始值
init_y: y的初始值
返回:
samples: 采样结果数组, shape=(n_samples, 2)
"""
total_iterations = n_samples + burn_in # 总迭代次数
samples = np.zeros((total_iterations, 2)) # 存储所有采样点
x, y = init_x, init_y # 初始化当前状态
samples[0] = [x, y] # 存储初始点
for i in range(1, total_iterations): # 迭代
# 步骤1：固定y，从条件分布 p(x|y) 采样新的x
x = conditional_x_given_y(y) # 从 p(x|y) 采样
# 步骤2：固定x（使用刚刚更新的x），从条件分布 p(y|x) 采样新的y
y = conditional_y_given_x(x) # 从 p(y|x) 采样
samples[i] = [x, y] # 存储新的采样点
# 丢弃burn-in期的样本
return samples[burn_in:] # 返回去除预热期的样本
np.random.seed(42) # 固定随机种子
# 运行Gibbs采样器
n_gibbs = 10000 # 需要保留的样本数
samples_gibbs = gibbs_sampler_bivariate(n_gibbs, burn_in=1000) # 运行采样器
print(f"Gibbs采样完成，样本数: {len(samples_gibbs)}") # 输出样本数
# 可视化
fig, axes = plt.subplots(2, 3, figsize=(18, 10))
# 子图1：采样散点图
axes[0, 0].scatter(samples_gibbs[:, 0], samples_gibbs[:, 1],
s=2, alpha=0.2, c='navy') # 绘制采样点
axes[0, 0].set_xlabel('x')
axes[0, 0].set_ylabel('y')
axes[0, 0].set_title('Gibbs采样散点图')
axes[0, 0].set_aspect('equal') # 等比例坐标轴
axes[0, 0].grid(True, alpha=0.3)
# 子图2：x分量的迹图
axes[0, 1].plot(samples_gibbs[:3000, 0], linewidth=0.3, color='blue')
axes[0, 1].set_xlabel('迭代次数')
axes[0, 1].set_ylabel('x值')
axes[0, 1].set_title('x分量的迹图（Trace Plot）')
axes[0, 1].axhline(y=mu_x, color='red', linestyle='--', label=f'E[x]={mu_x}')
axes[0, 1].legend()
# 子图3：y分量的迹图
axes[0, 2].plot(samples_gibbs[:3000, 1], linewidth=0.3, color='green')
axes[0, 2].set_xlabel('迭代次数')
axes[0, 2].set_ylabel('y值')
axes[0, 2].set_title('y分量的迹图（Trace Plot）')
axes[0, 2].axhline(y=mu_y, color='red', linestyle='--', label=f'E[y]={mu_y}')
axes[0, 2].legend()
# 子图4：x的边际分布
axes[1, 0].hist(samples_gibbs[:, 0], bins=80, density=True, alpha=0.7, color='steelblue',
label='Gibbs采样') # 采样直方图
x_range = np.linspace(-4, 4, 200)
axes[1, 0].plot(x_range, stats.norm.pdf(x_range, mu_x, sigma_x),
'r-', linewidth=2, label='真实边际分布') # 真实边际分布
axes[1, 0].set_title('x的边际分布')
axes[1, 0].legend()
# 子图5：y的边际分布
axes[1, 1].hist(samples_gibbs[:, 1], bins=80, density=True, alpha=0.7, color='coral',
label='Gibbs采样')
y_range = np.linspace(-5, 5, 200)
axes[1, 1].plot(y_range, stats.norm.pdf(y_range, mu_y, sigma_y),
'r-', linewidth=2, label='真实边际分布')
axes[1, 1].set_title('y的边际分布')
axes[1, 1].legend()
# 子图6：Gibbs采样路径图（前100步，展示锯齿形路径）
path = gibbs_sampler_bivariate(100, burn_in=0, init_x=3.0, init_y=-2.0)
axes[1, 2].plot(path[:, 0], path[:, 1], 'o-', markersize=3,
linewidth=0.5, alpha=0.7, color='purple') # 绘制采样路径
axes[1, 2].set_xlabel('x')
axes[1, 2].set_ylabel('y')
axes[1, 2].set_title('Gibbs采样路径（前100步）')
axes[1, 2].set_aspect('equal')
axes[1, 2].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# ============================================================
# 案例2：Gibbs采样用于贝叶斯线性回归
# ============================================================
# 模型：y = w0 + w1*x + noise, noise ~ N(0, sigma^2)
# 先验：w0, w1 ~ N(0, tau^2)
# 使用Gibbs采样对后验 p(w0, w1, sigma^2 | data) 进行采样
np.random.seed(42)
# 生成模拟数据
true_w0 = 2.0 # 真实截距
true_w1 = 0.5 # 真实斜率
true_sigma = 1.0 # 真实噪声标准差
n_data = 50 # 数据点数
x_data = np.random.uniform(-3, 3, n_data) # 生成自变量数据
y_data = true_w0 + true_w1 * x_data + np.random.normal(0, true_sigma, n_data) # 生成因变量
# 超参数
tau2 = 10.0 # 权重先验的方差（弱先验）
def gibbs_linear_regression(x_data, y_data, n_iter=5000, burn_in=1000):
"""
贝叶斯线性回归的Gibbs采样
条件分布推导：
- p(w | sigma^2, data) 是多元高斯
- p(sigma^2 | w, data) 是逆伽马分布
参数:
x_data: 自变量数组
y_data: 因变量数组
n_iter: 总迭代次数
burn_in: 预热期
返回:
w_samples: 权重样本 (n_iter - burn_in, 2)
sigma_samples: 噪声标准差样本 (n_iter - burn_in,)
"""
n = len(x_data) # 数据点数
# 构造设计矩阵 X = [1, x]
X = np.column_stack([np.ones(n), x_data]) # 设计矩阵，第一列全1
# 初始化参数
w = np.array([0.0, 0.0]) # 权重初始值 [w0, w1]
sigma2 = 1.0 # 噪声方差初始值
w_all = [] # 存储权重采样
sigma2_all = [] # 存储方差采样
for iteration in range(n_iter): # 迭代采样
# === 步骤1：采样 w | sigma^2, data ===
# 后验 w | sigma^2, data ~ N(mu_w, Sigma_w)
# Sigma_w = (X^T X / sigma^2 + I/tau^2)^{-1}
# mu_w = Sigma_w * X^T y / sigma^2
XtX = X.T @ X # X^T * X 矩阵乘法
Xty = X.T @ y_data # X^T * y 向量乘法
# 后验协方差矩阵
Sigma_w = np.linalg.inv(XtX / sigma2 + np.eye(2) / tau2) # 2x2协方差矩阵
# 后验均值
mu_w = Sigma_w @ (Xty / sigma2) # 2维均值向量
# 从多元高斯分布采样权重
w = np.random.multivariate_normal(mu_w, Sigma_w) # 采样 w ~ N(mu_w, Sigma_w)
# === 步骤2：采样 sigma^2 | w, data ===
# 后验 sigma^2 | w, data ~ InvGamma(a_n, b_n)
# a_n = a_0 + n/2
# b_n = b_0 + 0.5 * sum((y - Xw)^2)
residuals = y_data - X @ w # 计算残差 y - Xw
ssr = np.sum(residuals**2) # 残差平方和
a_n = 1.0 + n / 2 # 逆伽马分布的形状参数
b_n = 1.0 + ssr / 2 # 逆伽马分布的尺度参数
# 从逆伽马分布采样（通过Gamma分布的倒数实现）
sigma2 = 1.0 / np.random.gamma(a_n, 1.0 / b_n) # InvGamma = 1/Gamma
# 存储采样结果
if iteration >= burn_in: # 超过burn-in期后开始记录
w_all.append(w.copy()) # 存储权重样本
sigma2_all.append(sigma2) # 存储方差样本
return np.array(w_all), np.sqrt(np.array(sigma2_all)) # 返回权重和标准差样本
# 运行Gibbs采样
w_samples, sigma_samples = gibbs_linear_regression(x_data, y_data)
print(f"后验估计 w0: {np.mean(w_samples[:, 0]):.3f} ± {np.std(w_samples[:, 0]):.3f}")
print(f"后验估计 w1: {np.mean(w_samples[:, 1]):.3f} ± {np.std(w_samples[:, 1]):.3f}")
print(f"后验估计 sigma: {np.mean(sigma_samples):.3f} ± {np.std(sigma_samples):.3f}")
print(f"真实值: w0={true_w0}, w1={true_w1}, sigma={true_sigma}")
# 可视化贝叶斯线性回归结果
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
# 子图1：后验采样的回归线
x_plot = np.linspace(-3, 3, 100) # 绘图x范围
for i in range(0, len(w_samples), 50): # 每隔50条画一条回归线
axes[0].plot(x_plot, w_samples[i, 0] + w_samples[i, 1] * x_plot,
'b-', alpha=0.05, linewidth=1) # 后验采样的回归线
# 数据点
axes[0].scatter(x_data, y_data, c='red', s=20, zorder=5, label='数据点')
# 真实回归线
axes[0].plot(x_plot, true_w0 + true_w1 * x_plot, 'g--',
linewidth=2, label='真实回归线')
# 后验均值回归线
axes[0].plot(x_plot, np.mean(w_samples[:, 0]) + np.mean(w_samples[:, 1]) * x_plot,
'r-', linewidth=2, label='后验均值回归线')
axes[0].set_xlabel('x')
axes[0].set_ylabel('y')
axes[0].set_title('贝叶斯线性回归（Gibbs采样）')
axes[0].legend()
# 子图2：w0和w1的后验分布
axes[1].scatter(w_samples[:, 0], w_samples[:, 1], s=1, alpha=0.2, c='navy')
axes[1].axvline(x=true_w0, color='red', linestyle='--', label=f'真实w0={true_w0}')
axes[1].axhline(y=true_w1, color='green', linestyle='--', label=f'真实w1={true_w1}')
axes[1].set_xlabel('w0 (截距)')
axes[1].set_ylabel('w1 (斜率)')
axes[1].set_title('权重后验分布')
axes[1].legend()
# 子图3：sigma的后验分布
axes[2].hist(sigma_samples, bins=60, density=True, alpha=0.7, color='orange')
axes[2].axvline(x=true_sigma, color='red', linestyle='--',
linewidth=2, label=f'真实σ={true_sigma}')
axes[2].set_xlabel('σ')
axes[2].set_ylabel('概率密度')
axes[2].set_title('噪声标准差后验分布')
axes[2].legend()
plt.tight_layout()
plt.show()
# ============================================================
# 案例3：Gibbs采样用于图像去噪（Ising模型）
# ============================================================
def gibbs_image_denoise(noisy_image, n_iter=20, beta=2.0, eta=1.0):
"""
基于Ising模型的Gibbs采样图像去噪
模型：
- 隐变量 x_{i,j} ∈ {-1, +1}（干净图像像素）
- 观测 y_{i,j} ∈ {-1, +1}（含噪图像像素）
- 能量函数: E(x,y) = -beta * sum(x_i * x_j) - eta * sum(x_i * y_i)
- 第一项：相邻像素的平滑先验（鼓励相邻像素取相同值）
- 第二项：观测似然（鼓励与观测一致）
参数:
noisy_image: 含噪二值图像（值为-1或+1）
n_iter: 每个像素的Gibbs迭代次数
beta: 先验强度参数（越大越平滑）
eta: 似然强度参数（越大约束越强）
返回:
denoised: 去噪后的图像
"""
rows, cols = noisy_image.shape # 获取图像尺寸
x = noisy_image.copy() # 初始化隐变量为含噪图像
for iteration in range(n_iter): # 迭代去噪
for i in range(rows): # 遍历每一行
for j in range(cols): # 遍历每一列
# 计算四邻域的和
neighbor_sum = 0 # 初始化邻域和
if i > 0: # 上方邻居
neighbor_sum += x[i-1, j]
if i < rows - 1: # 下方邻居
neighbor_sum += x[i+1, j]
if j > 0: # 左方邻居
neighbor_sum += x[i, j-1]
if j < cols - 1: # 右方邻居
neighbor_sum += x[i, j+1]
# 计算 p(x_{i,j} = +1 | rest) 的未归一化对数概率
# log p(x=+1) = beta * neighbor_sum + eta * y_{i,j}
# log p(x=-1) = -beta * neighbor_sum - eta * y_{i,j}
log_p_pos = beta * neighbor_sum + eta * noisy_image[i, j] # x=+1的对数概率
log_p_neg = -beta * neighbor_sum - eta * noisy_image[i, j] # x=-1的对数概率
# 使用softmax计算概率（数值稳定版本）
max_log_p = max(log_p_pos, log_p_neg) # 取最大值防止溢出
p_pos = np.exp(log_p_pos - max_log_p) # 计算 x=+1 的（未归一化）概率
p_neg = np.exp(log_p_neg - max_log_p) # 计算 x=-1 的（未归一化）概率
prob_pos = p_pos / (p_pos + p_neg) # 归一化得到概率
# 从伯努利分布采样
if np.random.random() < prob_pos: # 以 prob_pos 的概率
x[i, j] = 1 # 设置为 +1
else: # 以 1 - prob_pos 的概率
x[i, j] = -1 # 设置为 -1
if (iteration + 1) % 5 == 0: # 每5次迭代打印进度
print(f" Gibbs去噪: 迭代 {iteration+1}/{n_iter}")
return x # 返回去噪后的图像
# 生成测试图像
np.random.seed(42)
img_size = 30 # 图像尺寸 30x30
# 创建简单的二值测试图像（条纹模式）
clean_image = np.ones((img_size, img_size)) # 全白背景
clean_image[:, 10:15] = -1 # 中间竖条为黑色
clean_image[10:15, :] = -1 # 中间横条为黑色
# 形成十字形图案
# 添加噪声：翻转约20%的像素
noise_prob = 0.2 # 噪声翻转概率
noise_mask = np.random.random((img_size, img_size)) < noise_prob # 生成噪声掩码
noisy_image = clean_image.copy() # 复制干净图像
noisy_image[noise_mask] *= -1 # 翻转被选中的像素
# 运行Gibbs采样去噪
print("开始Gibbs采样去噪...")
denoised_image = gibbs_image_denoise(
noisy_image, # 输入含噪图像
n_iter=10, # 迭代次数
beta=1.5, # 先验强度
eta=1.0 # 似然强度
)
# 计算去噪准确率
noise_rate_before = np.mean(noisy_image != clean_image) # 去噪前的错误率
noise_rate_after = np.mean(denoised_image != clean_image) # 去噪后的错误率
print(f"去噪前错误率: {noise_rate_before:.4f}")
print(f"去噪后错误率: {noise_rate_after:.4f}")
# 可视化
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
axes[0].imshow(clean_image, cmap='gray', vmin=-1, vmax=1) # 显示干净图像
axes[0].set_title('原始干净图像')
axes[0].axis('off')
axes[1].imshow(noisy_image, cmap='gray', vmin=-1, vmax=1) # 显示含噪图像
axes[1].set_title(f'含噪图像 (错误率={noise_rate_before:.2%})')
axes[1].axis('off')
axes[2].imshow(denoised_image, cmap='gray', vmin=-1, vmax=1) # 显示去噪结果
axes[2].set_title(f'Gibbs去噪结果 (错误率={noise_rate_after:.2%})')
axes[2].axis('off')
plt.tight_layout()
plt.show()
五、不同峰值之间的混合挑战（Mixing Challenge）
5.1 核心概念
混合（Mixing）问题 是MCMC方法面临的最核心挑战之一。当目标分布具有多个分离的峰值（模式）时，马尔可夫链可能会长时间被困在某一个峰值附近，难以在不同峰值之间转移，导致：
- 采样偏差：样本主要集中在初始峰值附近，不能正确反映完整的目标分布
- 自相关性高：连续样本高度相关，有效样本量远低于实际样本量
- 收敛缓慢：需要极长的链才能近似正确的分布
根本原因：
- Metropolis类算法使用局部提议（如高斯随机游走），从一个峰值到另一个峰值的路径可能经过概率极低的区域
- 提议点大概率落在低概率区域而被拒绝
- 即使峰值之间的距离不远，若中间有低概率的"谷底"，链也难以穿越
5.2 常见解决方案
| 方法 | 核心思想 |
|---|---|
| 退火重要采样（AIS） | 通过一系列逐渐变化的分布平滑过渡 |
| 并行回火（Parallel Tempering） | 同时运行不同"温度"的链，高温链混合好，通过交换促进混合 |
| 汉密尔顿蒙特卡罗（HMC） | 利用梯度信息引导采样，更有效地穿越低概率区域 |
| 混合提议 | 结合局部和全局提议 |
5.3 代码案例
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
# ============================================================
# 案例1：展示混合问题 —— 双峰分布的标准MH采样
# ============================================================
def bimodal_pdf(x):
"""
计算双峰分布的概率密度
p(x) = 0.4 * N(-4, 1) + 0.6 * N(4, 1.2)
参数:
x: 输入值
返回:
概率密度值
"""
comp1 = 0.4 * np.exp(-0.5 * (x + 4)**2) / np.sqrt(2 * np.pi) # 左侧峰值
comp2 = 0.6 * np.exp(-0.5 * (x - 4)**2 / 1.2**2) / (1.2 * np.sqrt(2 * np.pi)) # 右侧峰值
return comp1 + comp2 # 返回混合密度
def standard_mh_sampler(n_samples, proposal_std=0.5, init_x=-4.0):
"""
标准Metropolis-Hastings采样器（对称高斯提议）
参数:
n_samples: 采样数
proposal_std: 提议分布标准差
init_x: 初始点
返回:
samples: 采样链
"""
samples = np.zeros(n_samples) # 初始化样本数组
samples[0] = init_x # 设置初始点
n_accepted = 0 # 接受计数
for i in range(1, n_samples): # 迭代
current = samples[i-1] # 当前状态
# 对称高斯提议
proposal = current + np.random.normal(0, proposal_std) # 候选点
# MH比率
p_curr = bimodal_pdf(current) # 当前点的概率密度
p_prop = bimodal_pdf(proposal) # 候选点的概率密度
alpha = min(1.0, p_prop / (p_curr + 1e-300)) # 接受概率（避免除零）
if np.random.random() < alpha: # 接受判断
samples[i] = proposal # 接受候选
n_accepted += 1
else:
samples[i] = current # 保持当前状态
return samples, n_accepted / (n_samples - 1) # 返回样本和接受率
np.random.seed(42) # 固定随机种子
# 标准MH：从左侧峰值开始，proposal_std较小
samples_mh, acc_mh = standard_mh_sampler(
n_samples=50000, # 样本数
proposal_std=0.8, # 提议标准差
init_x=-4.0 # 从左侧峰值开始
)
print(f"标准MH 接受率: {acc_mh:.4f}") # 输出接受率
# 可视化混合问题
fig, axes = plt.subplots(2, 2, figsize=(16, 12))
x_range = np.linspace(-10, 10, 1000) # 绘图范围
true_pdf_vals = bimodal_pdf(x_range) # 真实PDF值
# 子图1：迹图 —— 链被困在左侧峰值
axes[0, 0].plot(samples_mh[:10000], linewidth=0.3, color='blue')
axes[0, 0].axhline(y=-4, color='red', linestyle='--', alpha=0.5, label='峰值1位置')
axes[0, 0].axhline(y=4, color='green', linestyle='--', alpha=0.5, label='峰值2位置')
axes[0, 0].set_xlabel('迭代次数')
axes[0, 0].set_ylabel('x值')
axes[0, 0].set_title('标准MH迹图 —— 链被困在单一峰值')
axes[0, 0].legend()
# 子图2：采样直方图 vs 真实分布 —— 严重偏差
axes[0, 1].hist(samples_mh, bins=200, density=True, alpha=0.7, color='steelblue',
label='标准MH采样结果')
axes[0, 1].plot(x_range, true_pdf_vals, 'r-', linewidth=2, label='真实目标分布')
axes[0, 1].set_xlabel('x')
axes[0, 1].set_ylabel('概率密度')
axes[0, 1].set_title('标准MH采样偏差 —— 几乎未探索第二个峰值')
axes[0, 1].legend()
# ============================================================
# 案例2：解决方案 —— 并行回火 (Parallel Tempering)
# ============================================================
def parallel_tempering_sampler(n_samples, n_temps=10, beta_max=1.0, proposal_std=1.0):
"""
并行回火采样器
核心思想：
- 同时运行多条链，每条链使用不同的"温度"
- 高温链（beta小）的目标分布更平坦，更容易在峰值间移动
- 低温链（beta大，特别是beta=1）精确采样目标分布
- 通过相邻温度链之间的交换操作促进信息传递
目标：在温度T下，采样 p(x)^{beta} ∝ exp(beta * log p(x))
其中 beta = 1/T 是逆温度
参数:
n_samples: 每条链的采样数
n_temps: 温度层数
beta_max: 最大逆温度（=1对应目标分布）
proposal_std: 提议标准差
返回:
all_chains: 所有温度链的采样结果
"""
# 设置温度序列（几何序列）
betas = np.linspace(0.1, beta_max, n_temps) # 逆温度序列 [0.1, ..., 1.0]
# 初始化所有链
chains = np.full(n_temps, -4.0) # 所有链都从左侧峰值开始
all_chains = np.zeros((n_samples, n_temps)) # 存储所有链的所有采样
all_chains[0] = chains.copy() # 存储初始状态
n_swaps = 0 # 交换计数
n_swap_attempts = 0 # 交换尝试计数
for i in range(1, n_samples): # 迭代采样
# === 步骤1：每条链独立进行一步MH采样 ===
for t in range(n_temps): # 遍历每个温度
current = chains[t] # 当前链的状态
proposal = current + np.random.normal(0, proposal_std) # 提议候选点
# 在温度beta_t下的接受概率
# alpha = min(1, [p(x')^beta_t] / [p(x)^beta_t])
# = min(1, p(x')/p(x))^beta_t
p_curr = bimodal_pdf(current) + 1e-300 # 当前点的密度
p_prop = bimodal_pdf(proposal) + 1e-300 # 候选点的密度
log_ratio = np.log(p_prop) - np.log(p_curr) # 对数比率
alpha = min(1.0, np.exp(betas[t] * log_ratio)) # 温度调制的接受概率
if np.random.random() < alpha: # 接受判断
chains[t] = proposal # 更新链状态
# === 步骤2：相邻温度链之间的交换操作 ===
# 随机选择一对相邻温度（奇偶交替以保证各态历经性）
if i % 2 == 0: # 偶数步：交换 (0,1), (2,3), ...
pairs = [(k, k+1) for k in range(0, n_temps-1, 2)]
else: # 奇数步：交换 (1,2), (3,4), ...
pairs = [(k, k+1) for k in range(1, n_temps-1, 2)]
for t1, t2 in pairs: # 遍历每对相邻温度
n_swap_attempts += 1 # 交换尝试计数
# 计算交换的接受概率
# alpha_swap = min(1, [p(x1)^beta2 * p(x2)^beta1] / [p(x1)^beta1 * p(x2)^beta2])
# = min(1, p(x1)^(beta2-beta1) * p(x2)^(beta1-beta2))
x1, x2 = chains[t1], chains[t2] # 获取两条链的状态
log_swap_ratio = (betas[t2] - betas[t1]) * (
np.log(bimodal_pdf(x1) + 1e-300) - np.log(bimodal_pdf(x2) + 1e-300)
) # 计算交换的对数接受比
alpha_swap = min(1.0, np.exp(log_swap_ratio)) # 交换接受概率
if np.random.random() < alpha_swap: # 若接受交换
chains[t1], chains[t2] = x2, x1 # 交换两条链的状态
n_swaps += 1 # 交换计数+1
all_chains[i] = chains.copy() # 存储当前所有链的状态
swap_rate = n_swaps / n_swap_attempts if n_swap_attempts > 0 else 0 # 交换率
print(f"并行回火 - 温度层数: {n_temps}, 交换率: {swap_rate:.4f}")
return all_chains, betas # 返回所有链和温度序列
# 运行并行回火
np.random.seed(42)
n_pt_samples = 50000
all_chains_pt, betas_pt = parallel_tempering_sampler(
n_samples=n_pt_samples,
n_temps=10, # 10个温度层
proposal_std=1.0 # 提议标准差
)
# 只取最低温（beta=1.0，即目标温度）的链
samples_pt = all_chains_pt[:, -1] # 最后一列对应beta=1的链
# 子图3：并行回火的迹图
axes[1, 0].plot(samples_pt[:10000], linewidth=0.3, color='purple')
axes[1, 0].axhline(y=-4, color='red', linestyle='--', alpha=0.5, label='峰值1')
axes[1, 0].axhline(y=4, color='green', linestyle='--', alpha=0.5, label='峰值2')
axes[1, 0].set_xlabel('迭代次数')
axes[1, 0].set_ylabel('x值')
axes[1, 0].set_title('并行回火迹图 —— 能在峰值间自由移动')
axes[1, 0].legend()
# 子图4：并行回火采样直方图 vs 真实分布
axes[1, 1].hist(samples_pt, bins=200, density=True, alpha=0.7, color='mediumpurple',
label='并行回火采样结果')
axes[1, 1].plot(x_range, true_pdf_vals, 'r-', linewidth=2, label='真实目标分布')
axes[1, 1].set_xlabel('x')
axes[1, 1].set_ylabel('概率密度')
axes[1, 1].set_title('并行回火采样 —— 正确覆盖两个峰值')
axes[1, 1].legend()
plt.tight_layout()
plt.show()
# ============================================================
# 案例3：不同温度链的可视化
# ============================================================
fig, axes = plt.subplots(3, 3, figsize=(18, 12))
temp_indices = [0, 4, 9] # 选择3个代表性温度
sample_range = slice(0, 5000) # 显示前5000个样本
for row, t_idx in enumerate(temp_indices): # 遍历选定的温度
chain = all_chains_pt[sample_range, t_idx] # 获取该温度的链
beta_val = betas_pt[t_idx] # 获取该温度的beta值
# 迹图
axes[row, 0].plot(chain, linewidth=0.3, color=['orange', 'blue', 'red'][row])
axes[row, 0].set_ylabel(f'β={beta_val:.2f}')
axes[row, 0].axhline(y=-4, color='gray', linestyle='--', alpha=0.3)
axes[row, 0].axhline(y=4, color='gray', linestyle='--', alpha=0.3)
if row == 0:
axes[row, 0].set_title('迹图')
# 直方图
axes[row, 1].hist(chain, bins=100, density=True, alpha=0.7,
color=['orange', 'blue', 'red'][row])
axes[row, 1].plot(x_range, true_pdf_vals * betas_pt[t_idx], 'r--',
linewidth=1, alpha=0.5)
if row == 0:
axes[row, 1].set_title('采样直方图')
# 自相关函数
acf_vals = compute_autocorrelation(chain, 200) if len(chain) > 200 else np.zeros(200)
axes[row, 2].bar(range(100), acf_vals[:100], color=['orange', 'blue', 'red'][row], alpha=0.7)
axes[row, 2].axhline(y=0, color='black', linestyle='-', linewidth=0.5)
if row == 0:
axes[row, 2].set_title('自相关函数')
axes[row, 0].set_xlabel('迭代')
axes[row, 1].set_xlabel('x')
axes[row, 2].set_xlabel('滞后')
plt.suptitle('并行回火中不同温度链的对比分析', fontsize=16, y=1.02)
plt.tight_layout()
plt.show()
# ============================================================
# 案例4：混合挑战的定量分析 —— 有效样本数比较
# ============================================================
def effective_sample_size(chain):
"""
计算有效样本数（ESS）
ESS = N / (1 + 2 * sum_{k=1}^{K} rho_k)
其中 rho_k 是滞后k的自相关系数，K取到自相关首次变为负值
参数:
chain: 一维采样链
返回:
ess: 有效样本数
"""
n = len(chain) # 链长度
mean = np.mean(chain) # 均值
var = np.var(chain) # 方差
if var < 1e-10: # 链几乎没有变化
return 1.0 # ESS极小
# 计算自相关函数直到第一次变为负
autocorr_sum = 0 # 自相关累积和
for lag in range(1, min(n // 2, 5000)): # 遍历滞后值
cov = np.mean((chain[:n-lag] - mean) * (chain[lag:] - mean)) # 自协方差
rho = cov / var # 自相关系数
if rho < 0: # 若自相关变为负
break # 停止累加
autocorr_sum += rho # 累加正的自相关
ess = n / (1 + 2 * autocorr_sum) # 计算ESS
return ess
# 比较不同方法的ESS
ess_mh = effective_sample_size(samples_mh) # 标准MH的ESS
ess_pt = effective_sample_size(samples_pt) # 并行回火的ESS
print("\n" + "=" * 60)
print("混合效果定量比较")
print("=" * 60)
print(f"标准MH: 总样本={len(samples_mh)}, ESS={ess_mh:.0f}, "
f"ESS/N={ess_mh/len(samples_mh):.4f}")
print(f"并行回火: 总样本={len(samples_pt)}, ESS={ess_pt:.0f}, "
f"ESS/N={ess_pt/len(samples_pt):.4f}")
# 比较在两个峰值附近的采样比例
left_peak_ratio_mh = np.mean(samples_mh < 0) # 标准MH中左侧峰值的采样比例
left_peak_ratio_pt = np.mean(samples_pt < 0) # 并行回火中左侧峰值的采样比例
# 真实比例
# p(x<0) = 0.4 * Phi(0; -4, 1) 积分...近似计算
true_left_ratio, _ = quad(lambda x: bimodal_pdf(x), -20, 0) # 数值积分计算左侧真实概率
true_left_ratio = true_left_ratio / quad(bimodal_pdf, -20, 20)[0] # 归一化
print(f"\n左侧峰值 (x<0) 的采样比例:")
print(f" 标准MH: {left_peak_ratio_mh:.4f}")
print(f" 并行回火: {left_peak_ratio_pt:.4f}")
print(f" 真实比例: {true_left_ratio:.4f}")
# ============================================================
# 案例5：三维多峰分布的混合挑战可视化
# ============================================================
# 使用简单的MH + 大跳跃提议来处理多峰问题
def multimodal_3d_pdf(x):
"""
三维三峰分布
三个峰值分别位于 (0,0,0), (5,5,5), (-5,5,0)
参数:
x: 形状为(3,)的数组
返回:
概率密度值
"""
peaks = [np.array([0, 0, 0]), np.array([5, 5, 5]), np.array([-5, 5, 0])] # 峰值位置
weights = [0.3, 0.4, 0.3] # 各峰权重
stds = [1.0, 1.2, 0.8] # 各峰标准差
density = 0 # 初始化密度
for mu, w, s in zip(peaks, weights, stds): # 遍历每个峰
dist = np.linalg.norm(x - mu) # 计算到峰值的距离
density += w * np.exp(-0.5 * (dist / s)**2) / (s**3 * (2*np.pi)**1.5)
return density # 返回混合密度
def mh_with_jump_proposals(n_samples, local_std=0.5, jump_prob=0.05, jump_std=5.0):
"""
带有大跳跃提议的MH采样器
以 (1-jump_prob) 的概率使用局部提议，以 jump_prob 的概率使用全局跳跃提议
参数:
n_samples: 采样数
local_std: 局部提议标准差
jump_prob: 大跳跃概率
jump_std: 大跳跃提议标准差
返回:
samples: 采样结果
accept_rate: 接受率
"""
samples = np.zeros((n_samples, 3)) # 初始化样本数组
samples[0] = np.array([0.0, 0.0, 0.0]) # 从原点开始
n_accepted = 0 # 接受计数
for i in range(1, n_samples): # 迭代
current = samples[i-1] # 当前状态
p_curr = multimodal_3d_pdf(current) + 1e-300 # 当前点密度
if np.random.random() < jump_prob: # 以小概率进行大跳跃
proposal = current + np.random.normal(0, jump_std, 3) # 大跳跃提议
else: # 大部分时间使用局部提议
proposal = current + np.random.normal(0, local_std, 3) # 局部随机游走
p_prop = multimodal_3d_pdf(proposal) + 1e-300 # 候选点密度
# 考虑非对称提议的修正（混合提议的MH比率需要特殊处理）
# 这里简化处理，因为局部和跳跃都是对称的
alpha = min(1.0, p_prop / p_curr) # 接受概率
if np.random.random() < alpha: # 接受判断
samples[i] = proposal # 接受
n_accepted += 1
else:
samples[i] = current # 拒绝
return samples, n_accepted / (n_samples - 1) # 返回样本和接受率
np.random.seed(42)
# 标准局部MH
samples_local, ar_local = mh_with_jump_proposals(
n_samples=30000, jump_prob=0.0 # 不使用跳跃
)
# 混合提议MH
samples_mixed, ar_mixed = mh_with_jump_proposals(
n_samples=30000, jump_prob=0.05 # 5%概率大跳跃
)
print(f"\n三维多峰采样:")
print(f" 纯局部MH 接受率: {ar_local:.4f}")
print(f" 混合提议MH 接受率: {ar_mixed:.4f}")
# 可视化（3D散点图投影到2D）
fig, axes = plt.subplots(1, 2, figsize=(16, 6))
# 纯局部MH
axes[0].scatter(samples_local[:, 0], samples_local[:, 1], s=1, alpha=0.1, c='blue')
axes[0].scatter([0, 5, -5], [0, 5, 5], c='red', s=100, marker='*', zorder=5,
label='真实峰值位置')
axes[0].set_xlabel('x')
axes[0].set_ylabel('y')
axes[0].set_title(f'纯局部MH (接受率={ar_local:.2f}) —— 仅探索初始峰值附近')
axes[0].legend()
axes[0].set_xlim(-10, 10)
axes[0].set_ylim(-5, 12)
# 混合提议MH
axes[1].scatter(samples_mixed[:, 0], samples_mixed[:, 1], s=1, alpha=0.1, c='purple')
axes[1].scatter([0, 5, -5], [0, 5, 5], c='red', s=100, marker='*', zorder=5,
label='真实峰值位置')
axes[1].set_xlabel('x')
axes[1].set_ylabel('y')
axes[1].set_title(f'混合提议MH (接受率={ar_mixed:.2f}) —— 成功探索所有峰值')
axes[1].legend()
axes[1].set_xlim(-10, 10)
axes[1].set_ylim(-5, 12)
plt.tight_layout()
plt.show()
六、总结对比表
+------------------+-------------------+-------------------+-------------------+-------------------+
| 方法 | 需要采样分布 | 接受率 | 适用场景 | 混合能力 |
+------------------+-------------------+-------------------+-------------------+-------------------+
| 直接采样 | 需要知道CDF逆 | 100% | 简单分布 | N/A |
| 拒绝采样 | 需要包络函数 | < 100%（依赖M） | 中等复杂分布 | 一般 |
| 重要采样 | 需要提议分布 | 100%（估计而非采样）| 期望估计 | N/A（非MCMC） |
| Metropolis-Hast. | 需要提议分布 | < 100% | 通用 | 取决于提议 |
| Gibbs采样 | 需要条件分布 | 100% | 条件分布已知 | 一般（局部更新） |
| 并行回火 | 需要提议分布 | < 100% + 交换率 | 多峰分布 | 强 |
+------------------+-------------------+-------------------+-------------------+-------------------+
更多推荐
1
0- 0
已为社区贡献49条内容
所有评论(0)