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.gal...@gmail.com |
Read me: www.talgalili.com (Hebrew) | www.biostatistics.co.il (Hebrew) |
www.r-statistics.com (English)
----------------------------------------------------------------------------------------------
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
