基于Metropolis Sampling拟合模型估算mu与sigma的技术疑问

首先，先还原你的问题背景和伪代码：

我正在阅读《Bayesian Analysis in Python》一书，该书侧重PyMC3包的使用，但对背后的理论讲解模糊，我对此十分困惑。现有数据：

data = np.array([51.06, 55.12, 53.73, 50.24, 52.05, 56.40, 48.45, 52.34, 55.65, 51.49, 51.86, 63.43, 53.00, 56.09, 51.93, 52.31, 52.33, 57.48, 57.44, 55.14, 53.93, 54.62, 56.09, 68.58, 51.36, 55.47, 50.73, 51.94, 54.95, 50.39, 52.91, 51.5, 52.68, 47.72, 49.73, 51.82, 54.99, 52.84, 53.19, 54.52, 51.46, 53.73, 51.61, 49.81, 52.42, 54.3, 53.84, 53.16])

我希望使用Metropolis Sampling拟合模型以估算mu和sigma，并写出了伪代码：

M, S = 50, 1
G = 1 # These are priors right?
mu = stats.norm(loc=M, scale=S)
sigma = stats.halfnorm(scale=G)
target = stats.norm
steps = 1000
mu_samples = [50]
sigma_samples = [1]
for i in range(steps):
    # proposed sample...
    mu_i, sigma_i = mu.rvs(), sigma.rvs()
    # Something happens here
    # How do I calculate the likelidhood??
    "..."
    # some evaluation of a likelihood ratio??
    a = "some"/"ratio"
    acceptance_bar = np.random.random()
    if a > acceptance_bar:
        mu_samples.append(mu_i)
        sigma_samples.append(sigma_i)

请问这段伪代码中缺失了哪些关键步骤？

好的，咱们来一步步拆解你的伪代码里漏掉的核心逻辑——本质上是没把Metropolis算法的核心规则落地，具体缺失的关键步骤如下：

1. 错误的候选样本生成逻辑 #

你现在直接用先验分布mu.rvs()和sigma.rvs()生成候选值，这完全不对！Metropolis算法的候选样本应该是基于当前样本的局部提议分布（比如正态分布），而不是直接从先验里抽取。正确的生成方式应该是：

# 先获取当前链的最后一个样本
current_mu = mu_samples[-1]
current_sigma = sigma_samples[-1]
# 用提议分布生成候选（步长可以根据接受率调整，这里给个示例值）
mu_candidate = np.random.normal(loc=current_mu, scale=1.0)
# sigma是正数，所以要确保候选值大于0，避免无效参数
sigma_candidate = np.random.normal(loc=current_sigma, scale=0.5)
sigma_candidate = max(sigma_candidate, 1e-6)  # 防止sigma为0或负数

2. 缺失似然函数的计算（后验概率的核心组成） #

你需要计算候选参数下的似然值，也就是给定mu_candidate和sigma_candidate时，观测数据出现的概率。因为数据是独立同分布的正态样本，用对数似然更稳定（避免乘积过小导致数值下溢）：

# 候选参数的对数似然：所有数据点的正态对数概率密度之和
log_likelihood_candidate = np.sum(stats.norm.logpdf(data, loc=mu_candidate, scale=sigma_candidate))
# 当前参数的对数似然，用来计算后验比值
log_likelihood_current = np.sum(stats.norm.logpdf(data, loc=current_mu, scale=current_sigma))

3. 缺失先验概率的计算（贝叶斯后验=似然×先验） #

后验概率是似然乘以先验，所以还要计算候选参数和当前参数的先验概率（同样用对数形式）：

# 候选参数的对数先验：mu的先验 + sigma的先验
log_prior_candidate = stats.norm.logpdf(mu_candidate, loc=M, scale=S) + stats.halfnorm.logpdf(sigma_candidate, scale=G)
# 当前参数的对数先验
log_prior_current = stats.norm.logpdf(current_mu, loc=M, scale=S) + stats.halfnorm.logpdf(current_sigma, scale=G)

4. 正确计算接受概率（Metropolis准则） #

接受概率a是后验概率的比值（候选参数的后验 / 当前参数的后验），用对数计算后再取指数更安全：

# 对数后验的差值：(候选似然+候选先验) - (当前似然+当前先验)
log_posterior_ratio = (log_likelihood_candidate + log_prior_candidate) - (log_likelihood_current + log_prior_current)
# 接受概率不能超过1，所以取指数后和1取最小值
acceptance_prob = min(1, np.exp(log_posterior_ratio))

5. 处理拒绝候选的情况 #

你现在只有接受候选时才添加样本，但拒绝时应该保留当前样本，否则采样链会断裂，无法反映后验分布的连续性：

acceptance_bar = np.random.random()
if acceptance_bar < acceptance_prob:
    # 接受候选，添加新样本
    mu_samples.append(mu_candidate)
    sigma_samples.append(sigma_candidate)
else:
    # 拒绝候选，重复当前样本
    mu_samples.append(current_mu)
    sigma_samples.append(current_sigma)

补充：可选但重要的优化步骤 #

• 燃烧期（burn-in）：采样结束后，建议丢弃前10%-20%的样本（比如前200步），因为初始样本可能还没收敛到真实的后验分布。
• 步长调整：如果采样的接受率太高（>50%）或太低（<20%），要调整提议分布的步长，让接受率维持在20%-50%之间，保证采样效率。

内容的提问来源于stack exchange，提问作者RSHAP