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

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

统计文档显示 here nrd0 是 bw 调整的 silverman 方法。

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

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

Answer 1

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

scipy docs say that Silverman's rule is 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 给出：

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 版本。

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

阅读 Silverman 论文，第 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 描述了这种方法：

While (3.28) will work well if the population really is normally distributed, it may oversmooth somewhat if the population is multimodal... as the mixture becomes more strongly bimodal the formula (3.28) will oversmooth more and more, relative to the optimal choice of smoothing parameter.

scipy 文档 reference this weakness，说明：

It includes automatic bandwidth determination. The estimation works best for a unimodal distribution; bimodal or multi-modal distributions tend to be oversmoothed.

但是，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 将此方法描述为：

The best of both possible worlds... In summary, the choice ([eqn] 3.31) for the smoothing parameter will do very well for a wide range of densities and is trivial to evaluate. For many purposes it will certainly be an adequate choice of window width, and for others it will be a good starting point for subsequent fine-tuning.

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

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

Why do stat_density (R; ggplot2) and gaussian_kde (Python; scipy) differ?

python

r

scipy

ggplot2

kernel-density