21 Metropolis-Hastings算法
在马尔可夫链蒙特卡洛(MCMC)方法中,Metropolis-Hastings算法
是一个重要的产生采样方法,广泛应用于贝叶斯学习和统计推断的领域。它是一种构建马尔可夫链的方法,通过从一个指定的概率分布中生成样本,来近似该分布。接下来,我们将详细探讨该算法的原理、步骤和应用案例。
算法原理
Metropolis-Hastings算法
是Gibbs采样
的推广,适用于更为复杂的联合分布。该算法的主要思想是通过在目标分布和建议分布之间建立一种接受-拒绝机制,来有效生成所需样本。
目标分布
假设我们要从一个复杂的概率分布 中采样。通常,目标分布为后验分布
,通常是难以直接采样的。
提议分布
为了生成样本,我们引入一个可行的提议分布 ,其用于从当前位置 中生成一个新样本 。典型的提议分布可以是对称的,如正态分布。
接受率
在生成新样本后,我们计算接受率 ,决定是否接受新样本。接受率的计算公式为:
- 若 ,则无条件接受新样本。
- 否则,以概率 接受新样本。
算法步骤
- 选择初始点 。
- 迭代 次:
- 从提议分布 中生成新样本 。
- 计算接受率 。
- 生成一个均匀随机数 。
- 如果 ,则接受 , 令 ;否则,保留旧样本,令 。
算法示例
让我们通过一个具体的案例,使用Metropolis-Hastings算法
从一个目标分布中采样。
目标分布示例
假设我们的目标是从一个一维的正态分布中采样,设其均值为 ,标准差为 。我们设定提议分布为当前点附近的正态分布,例如 。
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算法
,并对生成的样本进行了可视化。你可以看到,样本的分布逐渐趋近于目标分布。
结论
Metropolis-Hastings算法
是MCMC的一个强大工具,能够在无法直接采样的情况下,从复杂的后验分布中生成样本。本节中,我们不仅讨论了算法的原理及其步骤,还提供了实际代码的示例,帮助你在实际应用中理解和实现这一算法。
在下一篇中,我们将探讨贝叶斯学习
在实际中的应用案例,深入理解这一理论如何通过实际问题得以体现。