Расхождения в ядерной оценке плотности () по сравнению с расчетами с нуля

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

Базовый расчет был построен в соответствии с деталями, изложенными в разделе 3.3.6 в «Статистических методах в науках об атмосфере», 3-е издание, Дэниела Уилкса. Уравнение 3.13 из учебника Уилкса с ядром Гаусса, установленным как введите здесь описание изображения и t введите здесь описание изображения

Впрочем, и вот моя проблема.

Затем я рисую их вместе...

plot(y=fx0,x=x0, type="l", ylim=c(0,0.07))
lines(x=auto.dens$x, y=auto.dens$y, col="red")

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

!Эти два расчета явно различаются!

Я не понял, как работает функция плотности? Почему я не могу вычислить те же результаты с нуля? Почему мой оценщик ядра дает разные результаты? Почему мои результаты менее гладкие?

Мне нужно построить и применить более гладкое ядро ​​(не только плотность) к гораздо более сложному набору данных, и я сделал этот небольшой пример только для того, чтобы убедиться, что я делаю то же самое, что и автоматизированные функции, и действительно не ожидал, что Эта проблема. Я пробовал все виды вещей, и просто не понимаю, почему я получаю другой результат.

Всем заранее спасибо за прочтение и любые комментарии, маленькие или большие.

Редактировать: 13:40 29/11/2016 Решение, подробно описанное в ответе ниже введите здесь описание изображения


r gaussian probability-density kernel-density
person Kate2808    schedule 29.11.2016    source источник

Ответы (1)


arrow_upward
2
arrow_downward

Вам не нужно 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")
person Kota Mori    schedule 29.11.2016
comment
Спасибо! Это устраняет проблему, и, очевидно, это просто мое понимание математики из учебника, которое не соответствует коду. Я так рад, что это такая простая проблема! - person Kate2808; 29.11.2016