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

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.

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.

• +1 for the reference to Bruce Hansen's book alone!
– jhin
Aug 20 at 16:49

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

• Many thanks dear @Gi_F. it is what I was looking for. I am going to accept your answer. Just a small question if x1 were from a exponential distribution, I should use dexp to compute splxy object, right?
– Duck
Feb 10 '19 at 11:00
• I think that would work. On the other hand give a close look to the smoothed density obtained with the kde before computing the ISE. Boundaries could be an issue . Feb 10 '19 at 16:48