19 马尔可夫链蒙特卡洛方法之MCMC方法的基础

在前一篇文章中，我们讨论了贝叶斯分类中的模型评估与改进，强调了如何通过有效的模型选择和验证来提升分类器的表现。这次，我们将介绍马尔可夫链蒙特卡洛（MCMC）方法的基础知识，作为下一篇关于Gibbs采样的铺垫。

MCMC方法简介

马尔可夫链蒙特卡洛（MCMC）是一种利用马尔可夫链来进行概率分布采样的统计方法。它非常适用于那些难以直接从中采样的复杂分布，尤其是在贝叶斯统计中，MCMC方法被广泛应用于后验分布的估计。

MCMC的基本原理

在MCMC中，我们需要构造一个马尔可夫链，使其稳定分布为我们想要的目标分布。该马尔可夫链通过一系列状态转移在状态空间中进行随机游走。在足够的时间后，这个链将收敛到目标分布上。

一个常见的MCMC方法是“Metropolis-Hastings”算法，其核心步骤包括：

1. 初始化：选择一个初始状态$x_0$。
2. 迭代：
   • 从提议分布$q(x’|x_t)$中生成一个候选状态$x’$。
   • 计算接受率：
     $$ \alpha = \min\left(1, \frac{\pi(x’) q(x_t|x’)}{\pi(x_t) q(x’|x_t)}\right) $$
   • 以概率$\alpha$接受状态$x’$，否则保持当前状态$x_t$。
3. 重复迭代直到达到所需的样本量。

案例分析：MCMC在贝叶斯回归中的应用

我们通过一个简单的案例来演示MCMC方法的应用。在贝叶斯线性回归中，我们假设模型为：

$$ y = \beta_0 + \beta_1 x + \epsilon $$

其中，$\epsilon \sim N(0, \sigma^2)$。我们希望从后验分布中抽样，后验分布的形式为：

$$ p(\beta | y, x) \propto p(y | x, \beta) p(\beta) $$

我们选择一个简单的先验分布，比如$\beta \sim N(0, 10^2)$。

使用Python实现MCMC

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

# 生成模拟数据
np.random.seed(42)
x = np.random.rand(100)
true_beta = [2, 3]
y = true_beta[0] + true_beta[1] * x + np.random.normal(0, 0.5, x.shape)

# MCMC实现
def mcmc_bayesian_linear_regression(x, y, iterations=10000):
    n = len(y)
    beta_samples = np.zeros((iterations, 2))
    beta = np.random.randn(2)  # 初始参数

    for i in range(iterations):
        # 生成新样本
        beta_new = beta + np.random.normal(0, 0.5, 2)

        # 计算后验概率比例
        likelihood_current = np.sum(np.log(np.exp(-0.5 * ((y - (beta[0] + beta[1] * x)) / 0.5)**2)))
        likelihood_new = np.sum(np.log(np.exp(-0.5 * ((y - (beta_new[0] + beta_new[1] * x)) / 0.5)**2)))
        
        prior_current = np.exp(-0.5 * (beta[0]**2 / 10**2 + beta[1]**2 / 10**2))
        prior_new = np.exp(-0.5 * (beta_new[0]**2 / 10**2 + beta_new[1]**2 / 10**2))
        
        acceptance_ratio = (likelihood_new * prior_new) / (likelihood_current * prior_current)
        if np.random.rand() < acceptance_ratio:
            beta = beta_new
        
        beta_samples[i] = beta
    
    return beta_samples

# 运行MCMC
beta_samples = mcmc_bayesian_linear_regression(x, y)

# 绘制结果
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.plot(beta_samples[:, 0], label=r'$\beta_0$ samples')
plt.title('Trace plot of $\beta_0$')
plt.subplot(1, 2, 2)
plt.plot(beta_samples[:, 1], label=r'$\beta_1$ samples', color='orange')
plt.title('Trace plot of $\beta_1$')
plt.legend()
plt.show()

# 可视化后验分布
sns.kdeplot(beta_samples[:, 0], label='Posterior of $\\beta_0$', fill=True)
sns.kdeplot(beta_samples[:, 1], label='Posterior of $\\beta_1$', fill=True, color='orange')
plt.title('Posterior distributions')
plt.legend()
plt.show()

在上面的代码中，我们实现了一个简单的贝叶斯线性回归模型的MCMC采样。我们通过迭代生成样本，并绘制了参数$\beta_0$和$\beta_1$的轨迹图和后验分布。这让我们直观地看到MCMC方法如何从复杂的后验分布中提取样本。

结语

MCMC方法为我们提供了一种强大灵活的工具，在复杂的贝叶斯模型中进行参数推断。在下一个讲解中，我们将深入探讨更具体的MCMC变体——Gibbs采样，这是一种在某些条件下更为高效的采样方法。在应用MCMC方法时，理解其原理和实现是至关重要的，这将为您在贝叶斯统计的旅程中打下坚定的基础。