# Area under the "pdf" in kernel density estimation in R

I am trying to use the 'density' function in R to do kernel density estimates. I am having some difficulty interpreting the results and comparing various datasets as it seems the area under the curve is not necessarily 1. For any probability density function (pdf) $\phi(x)$, we need to have the area $\int_{-\infty}^\infty \phi(x) dx = 1$. I am assuming that the kernel density estimate reports the pdf. I am using integrate.xy from sfsmisc to estimate the area under the curve.

> # generate some data
> xx<-rnorm(10000)
> # get density
> xy <- density(xx)
> # plot it
> plot(xy)


> # load the library
> library(sfsmisc)
> integrate.xy(xy$x,xy$y)
[1] 1.000978
> # fair enough, area close to 1
> # use another bw
> xy <- density(xx,bw=.001)
> plot(xy)


> integrate.xy(xy$x,xy$y)
[1] 6.518703
> xy <- density(xx,bw=1)
> integrate.xy(xy$x,xy$y)
[1] 1.000977
> plot(xy)


> xy <- density(xx,bw=1e-6)
> integrate.xy(xy$x,xy$y)
[1] 6507.451
> plot(xy)


Shouldn't the area under the curve always be 1? It seems small bandwidths are a problem, but sometimes you want to show the details etc. in the tails and small bandwidths are needed.

It seems that the answer below about the overestimation in convex regions is correct as increasing the number of integration points seems to lessen the problem (I didn't try to use more than $2^{20}$ points.)

> xy <- density(xx,n=2^15,bw=.001)
> plot(xy)


> integrate.xy(xy$x,xy$y)
[1] 1.000015
> xy <- density(xx,n=2^20,bw=1e-6)
> integrate.xy(xy$x,xy$y)
[1] 2.812398


• This looks like a floating point limitation in density(): in using a bandwidth of 1e-6, you are creating (in theory) a collection of 10,000 spikes, each of total mass 1/10000. Those spikes end up being represented primarily by their peaks, without the gaps being adequately characterized. You're merely pushing density() beyond its limits.
– whuber
Commented Aug 10, 2011 at 1:58
• @whuber, by floating point limitation, do you mean limits of the precision, as in using floats would lead to greater overestimation of the error compared to using doubles. I don't think I see how that would happen but would like to see some evidence. Commented Aug 10, 2011 at 14:02
• Your update demonstrates that convexity is not the issue; the issue lies in using too small a value of $n$ in the density calculation.
– whuber
Commented Aug 10, 2011 at 15:04
• Shouldn't the integral value of a proper density estimate be $1$? Commented Sep 17, 2013 at 13:18
• @Anony-Mousse, yes, that is what this question is asking. Why is it not evaluating to 1? Commented Sep 18, 2013 at 20:26

Think about the trapezoid rule integrate.xy() uses. For the normal distribution, it will underestimate the area under the curve in the interval (-1,1) where the density is concave (and hence the linear interpolation is below the true density), and overestimate it elsewhere (as the linear interpolation goes on top of the true density). Since the latter region is larger (in Lesbegue measure, if you like), the trapezoid rule tends to overestimate the integral. Now, as you move to smaller bandwidths, pretty much all of your estimate is piecewise convex, with a lot of narrow spikes corresponding to the data points, and valleys between them. That's where the trapezoid rule breaks down especially badly.
• This explanation does not hold water. The problem is that the density is inadequately discretized, not that the trapezoid rule breaks down badly. integrate() is helpless to get a correct answer because density() does not produce a correct representation. To see this, just inspect xy$x: it has only 512 values intended to represent 10,000 narrow spikes! – whuber Commented Aug 10, 2011 at 2:02 • @whuber, that's what the answer said. The point is you need to use the trapezoidal rule for finite number of samples, and it overestimates the area compared to the true density on a continuous axis according to the kernels. My update at the end of the question expands on it. Commented Aug 10, 2011 at 2:07 • @high No; the trapezoidal rule is working fine. The problem is that it is working with an incorrect discretization of the integrand. You can't possibly have "a lot of narrow spikes corresponding to the data points" when there are 10,000 data points and only 512 values in the density array! – whuber Commented Aug 10, 2011 at 2:10 • Looking at these graphs, I am now thinking that the problem is with density rather than with integrate.xy. With N=10000 and bw=1e-6, you would have to see a comb with a height of each tooth of about 1e6, and the teeth being denser around 0. Instead, you still see a recognizable bell-shaped curve. So density is cheating on you, or at least it should be used differently with tiny bandwidths: n should be about (range of data)/(bw) rather than the default n=512. The intergrator must be picking up one of these huge values that density returns by an unhappy coincidence. Commented Aug 10, 2011 at 16:10 Ive been integrating using the automatic bandwidth and getting answers ranging from .1 to 7 , not particularly trustworty. #### Use the KDE function in the utilties package instead Note: This answer is adapted from a silmiar answer here. The problem you are encountering is that the density function in the stats package gives you an output that is just a vector of density values at a finite number of points --- it is not actually giving you a full density function that you can integrate effectively. This is a reasonable way to give a density output if you just want to plot the density reasonable well, but as you can see, it leads to a value that is slightly different to one when you try to integrate it. To get around this deficiency, you can use the KDE function in the utilties package. This function generates the KDE in the same way as the density function in R,$$^\dagger$$ but instead of producing an output computed over a relatively small set of points, it produces an output that includes the probability functions for the KDE. You can also call the function in such a way that it loads those probability functions directly to the global environment, so that you can easily call them just like any other density function in R. This will give you a full density function that you can call at any set of points, and also a cumulative distribution function where the integration of the density has already been done for you. Below I give an example of how to generate the KDE using this function, and how to call the cumulative distribution function over an arbitrary set of values. As you can see, the KDE function produces a set of probability functions (dkde, pkde, qkde, and rkde) that can be called just like the probability functions for any of the pre-programmed families of distributions. This allows you to compute the cumulative distribution from pkde at any point you want, including points that are far outside the data range used to generate the KDE. #Load the package library(utilities) #Generate some mock data set.seed(1) DATA <- rnorm(40) #Create a KDE and show its output MY_KDE <- KDE(DATA, to.environment = TRUE) MY_KDE Kernel Density Estimator (KDE) Computed from 40 data points in the input 'DATA' Estimated bandwidth = 0.367412 Input degrees-of-freedom = Inf Probability functions for the KDE are the following: Density function: dkde * Distribution function: pkde * Quantile function: qkde * Random generation function: rkde * * This function is presently loaded in the global environment #Call the CDF over a set of points (including points far in the tails) POINTS <- -10:10 pkde(POINTS) [1] 1.489573e-101 4.685332e-78 9.132757e-58 1.112228e-40 8.584043e-27 4.322183e-16 [7] 1.529301e-08 4.819333e-04 3.326124e-02 1.236576e-01 4.251039e-01 8.352227e-01 [13] 9.927698e-01 9.999976e-01 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 [19] 1.000000e+00 1.000000e+00 1.000000e+00  From the outputs of the CDF you can see that the density function integrates to one over the full range. (If you're unsure, just call pkde(Inf) to see the CDF value integrating over the whole range.) $$^\dagger$$ The KDE function has the advantage of giving a more useful output (in my opinion) but it is not as general as the density function in the base package. It does not accomodate as wide a range of kernel types or bandwidth estimation methods. Both functions can produce a KDE using the normal kernel. That's okay, you can fix it shifting and scaling; add the smallest number such that the density is non-negative, then multiply the whole thing by a constant such that the area is unity. This is the easy way. The optimal$L_2$solution involves "water pouring"; finding the constant$c$such that$\left[\phi(x)-c\right]^+\$ integrates to unity.
• Notice that the question is rather on why the density function does not produce the "proper" density that integrates to 1 - rather then on how to fix it.