为什么 stat_density (R; ggplot2) 和 gaussian_kde 不同？

Question

我试图在一系列可能不是正态分布的分布上生成基于 KDE 的 PDF 估计。我喜欢 R 中 ggplot 的 stat_density 似乎可以识别频率中每个增量颠簸的方式，但无法通过 Python 的 scipy-stats-gaussian_kde 方法复制这一点，这似乎过于平滑。我已经按如下方式设置了我的 R 代码：ggplot(test, aes(x=Val, color = as.factor(Class), group=as.factor(Class))) +             stat_density(geom='line',kernel='gaussian',bw='nrd0'                                                             #nrd0='Silverman'                                                            ,size=1,position='identity')我的python代码是：kde = stats.gaussian_kde(data.ravel())kde.set_bandwidth(bw_method='silverman')统计文档在这里显示 nrd0 是用于 bw 调整的silverman 方法。基于上面的代码，我使用相同的内核（高斯）和带宽方法（西尔弗曼）。谁能解释为什么结果如此不同？

Answer 1

关于什么是西尔弗曼规则似乎存在分歧。TL; DR - scipy 使用更糟糕的规则版本，该规则仅适用于正态分布的单峰数据。R 使用更好的版本，它是“两全其美”并且“适用于广泛的密度”。

该SciPy的文档说，西尔弗曼的规则是为实现：

def silverman_factor(self):

    return power(self.neff*(self.d+2.0)/4.0, -1./(self.d+4))

哪里d是维度数（在您的情况下neff为1），是有效样本大小（点数，假设没有权重）。所以 scipy 带宽是(n * 3 / 4) ^ (-1 / 5)（乘以标准偏差，以不同的方法计算）。

相比之下，R的stats封装文档描述西尔弗曼的方法为“最小0.9倍的标准偏差的和四分位范围由1.34倍样品量划分到负五分之一功率”，这也可作为R代码验证，打字bw.nrd0在控制台给出：

function (x) 

{

    if (length(x) < 2L) 

        stop("need at least 2 data points")

    hi <- sd(x)

    if (!(lo <- min(hi, IQR(x)/1.34))) 

        (lo <- hi) || (lo <- abs(x[1L])) || (lo <- 1)

    0.9 * lo * length(x)^(-0.2)

}

另一方面，维基百科将“西尔弗曼的经验法则”作为估计器的众多可能名称之一：

1.06 * sigma * n ^ (-1 / 5)

维基百科版本相当于 scipy 版本。

所有三个来源（scipy 文档、维基百科和 R 文档）都引用了相同的原始参考资料：Silverman, BW (1986)。统计和数据分析的密度估计。伦敦：查普曼和霍尔/CRC。p. 48. ISBN 978-0-412-24620-3。维基百科和 R 专门引用了第 48 页，而 scipy 的文档没有提到页码。（我已向 Wikipedia 提交了编辑以将其页面引用更新为第 45 页，见下文。）

读西尔弗曼纸，45页，方程3.28上是什么是维基百科的文章中使用：(4 / 3) ^ (1 / 5) * sigma * n ^ (-1 / 5) ~= 1.06 * sigma * n ^ (-1 / 5)。Scipy 使用相同的方法，重写(4 / 3) ^ (1 / 5)为等效的(3 / 4) ^ (-1 / 5). Silverman 描述了这种方法：

虽然（3.28）在群体确实是正态分布的情况下会很好地工作，但如果群体是多峰的，它可能会有些过度平滑……随着混合物变得越来越强双峰，相对于最佳选择，公式（3.28）将越来越过平滑平滑参数。

scipy 文档引用了这个弱点，指出：

它包括自动带宽确定。该估计最适用于单峰分布；双峰或多峰分布往往过于平滑。

然而，Silverman 的文章继续，改进了 scipy 使用的方法，以获得 R 和 Stata 使用的方法。在第 48 页，我们得到方程 3.31：

h = 0.9 * A * n ^ (-1 / 5)

# A defined on previous page, eqn 3.30

A = min(standard deviation, interquartile range / 1.34)

Silverman 将这种方法描述为：

两全其美……总而言之，平滑参数的选择 ([eqn] 3.31) 将适用于广泛的密度范围，并且易于评估。对于许多目的，它肯定是窗口宽度的充分选择，对于其他人来说，它将是后续微调的良好起点。

因此，Wikipedia 和 Scipy 似乎使用了 Silverman 提出的具有已知弱点的估计器的简单版本。R 和 Stata 使用更好的版本。