R help - Oct 2013 - Reproducing density(..., kernel='rectangular', width=1) using convolve? (for educational purposes)

Hello dear R-help members,

I am currently teaching students on kernel density estimators, and I was
hoping to show them an animation of the convolution which is taking place
in the "density" function. I wrote a (somewhat clunky) piece of code
to
reproduce the density function, but I seem to get rather different results
(which means I am missing something).

If any of you are in the mood to help and see what is wrong (or if you have
your own piece of code which already does this), it would be lovely for
future students who are studying this method.

Here is the code I got thus far:


#### The function to reproduce
set.seed(13413)
some_data <- c(round(rnorm(100)))
plot(density(some_data, kernel='rectangular', width=1, n = 10**5), lwd =
3)

#### The attempt
convolve_new <- function(data_range, kernel_range, data_fun, kernel_fun,
normalize = TRUE) {
   # data_range, kernel_range - MUST have the same length

   # if we have the range of a kernel function and of the data function
   # we would like to have one move over the other without any problems.
   y_data_fun <- data_fun(data_range)
   y_kernel_fun <- kernel_fun(kernel_range)
   n_kernel <- length(y_kernel_fun)
   n_data <- length(y_kernel_fun)
   conv_vec <- convolve(y_data_fun, rev(y_kernel_fun), type = "o")

   items_to_remove <- floor(n_kernel/2)
   conv_vec <- head(conv_vec,-items_to_remove)
   conv_vec <- tail(conv_vec,-items_to_remove)

   dx <- ifelse(normalize, mean(diff(data_range)), 1)# this is dx
   conv_vec  * dx
}


w_R_01 <- function(u) {
   ifelse(abs(u)<=.5 & u >=0,1,0)
}


# fun1 <- dnorm # kernel
fun1 <- w_R # kernel
# fun2 <- w_R_01 # data
fun2 <-  # data
   function(x) {
      sapply(x, function(y) sum(some_data %in% y))
   }
# fun2(1)


# the_data_range <- seq(-2,2, by = .01)
the_data_range <- sort(c(seq(min(some_data),max(some_data), by = .01),
some_data))
the_kernel_range <- the_data_range
convo_the_range <- convolve_new(the_data_range, the_kernel_range, fun2,
fun1)
the_animation_range <- seq(-2,2, by = .2)
# plot(convo_the_range/sum(convo_the_range)~the_range, type = "l")

for(i in the_animation_range) {
#    curve(fun2,
#          from = min(the_data_range), to = max(the_data_range), ylim c(0,3),
#          n=10**3)
   plot(fun2(the_data_range) ~ the_data_range, type= "h")
   abline(v = i, lty = 2, col = "purple", lwd = 2)
   temp_fun1 <- function(x) {fun1(x-i)}
   curve(temp_fun1,
         from = min(the_animation_range), to = max(the_animation_range),
add = TRUE,
         n=10**3, lwd = 3, col = 2)
#    temp_convolution_range <- function(x) {convolution_range(x, fun1
temp_fun1, fun2 = fun2)}
#    #    curve(convolution_range, from = -3, to =3, add = TRUE)
   ss<- the_data_range<=i
   lines(convo_the_range[ss]~the_data_range[ss], type = "l", col = 3,
lwd 3)
   rug(some_data)
}



Thanks,
Tal


----------------Contact
Details:-------------------------------------------------------
Contact me: Tal.Galili@gmail.com |
Read me: www.talgalili.com (Hebrew) | www.biostatistics.co.il (Hebrew) |
www.r-statistics.com (English)
----------------------------------------------------------------------------------------------

	[[alternative HTML version deleted]]