Получение значений из оценки плотности ядра в R

Я пытаюсь получить оценки плотности для журнала цен на акции в 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)

Любая помощь будет высоко оценена. Это было в течение нескольких дней!

Ответы

Ответ 1

Вы можете вывести значения непосредственно из функции density:

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

В качестве альтернативы, если вы действительно хотите написать свою собственную функцию плотности ядра, здесь приведен код для запуска:

  • Задайте точки z и x:

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

  • Теперь мы добавим точки в график

    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)

  • Поместите нормальную плотность вокруг каждой точки:

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

  • Наконец, вычислим полную оценку:

    ## Just as a histogram is the sum of the boxes, 
## the kernel density estimate is just the sum of the bumps. 
## All that 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)

Приведенные выше графики дают:

enter image description here