How to compute Integrated Squared Error for kernel density estimation in R

Question

I am working on R for kernel density estimation. I am testing different kernels and I need to evaluate them. I use next code for density estimation:

set.seed(1)
x1 <- rnorm(250) 
#Computing of bandwidth
h1 <- bw.ucv(x1)
#Estimate
est <- density(x1,kernel = "gaussian",bw = h1)
#Values
est$y

I am using the function density with gaussian kernel and bandwidth computed from unbiased crossvalidation. The est$y values are the estimations from kernel. I have to evaluate these estimations with Integrated Squared Error. This is expressed as $ISE=\int\{\hat{f}(x)-f(x)\}^2$. I know $\hat{f}(x)$ is y from est but I have to compute all $ISE$. Is there any way to compute this on R or any package? I have to test this for multiple kernels. Many thanks for your help.

Answer 1

Since you designed $f$ on spec and therefore know what is exactly, you can just integrate $(\hat{f}(x)-f(x))^2$ numerically using quadrature. Your values of the limits and $dx$ should agree with your set values in the density() call, or with the default values. Either way, the values can be extracted from your density() object, which you call est.

I think you will find that the various kernels perform very similarly. For details on KDE, the best reference is Bruce Hansen's unfinished opus. 17 years and counting! Most recent change was last month! (Go, Bruce, go!) See chapter 17, in particular 17.9 on kernel selection. https://www.ssc.wisc.edu/~bhansen/econometrics/Econometrics.pdf

N.B. Because density estimates have bias, and because bias is proportional to the curvature of $f$ ( $f"$, the second derivative) you should calculate the KDE of the derivatives in order to make a bias-corrected estimate. And also note that the bias=0 at inflection points, so if there is something you can do only using inflection points, then that would be a computationally cheap/easy thing to do.

Answer 2

As suggested in a previous answer, numerical integration is a solution in this case. Below you will find 2 ways of achieving this using R. The first solution uses an approximation of your kernel density estimates using spline interpolation. A second alternative is to use a trapezoid approximation (one easy to use function can be find in the R pkg sfsmisc).
The term $f(x)$ is computed for the "x"-values in the est object i.e. using dnorm(est$x)

Hope this helps.

set.seed(1)
x1 = rnorm(250) 
#Computing of bandwidth
h1 = bw.ucv(x1)
#Estimate
est = density(x1,kernel = "gaussian",bw = h1)

# Create approx func obj and integrate()
splxy = splinefun(est$x, (est$y - dnorm(est$x))^2)
integrate(splxy, lower = min(x1), upper = max(x1)) # res 0.0006681713

# Simple numerical integration in pkg sfsminsc
sfsmisc::integrate.xy(x = est$x, (est$y - dnorm(est$x))^2) # res 0.0006682952