density() kernel estimator 与 scratch 计算的差异
Discrepancies in the density() kernel estimator compared to calculations by scratch
我正在尝试计算高斯核密度,为了测试我对 density()
函数的了解,我决定从头开始计算它并比较两个结果。
然而,他们没有提供相同的答案。
我从现有数据集开始
xi <- mtcars$mpg
并且可以绘制出这个数据的核密度,如下
plot(density(xi, kernel = "gaussian"))
它提供了这个...
然后我从这个计算中抓取一些细节,使我的计算是一致的。
auto.dens <- density(xi, kernel = "gaussian")
h <- auto.dens$bw # bandwidth for kernel
x0 <- auto.dens$x # points for prediction
然后我自己计算高斯核密度,我有
在循环中完成此操作,以便阅读更清晰。
fx0 <- NULL
for (j in 1:length(x0)){
t <- abs(x0[j]-xi)/h
K <- (1/sqrt(2*pi))*exp(-(t^2)/2)
fx0 <- c(fx0,sum(K*t)/(length(t)*h))
}
基本计算是按照 Daniel Wilks 的《大气科学统计方法》第 3 版第 3.3.6 节中的详细信息构建的。
高斯核设置为 and t being
但是,这是我的问题。
然后我把两者放在一起...
plot(y=fx0,x=x0, type="l", ylim=c(0,0.07))
lines(x=auto.dens$x, y=auto.dens$y, col="red")
密度函数的输出(红色)和我的计算(黑色),我得到
!这两个计算明显不同!
我是否没有理解密度函数的工作原理?为什么我不能设法从头开始计算相同的结果?为什么我的核估计器提供不同的结果?为什么我的结果不太流畅?
我需要构建内核平滑器(不仅仅是密度)并将其应用于更复杂的数据集,并且只做了这个小示例以确保我在做与自动化功能相同的事情,而且确实如此'期待有这个问题。我已经尝试了各种各样的事情,只是不明白为什么我会得到不同的结果。
提前感谢大家的阅读和任何意见,无论大小。
编辑:13:402016 年 11 月 29 日
解决方案详见下方答案
您不需要sum(K*t)
,只需sum(K)
。
xi <- mtcars$mpg
plot(density(xi, kernel = "gaussian"), lwd = 2)
auto.dens <- density(xi, kernel = "gaussian")
h <- auto.dens$bw # bandwidth for kernel
x0 <- auto.dens$x # points for prediction
fx0 <- NULL
for (j in 1:length(x0)) {
t <- abs(x0[j]-xi)/h
K <- (1/sqrt(2*pi))*exp(-(t^2)/2)
fx0 <- c(fx0, sum(K)/(length(t)*h))
}
lines(x0, fx0, col = "red", lty = "dotted")
我正在尝试计算高斯核密度,为了测试我对 density()
函数的了解,我决定从头开始计算它并比较两个结果。
然而,他们没有提供相同的答案。
我从现有数据集开始
xi <- mtcars$mpg
并且可以绘制出这个数据的核密度,如下
plot(density(xi, kernel = "gaussian"))
它提供了这个...
然后我从这个计算中抓取一些细节,使我的计算是一致的。
auto.dens <- density(xi, kernel = "gaussian")
h <- auto.dens$bw # bandwidth for kernel
x0 <- auto.dens$x # points for prediction
然后我自己计算高斯核密度,我有 在循环中完成此操作,以便阅读更清晰。
fx0 <- NULL
for (j in 1:length(x0)){
t <- abs(x0[j]-xi)/h
K <- (1/sqrt(2*pi))*exp(-(t^2)/2)
fx0 <- c(fx0,sum(K*t)/(length(t)*h))
}
基本计算是按照 Daniel Wilks 的《大气科学统计方法》第 3 版第 3.3.6 节中的详细信息构建的。
但是,这是我的问题。
然后我把两者放在一起...
plot(y=fx0,x=x0, type="l", ylim=c(0,0.07))
lines(x=auto.dens$x, y=auto.dens$y, col="red")
密度函数的输出(红色)和我的计算(黑色),我得到
!这两个计算明显不同!
我是否没有理解密度函数的工作原理?为什么我不能设法从头开始计算相同的结果?为什么我的核估计器提供不同的结果?为什么我的结果不太流畅?
我需要构建内核平滑器(不仅仅是密度)并将其应用于更复杂的数据集,并且只做了这个小示例以确保我在做与自动化功能相同的事情,而且确实如此'期待有这个问题。我已经尝试了各种各样的事情,只是不明白为什么我会得到不同的结果。
提前感谢大家的阅读和任何意见,无论大小。
编辑:13:402016 年 11 月 29 日
解决方案详见下方答案
您不需要sum(K*t)
,只需sum(K)
。
xi <- mtcars$mpg
plot(density(xi, kernel = "gaussian"), lwd = 2)
auto.dens <- density(xi, kernel = "gaussian")
h <- auto.dens$bw # bandwidth for kernel
x0 <- auto.dens$x # points for prediction
fx0 <- NULL
for (j in 1:length(x0)) {
t <- abs(x0[j]-xi)/h
K <- (1/sqrt(2*pi))*exp(-(t^2)/2)
fx0 <- c(fx0, sum(K)/(length(t)*h))
}
lines(x0, fx0, col = "red", lty = "dotted")