我试图对一系列可能不会正态分布的分布产生基于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 - 为什么stat_density(R; ggplot2)和gaussian_kde(Python; scipy)不同?-LMLPHP

我的python代码是:

kde = stats.gaussian_kde(data.ravel())
kde.set_bandwidth(bw_method='silverman')


python - 为什么stat_density(R; ggplot2)和gaussian_kde(Python; scipy)不同?-LMLPHP

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

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

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

最佳答案

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

scipy docs说Silverman的规则是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版本。

所有这三个来源都引用了相同的参考文献:Silverman,B.W. (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 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)在各种密度下都非常适用,并且评估起来很简单。对于许多目的来说,这当然是窗口宽度的适当选择,对于其他目的,这将是后续微调的良好起点。


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

关于python - 为什么stat_density(R; ggplot2)和gaussian_kde(Python; scipy)不同?,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/55366188/

10-12 18:52