I have entered into the R code of the density function, and I have noticed the following strange lines:
lo <- from - 4 * bw
up <- to + 4 * bw
To my understanding, they mean that the density is estimated on the interval [from - 4*bw, to + 4*bw] instead of [from, to].
To really understand what is going on, I have created a densityfun function, which copy-pastes the code of density, except at the end in order to return a function (you can find the R code on GitHub here). Then I get the following:
set.seed(1)
x <- rbeta(10000, 0.5, 0.5)
f <- densityfun(x, from = 0, to = 1)
f(-0.01) # 1.135904
Surprise: f(-0.01) is nonzero!
It also implies that the integral of f on [0, 1] is not 1:
integrate(f, 0, 1) # 0.8787954
integrate(f, -0.1, 1.1) # 0.997002
So why is the density function written that way (is it a bug?), and what could I do to avoid this behaviour (to have f(-0.01) = 0 in this example) without loosing any mass on f (to have integrate(f, 0, 1) approximately equal to 1 in this example)?
Thanks!
EDIT: I have changed a bit the values used in the example.