为什么stat_density（R; ggplot2）和gaussian_kde（Python; scipy）有所不同？

Question

我试图在一系列可能无法正常分发的发行版上生成基于KDE的PDF估计。

我喜欢Rg中的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')

统计文档显示here，nrd0是bw adjust的silverman方法。

基于上面的代码，我使用相同的kernals（高斯）和带宽方法（Silverman）。

谁能解释为什么结果如此不同？

Answer 1

关于西尔弗曼的规则似乎存在分歧。

scipy docs说西尔弗曼的规则是implemented as：

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

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

相比之下，R的stats package docs将Silverman的方法描述为“标准偏差的最小值的0.9倍，四分位数范围除以样本大小的1.34倍到负的五分之一”，这也可以在R代码中验证，输入bw.nrd0 in控制台给出：

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)
}

另一方面，Wikipedia将“西尔弗曼的经验法则”作为估算器的许多可能名称之一：

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

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

所有三个来源引用相同的参考：Silverman，B.W。 （1986年）。统计和数据分析的密度估计。伦敦：Chapman＆Hall / CRC。页。 48. ISBN 978-0-412-24620-3。维基百科和R特别引用了第48页，而scipy的文档没有提到页码。 （我已经向维基百科提交了一个编辑，以更新其页面引用到第45页，见下文。）

---

附录

我找到了Silverman参考文献的PDF。

在第45页，公式3.28是维基百科文章中使用的：(4 / 3) ^ (1 / 5) * sigma * n ^ (-1 / 5) ~= 1.06 * sigma * n ^ (-1 / 5)。 Scipy使用相同的方法，重写(3 / 4) ^ (-1 / 5)作为等效的(4 / 3) ^ (1 / 5)。 Silverman将此方法描述为：

如果人口真正是正态分布的话，（3.28）会很好地工作，如果人口是多模态的话，它可能会有些过度...因为混合物变得更加强烈双峰，公式（3.28）会越来越过分，相对于最优选择平滑参数。

scipy docs reference this weakness，陈述：

它包括自动带宽确定。估算最适合单峰分布;双峰或多峰分布往往过度平滑。

Silverman文章继续激励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）对于各种密度都会很好，并且很容易评估。出于许多目的，它肯定是窗口宽度的适当选择，对于其他目的，它将是后续微调的良好起点。

因此，似乎维基百科和Scipy使用了Silverman提出的简单估计器。 R和Stata使用更精致的版本。