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

Question

我试图对一系列可能不会正态分布的分布产生基于KDE的PDF估计。

我喜欢ggplot在R中的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是西尔弗曼方法为体重调整。

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

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

Answer 1

关于西尔弗曼法则，似乎存在分歧。

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

另一方面，Wikipedia将“ Silverman的经验法则”作为估计器的许多可能名称之一：

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

维基百科版本等同于scipy版本。

所有这三个来源都引用了相同的参考文献：Silverman，BW（1986）。统计和数据分析的密度估计。伦敦：查普曼和霍尔/ CRC。p。48. ISBN 978-0-412-24620-3。Wikipedia和R特别引用了第48页，而scipy的文档未提及页码。（我已向Wikipedia提交了修改，以将其页面引用更新为第45页，请参见下文。）

附录

我找到了Silverman参考资料的PDF。

在第45页上，等式3.28是Wikipedia文章中使用的：(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）可能会在某种程度上过度平滑...随着混合物变得更强的双峰，相对于最佳选择，公式（3.28）将会越来越平滑平滑参数。

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

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

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）在各种密度下都非常适用，并且评估起来很简单。对于许多目的来说，这当然是窗口宽度的适当选择，对于其他目的，这将是后续微调的良好起点。

因此，似乎Wikipedia和Scipy使用Silverman提出的估算器的简单版本。R和Stata使用更精致的版本。