21 Metropolis-Hastings算法

查看大图

贝叶斯学习的重点是把已有判断和新证据合在一起，并明确表达不确定性。阅读时可以按「算法原理 -> 目标分布 -> 提议分布 -> 接受率」建立结构，再回到正文里的代码、案例或指标做验证。

查看大图

读完后，用一个真实小任务复查：输入是什么，处理环节在哪里，输出是否可验收；失败时先查「算法原理」，再查「目标分布」。

在马尔可夫链蒙特卡洛（MCMC）方法中，Metropolis-Hastings算法是一个重要的产生采样方法，广泛应用于贝叶斯学习和统计推断的领域。它是一种构建马尔可夫链的方法，通过从一个指定的概率分布中生成样本，来近似该分布。接下来，我们将详细探讨该算法的原理、步骤和应用案例。

算法原理

Metropolis-Hastings算法是Gibbs采样的推广，适用于更为复杂的联合分布。该算法的主要思想是通过在目标分布和建议分布之间建立一种接受-拒绝机制，来有效生成所需样本。

查看大图

学习 Metropolis-Hastings 时，先看目标分布、候选分布、接受概率、链更新、 burn-in 和样本混合情况。

目标分布

假设我们要从一个复杂的概率分布 $p(x)$ 中采样。通常，目标分布为后验分布，通常是难以直接采样的。

提议分布

为了生成样本，我们引入一个可行的提议分布 $q(x'|x)$，其用于从当前位置 $x$ 中生成一个新样本 $x'$。典型的提议分布可以是对称的，如正态分布。

接受率

在生成新样本后，我们计算接受率 $α$，决定是否接受新样本。接受率的计算公式为：

$$\alpha = \min\left(1, \frac{p(x')q(x|x')}{p(x)q(x'|x)}\right)$$

• 若 $α \geq 1$，则无条件接受新样本。
• 否则，以概率 $α$ 接受新样本。

算法步骤

1. 选择初始点 $x_0$。
2. 迭代 $n$ 次：
   • 从提议分布 $q(x'|x)$ 中生成新样本 $x'$。
   • 计算接受率 $α$。
   • 生成一个均匀随机数 $u \sim \text{Uniform}(0, 1)$。
   • 如果 $u < α$，则接受 $x'$, 令 $x_{t+1} = x'$；否则，保留旧样本，令 $x_{t+1} = x$。

算法示例

让我们通过一个具体的案例，使用Metropolis-Hastings算法从一个目标分布中采样。

查看大图

进入《Metropolis-Hastings算法》正文前，可以先扫一遍配图：它在问什么、要分清哪些概念、哪一步值得动手、最后用什么标准验收。

目标分布示例

假设我们的目标是从一个一维的正态分布中采样，设其均值为 $\mu=0$，标准差为 $\sigma=1$。我们设定提议分布为当前点附近的正态分布，例如 $q(x'|x) \sim \mathcal{N}(x, 1)$。

Python实现

下面是使用Python实现Metropolis-Hastings算法的一个示例代码：

import numpy as np
import matplotlib.pyplot as plt

# 目标分布: 标准正态分布
def target_distribution(x):
    return np.exp(-0.5 * x**2) / np.sqrt(2 * np.pi)

# 提议分布: 正态
def proposal_distribution(x):
    return np.random.normal(x, 1)

# Metropolis-Hastings算法
def metropolis_hastings(iterations):
    samples = []
    x = 0  # 初始化
    for _ in range(iterations):
        x_new = proposal_distribution(x)
        acceptance_ratio = target_distribution(x_new) / target_distribution(x)
        
        if np.random.rand() < acceptance_ratio:
            x = x_new
        
        samples.append(x)
    return np.array(samples)

# 生成样本
n_samples = 10000
samples = metropolis_hastings(n_samples)

# 绘图
plt.figure(figsize=(12, 6))
plt.hist(samples, bins=30, density=True, alpha=0.6, color='g')

# 绘制目标分布
x = np.linspace(-4, 4, 1000)
plt.plot(x, target_distribution(x), 'r', lw=2)
plt.title('Metropolis-Hastings Sampling')
plt.xlabel('x')
plt.ylabel('Density')
plt.legend(['Target Distribution', 'Samples'])
plt.show()

这段代码首先定义了目标分布和提议分布。然后实施了Metropolis-Hastings算法，并对生成的样本进行了可视化。你可以看到，样本的分布逐渐趋近于目标分布。

查看大图

复习《MetropolisHastings算法》时，建议把关键概念、操作步骤和可见结果放在同一页里回看。

查看大图

练习《MetropolisHastings算法》时，建议把输入条件、处理动作和可见结果写在一起，方便下次复查。

结论

Metropolis-Hastings算法是MCMC的一个强大工具，能够在无法直接采样的情况下，从复杂的后验分布中生成样本。本节中，我们不仅讨论了算法的原理及其步骤，还提供了实际代码的示例，帮助你在实际应用中理解和实现这一算法。

在下一篇中，我们将探讨贝叶斯学习在实际中的应用案例，深入理解这一理论如何通过实际问题得以体现。