1. 物理学家的酒吧对话：Metropolis算法的诞生
1953年，洛斯阿拉莫斯国家实验室的五位物理学家（Nicholas Metropolis、Arianna Rosenbluth、Marshall Rosenbluth、Augusta Teller和Edward Teller）在当地的酒吧里讨论一个棘手的问题：如何模拟液体分子的复杂运动？ 传统的蒙特卡罗方法在高维空间中效率极低，无法有效采样。
在几杯鸡尾酒的启发下，他们提出了一个革命性的想法：构造一个马尔可夫链，使其平稳分布等于目标分布。这个后来被称为Metropolis算法的突破，成为了统计物理学和计算统计学的里程碑。
"我们意识到，"Metropolis后来回忆道，"不需要直接从复杂分布中采样，而是可以构造一个随机游走，使得游走的路径最终会服从目标分布。这就像醉汉走路，虽然每一步都随机，但长期来看会在某些区域停留更久。"
2. MCMC的数学心脏：细致平衡条件
MCMC方法的核心是满足细致平衡条件（Detailed Balance）：
π(x)T(x→x′)=π(x′)T(x′→x)
其中π是目标分布，T是转移核。这个条件保证了马尔可夫链的平稳分布就是目标分布。
让我们用Python实现经典的Metropolis-Hastings算法：
import numpy as np
import matplotlib.pyplot as plt
def metropolis_hastings(target_density, proposal, initial, iterations=10000, burn_in=1000):
"""
target_density: 目标分布(非归一化)
proposal: 建议分布函数 q(x'|x)
initial: 初始状态
iterations: 迭代次数
burn_in: 预烧期
"""
samples = []
current = initial
accepted = 0
for i in range(iterations):
# 生成候选样本
candidate = proposal(current)
# 计算接受概率
acceptance_ratio = (target_density(candidate) / target_density(current)) * \
(proposal(current, candidate) / proposal(candidate, current))
acceptance_prob = min(1, acceptance_ratio)
# 决定是否接受
if np.random.rand() < acceptance_prob:
current = candidate
accepted += 1
# 记录样本(跳过预烧期)
if i >= burn_in:
samples.append(current)
acceptance_rate = accepted / iterations
return np.array(samples), acceptance_rate
# 示例：采样双峰分布
def target_density(x):
"""混合高斯分布"""
return 0.3 * np.exp(-0.2 * (x - 3)**2) + 0.7 * np.exp(-0.2 * (x + 3)**2)
def proposal(x, sigma=2):
"""对称建议分布：正态分布"""
return np.random.normal(x, sigma)
# 运行采样
samples, acc_rate = metropolis_hastings(target_density, proposal, 0, 15000)
print(f"接受率: {acc_rate:.2%}")
# 可视化
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.plot(samples)
plt.title('MCMC采样路径')
plt.subplot(1, 2, 2)
plt.hist(samples, bins=50, density=True, alpha=0.6)
x = np.linspace(-6, 6, 100)
plt.plot(x, target_density(x) / np.trapz(target_density(x), x), 'r', linewidth=2)
plt.title('采样分布 vs 真实分布')
plt.show()
3. Gibbs抽样：高维问题的优雅解法
Gibbs抽样是MCMC的一种特殊情形，由统计学家Stuart Geman和Donald Geman于1984年提出，用于图像处理。它特别适合高维分布的采样，通过轮流更新每个维度来实现：
def gibbs_sampler(conditional_dists, initial, iterations=10000, burn_in=1000):
"""
conditional_dists: 各变量的条件分布函数列表
initial: 初始状态
iterations: 迭代次数
burn_in: 预烧期
"""
dim = len(initial)
samples = []
current = initial.copy()
for i in range(iterations):
# 依次更新每个维度
for d in range(dim):
current[d] = conditional_dists[d](current)
if i >= burn_in:
samples.append(current.copy())
return np.array(samples)
# 示例：二元正态分布
def cond_x1(x):
"""给定x2，x1的条件分布 ~ N(0.5*x[1], sqrt(0.75))"""
return np.random.normal(0.5 * x[1], np.sqrt(0.75))
def cond_x2(x):
"""给定x1，x2的条件分布 ~ N(0.5*x[0], sqrt(0.75))"""
return np.random.normal(0.5 * x[0], np.sqrt(0.75))
# 运行Gibbs抽样
samples = gibbs_sampler([cond_x1, cond_x2], [0, 0], 5000)
# 可视化
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.scatter(samples[:, 0], samples[:, 1], alpha=0.5)
plt.title('采样点分布')
plt.subplot(1, 2, 2)
plt.plot(samples[:200, 0], label='x1')
plt.plot(samples[:200, 1], label='x2')
plt.title('前200次迭代的采样轨迹')
plt.legend()
plt.show()
4. MCMC在贝叶斯统计中的革命
贝叶斯统计的核心是计算后验分布：
p(θ∣D)=p(D)p(D∣θ)p(θ)
其中p(D)（边缘似然）通常难以计算。MCMC通过直接从后验分布中采样，绕过了这个难题。
4.1 贝叶斯线性回归案例
import numpy as np
from scipy.stats import multivariate_normal, invgamma
def bayesian_linear_regression(X, y, n_samples=10000):
"""
X: 设计矩阵 (n_samples, n_features)
y: 响应变量 (n_samples,)
返回: beta和sigma2的样本
"""
n, p = X.shape
XtX = X.T @ X
# 初始化
beta = np.zeros(p)
sigma2 = 1
samples = []
# 先验参数
mu0 = np.zeros(p)
Sigma0 = np.eye(p)
a0 = 1
b0 = 1
for _ in range(n_samples):
# 1. 采样beta | sigma2, y
Sigma1 = np.linalg.inv(XtX / sigma2 + np.linalg.inv(Sigma0))
mu1 = Sigma1 @ (X.T @ y / sigma2 + np.linalg.inv(Sigma0) @ mu0)
beta = multivariate_normal.rvs(mu1, Sigma1)
# 2. 采样sigma2 | beta, y
residuals = y - X @ beta
a1 = a0 + n / 2
b1 = b0 + (residuals @ residuals) / 2
sigma2 = invgamma.rvs(a1, scale=b1)
samples.append((beta.copy(), sigma2))
return samples
# 生成模拟数据
np.random.seed(42)
n = 100
X = np.column_stack([np.ones(n), np.random.normal(size=(n, 2))])
true_beta = np.array([1, 2, -1])
true_sigma = 1
y = X @ true_beta + np.random.normal(0, true_sigma, n)
# 运行MCMC
samples = bayesian_linear_regression(X, y, 10000)
beta_samples = np.array([s[0] for s in samples])
sigma_samples = np.array([s[1] for s in samples])
# 分析结果
print("参数估计(后验均值):")
print(f"β0: {np.mean(beta_samples[:, 0]):.3f} (真实值: {true_beta[0]})")
print(f"β1: {np.mean(beta_samples[:, 1]):.3f} (真实值: {true_beta[1]})")
print(f"β2: {np.mean(beta_samples[:, 2]):.3f} (真实值: {true_beta[2]})")
print(f"σ²: {np.mean(sigma_samples):.3f} (真实值: {true_sigma**2})")
# 绘制轨迹图
plt.figure(figsize=(12, 8))
for i in range(3):
plt.subplot(2, 2, i+1)
plt.plot(beta_samples[:, i])
plt.axhline(true_beta[i], color='r', linestyle='--')
plt.title(f'β{i}轨迹')
plt.subplot(2, 2, 4)
plt.plot(sigma_samples)
plt.axhline(true_sigma**2, color='r', linestyle='--')
plt.title('σ²轨迹')
plt.tight_layout()
plt.show()
5. 现代变种：Hamiltonian Monte Carlo (HMC)
HMC是MCMC的高级变种，由物理学家Simon Duane等人在1987年提出，利用哈密顿动力学实现更高效的采样：
import numpy as np
from autograd import grad
def hamiltonian_monte_carlo(target_log_prob, initial, step_size, n_steps, iterations=10000):
"""
target_log_prob: 目标分布的对数概率函数
initial: 初始状态
step_size: 积分步长
n_steps: 积分步数
iterations: 采样次数
"""
dim = len(initial)
samples = [initial]
current = initial.copy()
grad_log_prob = grad(target_log_prob)
for _ in range(iterations):
# 1. 随机生成动量
p = np.random.normal(size=dim)
current_p = p.copy()
# 2. 模拟哈密顿动力学
q = current.copy()
q_grad = grad_log_prob(q)
# 蛙跳积分
for _ in range(n_steps):
# 更新动量半步
p += 0.5 * step_size * q_grad
# 更新位置全步
q += step_size * p
# 更新梯度
q_grad_new = grad_log_prob(q)
# 更新动量半步
p += 0.5 * step_size * q_grad_new
q_grad = q_grad_new
# 3. Metropolis接受/拒绝
current_U = -target_log_prob(current)
current_K = 0.5 * np.sum(current_p**2)
proposed_U = -target_log_prob(q)
proposed_K = 0.5 * np.sum(p**2)
if np.random.rand() < np.exp(current_U + current_K - proposed_U - proposed_K):
current = q
samples.append(current.copy())
return np.array(samples)
# 示例：采样二维正态分布
def target_log_prob(x):
"""对数概率密度：二元正态分布"""
return -0.5 * (x[0]**2 + (x[1] - 0.5 * x[0])**2)
# 运行HMC
samples = hamiltonian_monte_carlo(target_log_prob, [0, 0], 0.1, 10, 1000)
# 可视化
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.scatter(samples[:, 0], samples[:, 1], alpha=0.5)
plt.title('HMC采样点')
plt.subplot(1, 2, 2)
plt.plot(samples[:, 0], label='x1')
plt.plot(samples[:, 1], label='x2')
plt.title('采样轨迹')
plt.legend()
plt.show()
6. MCMC在深度学习中的应用
6.1 贝叶斯神经网络
import torch
import torch.nn as nn
import torch.distributions as dist
class BayesianNN(nn.Module):
def __init__(self, input_dim, hidden_dim, output_dim):
super().__init__()
self.fc1 = nn.Linear(input_dim, hidden_dim)
self.fc2 = nn.Linear(hidden_dim, output_dim)
# 先验参数
self.prior_sigma = 1.0
def forward(self, x):
x = torch.relu(self.fc1(x))
return self.fc2(x)
def log_prior(self):
"""计算权重的对数先验"""
log_prior = 0
for param in self.parameters():
log_prior += dist.Normal(0, self.prior_sigma).log_prob(param).sum()
return log_prior
def log_likelihood(self, x, y):
"""计算对数似然"""
pred = self.forward(x)
return dist.Normal(pred, 0.1).log_prob(y).sum()
def log_posterior(self, x, y):
"""计算对数后验(非归一化)"""
return self.log_prior() + self.log_likelihood(x, y)
def mcmc_sampling(model, x, y, n_samples, step_size):
"""使用Metropolis算法采样神经网络权重"""
current_log_prob = model.log_posterior(x, y)
samples = []
for _ in range(n_samples):
# 保存当前状态
current_state = {k: v.clone() for k, v in model.state_dict().items()}
# 生成候选样本 - 随机游走
with torch.no_grad():
for param in model.parameters():
param.add_(torch.randn_like(param) * step_size)
# 计算接受概率
candidate_log_prob = model.log_posterior(x, y)
acceptance = min(0, candidate_log_prob - current_log_prob)
# 决定是否接受
if torch.log(torch.rand(1)) < acceptance:
current_log_prob = candidate_log_prob
else:
model.load_state_dict(current_state)
samples.append({k: v.clone() for k, v in model.state_dict().items()})
return samples
# 示例使用
model = BayesianNN(2, 10, 1)
x = torch.randn(100, 2)
y = torch.sin(x[:, 0]) + torch.cos(x[:, 1]) + 0.1 * torch.randn(100)
samples = mcmc_sampling(model, x, y, n_samples=1000, step_size=0.01)
6.2 变分自编码器(VAE)中的MCMC
import torch
import torch.nn as nn
import torch.nn.functional as F
class VAE(nn.Module):
def __init__(self, input_dim, hidden_dim, latent_dim):
super().__init__()
# 编码器
self.encoder = nn.Sequential(
nn.Linear(input_dim, hidden_dim),
nn.ReLU(),
nn.Linear(hidden_dim, hidden_dim),
nn.ReLU()
)
# 潜在空间参数
self.fc_mu = nn.Linear(hidden_dim, latent_dim)
self.fc_var = nn.Linear(hidden_dim, latent_dim)
# 解码器
self.decoder = nn.Sequential(
nn.Linear(latent_dim, hidden_dim),
nn.ReLU(),
nn.Linear(hidden_dim, hidden_dim),
nn.ReLU(),
nn.Linear(hidden_dim, input_dim),
nn.Sigmoid()
)
def encode(self, x):
h = self.encoder(x)
return self.fc_mu(h), self.fc_var(h)
def reparameterize(self, mu, logvar):
std = torch.exp(0.5 * logvar)
eps = torch.randn_like(std)
return mu + eps * std
def decode(self, z):
return self.decoder(z)
def forward(self, x):
mu, logvar = self.encode(x)
z = self.reparameterize(mu, logvar)
return self.decode(z), mu, logvar
def mcmc_vae(vae, x, n_steps=10, step_size=0.1):
"""在VAE潜在空间中使用MCMC采样"""
mu, logvar = vae.encode(x)
z = vae.reparameterize(mu, logvar)
# 定义目标分布(对数后验)
def log_posterior(z_sample):
recon = vae.decode(z_sample)
recon_loss = F.binary_cross_entropy(recon, x, reduction='sum')
kl_div = -0.5 * torch.sum(1 + logvar - mu.pow(2) - logvar.exp())
return -(recon_loss + kl_div)
# 运行MCMC
for _ in range(n_steps):
# 生成候选样本
z_candidate = z + torch.randn_like(z) * step_size
# 计算接受概率
log_alpha = log_posterior(z_candidate) - log_posterior(z)
if torch.log(torch.rand(1)) < log_alpha:
z = z_candidate
return vae.decode(z)
# 示例使用
vae = VAE(input_dim=784, hidden_dim=400, latent_dim=20)
x_sample = torch.randn(1, 784) # 模拟MNIST图像
recon_mcmc = mcmc_vae(vae, x_sample)
7. 结语：MCMC的过去、现在与未来
从Metropolis等人在洛斯阿拉莫斯酒吧的灵感，到如今支撑着贝叶斯统计、深度学习和科学计算的各个领域，MCMC走过了70余年的辉煌历程。它解决了传统方法无法处理的高维积分、复杂分布采样等难题，成为了现代统计学不可或缺的工具。
"MCMC最令人惊叹的是，"著名统计学家Persi Diaconis评价道，"它将物理学的直觉与数学的严谨完美结合，创造了一种既实用又深刻的方法论。"
未来，随着量子MCMC、神经MCMC等新兴技术的发展，这一领域仍将保持旺盛的生命力。正如Metropolis本人所说："好的算法就像好的理论，它们会不断进化，超越最初设计者的想象。"
扫一扫
5977
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
