从 R 中的核密度估计中获取值

Question

我正在尝试获取 R 中股票价格对数的密度估计。我知道我可以使用

plot(density(x))

绘制它。 然而，我实际上想要该函数的值。

我正在尝试实现核密度估计公式。 这是我到目前为止所拥有的：

a <- read.csv("boi_new.csv", header=FALSE)
S = a[,3] # takes column of increments in stock prices
dS=S[!is.na(S)] # omits first empty field

N = length(dS)                  # Sample size
rseed = 0                       # Random seed
x = rep(c(1:5),N/5)             # Inputted data

set.seed(rseed)   # Sets random seed for reproducibility

QL <- function(dS){
    h = density(dS)$bandwidth
    r = log(dS^2)
    f = 0*x
    for(i in 1:N){
        f[i] = 1/(N*h) * sum(dnorm((x-r[i])/h))
    }
    return(f)
}

QL(dS)

任何帮助将不胜感激。 已经这样好几天了！

Answer 1

您可以直接从

density

函数中提取值：

x = rnorm(100)
d = density(x, from=-5, to = 5, n = 1000)
d$x
d$y

或者，如果您确实想编写自己的核密度函数，这里有一些代码可以帮助您入门：

1. 设置点

   z

   和

   x

   范围：

   z = c(-2, -1, 2)
   x = seq(-5, 5, 0.01)
   

2. 现在我们将把点添加到图表中

   plot(0, 0, xlim=c(-5, 5), ylim=c(-0.02, 0.8), 
        pch=NA, ylab="", xlab="z")
   for(i in 1:length(z)) {
      points(z[i], 0, pch="X", col=2)
   }
    abline(h=0)
   

3. 在每个点周围放置法线密度：

   ## Now we combine the kernels,
   x_total = numeric(length(x))
   for(i in 1:length(x_total)) {
     for(j in 1:length(z)) {
       x_total[i] = x_total[i] + 
         dnorm(x[i], z[j], sd=1)
     }
   }
   

   并将曲线添加到图中：

   lines(x, x_total, col=4, lty=2)
   

4. 最后，计算完整的估算：

   ## Just as a histogram is the sum of the boxes, 
   ## the kernel density estimate is just the sum of the bumps. 
   ## All that's left to do, is ensure that the estimate has the
   ## correct area, i.e. in this case we divide by $n=3$:
   
   plot(x, x_total/3, 
          xlim=c(-5, 5), ylim=c(-0.02, 0.8), 
          ylab="", xlab="z", type="l")
   abline(h=0)
   

   这对应于

   density(z, adjust=1, bw=1)
   

上面的图给出：

Answer 2

这里有一个评估示例 无需从头开始创建自己的密度函数。假设你 有一些数据

X

，并且想要在新的 x 中评估核密度估计器，请调用

xnew

X=c(-1/2,0,1/2)
xnew=5
fdensingle=function(xnew){ density(x=X,
          kernel="epanechnikov",
          bw=.5,
          #using parameter to evaluate in xnew 
          from=xnew,
          to=xnew,
          n=1   
          )$y
}

技巧是对变量

from

和

to

使用相同的限制，加上

n =1

参数。

此外，您可以将此函数转换为其矢量化形式：

fden=function(xnew){
    unlist(lapply(xnew,fdensingle))
}

最后，让我们检查一下结果是否相同（至少从观看情节来看）。

# check the results are the same 
plot(seq(-3,3,by=0.1),fden(seq(-3,3,by=0.1)))
lines(density(x=X,
      kernel="epanechnikov",
      bw=.5)
)