John Kruschke 在他的《做贝叶斯数据分析》一书中指出，使用 R 中的 JAGS
...来自 MCMC 样本的模式估计可能相当不稳定，因为该估计基于一种平滑算法，该算法可能对 MCMC 样本中的随机颠簸和波纹敏感。（进行贝叶斯数据分析，第 205 页，第 8.2.5.1 节）
虽然我掌握了 Metropolis 算法和 Gibbs 采样等精确形式，但我不熟悉也提到的平滑算法，以及为什么这意味着从 MCMC 样本中对模式的估计是不稳定的。有没有人能够直观地了解平滑算法正在做什么以及为什么它会使模式的估计不稳定？
我手头没有这本书，所以我不确定 Kruschke 使用什么平滑方法，但出于直觉，请考虑这个由标准法线中的 100 个样本组成的图，以及使用 0.1 到 1.0 的各种带宽的高斯核密度估计。（简而言之，高斯 KDE 是一种平滑直方图：它们通过为每个数据点添加高斯来估计密度，平均值为观察值。）
您可以看到，即使平滑创建了单峰分布，其众数通常也低于已知值 0。
此外，这里是使用相同样本的用于估计密度的内核带宽的估计模式（y 轴）图。希望这可以直观地了解估计值如何随平滑参数的变化而变化。
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Feb 1 09:35:51 2017
@author: seaneaster
"""
import numpy as np
from matplotlib import pylab as plt
from sklearn.neighbors import KernelDensity
REAL_MODE = 0
np.random.seed(123)
def estimate_mode(X, bandwidth = 0.75):
kde = KernelDensity(kernel = 'gaussian', bandwidth = bandwidth).fit(X)
u = np.linspace(-3,3,num=1000)[:, np.newaxis]
log_density = kde.score_samples(u)
return u[np.argmax(log_density)]
X = np.random.normal(REAL_MODE, size = 100)[:, np.newaxis] # keeping to standard normal
bandwidths = np.linspace(0.1, 1., num = 8)
plt.figure(0)
plt.hist(X, bins = 100, normed = True, alpha = 0.25)
for bandwidth in bandwidths:
kde = KernelDensity(kernel = 'gaussian', bandwidth = bandwidth).fit(X)
u = np.linspace(-3,3,num=1000)[:, np.newaxis]
log_density = kde.score_samples(u)
plt.plot(u, np.exp(log_density))
bandwidths = np.linspace(0.1, 3., num = 100)
modes = [estimate_mode(X, bandwidth) for bandwidth in bandwidths]
plt.figure(1)
plt.plot(bandwidths, np.array(modes))
Sean Easter 提供了一个很好的答案；这是 Kruschke 书中附带的 R 脚本实际上是如何完成的。该plotPost()函数在名为 的 R 脚本中定义DBDA2E-utilities.R。它显示估计的模式。在函数定义中，有以下两行：
mcmcDensity = density(paramSampleVec)
mo = mcmcDensity$x[which.max(mcmcDensity$y)]
该density()函数带有 R 的基本 stats 包，并实现了 Sean Easter 所描述的那种内核密度过滤器。它具有平滑内核的带宽和要使用的内核类型的可选参数。它默认为高斯内核，并且有一些内部魔法可以找到合适的带宽。该density()函数返回一个对象，该对象具有一个名为的组件，该组件y具有各种值的平滑密度x。上面的第二行代码只是找到最大值的x值y。