Calculate tail probabilities from density() call in R This question concerns how to implement the following problem in R. 
x = rnorm(1000)
hist(x,freq=FALSE)
lines(density(x))

How would you calculate the upper (or lower) tail probability for a given cutoff (e.g. +1) given the density estimate above? NOTE: the following solution isn't good enough. I need the calculation based on the smoothed density curve, i.e. an integral of the curve, not the empirical histogram. 
sum(x>1)/length(x)

Also please do not suggest the use any of the standard pnorm functions because they are only correct if the underlying distribution is correctly specified. Thanks!
 A: The density function returns an object with various properties. You can access the $x$ and $y$ values using density(x)$x and density(x)$y. So you can do it like this:
set.seed(100)
x <- rnorm(1000)
d <- density(x)
x0 <- 1
idx <- which(abs(d$x-x0)==min(abs(d$x-x0)))
approx.tail.prob <- sum(d$y[idx:length(d$x)] * diff(d$x)[1])

This is just an approximation based on a Riemann sum. You could get a better approximation using another numerical technique, such as the Trapezium Rule or Simpson's Rule. But once you know how to get at density(x)$x and density(x)$y, it's straightforward to work out how to do these.
You could even use the R integrate function, maybe like this:
f <- function(x0) d$y[which(abs(d$x-x0) == min(abs(d$x-x0)))[1]]

and then:
integrate(Vectorize(f), 1, max(d$x))

A: By definition, given a "bandwidth" $h$ and a kernel density $k,$ the KDE of a data vector $\mathbf{x} = (x_1, x_2, \ldots, x_n)$ is
$$f(x; \mathbf{x}, h, k) = \frac{1}{nh}\sum_{i=1}^n k\left(\frac{x - x_i}{h}\right).$$
Consequently the distribution function (left tail probability function) is its integral,

$$F(x; \mathbf{x}, h, K) = \frac{1}{nh}\sum_{i=1}^n \int_{-\infty}^x k\left(\frac{t - x_i}{h}\right)\,\mathrm{d}t = \frac{1}{n}\sum_{i=1}^n K(t-x_i; h)$$

where $K$ is the integral of $k,$
$$K(x) = \frac{1}{h}\int_{-\infty}^x k\left(\frac{t}{h}\right)\,\mathrm{d}t.$$
By default, density uses a Gaussian (Normal) kernel: that is, $k$ is implemented as dnorm and $K$ as pnorm.  This leads to an extremely compact and precise solution:
pkde <- Vectorize(function(x, data, bw) mean(pnorm(x, data, bw)), "x")

The arguments of  pkde are the point x of evaluation, the data vector data $=\mathbf{x},$ and the bandwidth bw $=h.$  For instance, the right tail probability at the value $x=1$ in the example is obtained by storing the result of density in order to fetch bw:
obj <- density(x)
1 - pkde(1, x, obj$bw)

The answer (for the question's sample data) will be close to $0.16,$ depending on the specific values that were generated.

Some comments about the use of pnorm may be of interest. pnorm appears here because it describes your kernel.  Any other solution will be equivalent to this one.  This one is the most precise of any possible solution because it does not require interpolation.  It will be inferior in terms of computational complexity when $n$ is large, because it requires one evaluation of pnorm for each data point.  It will be greatly superior to any conceivable alternative when n is small.  Indeed, it doesn't even require the KDE to be computed: it only needs you to determine what bandwidth you want to use.
Here is a more interesting example, showing how pkde can be used to plot the entire cumulative distribution of the KDE for a decidedly non-Normal dataset.
set.seed(17)
X <- c(rexp(500), rgamma(500, 30))
obj <- density(X)
curve(pkde(x, X, obj$bw), min(X), max(X), lwd=2, main="Left tail kernel probability")
hist(X, breaks=100)


A: I would take the same approach as @Flounderer, but exploit another feature of R's density() function; namely the from and to arguments, which restrict the density estimation to the region enclosed by the two arguments. This results in the same density estimates as running the function without from and/or to, but by restricting the range of the density estimate to the region of interest, we focus all of the n evaluation points on the region of interest.
set.seed(1)
x <-rnorm(1000)
hist(x,freq=FALSE)
lines( dens <- density(x) )
lines( dens2 <- density(x, from = 1, n = 1024), col = "red", lwd = 2)

This produces

The red line is to illustrate that the density estimates in dens and dens2 are the same for the region of interest.
Then you can follow the approach @Flounderer used to evaluate the tail probability:
> with(dens2, sum(y * diff(x)[1]))
[1] 0.1680759

The advantage of this approach is to expend the n observations at which density() evaluates the KDE all on the region of interest. The larger n the higher the resolution that you have in evaluating the tail probability.
Note from ?density that given the FFT used in the implementation, having n as a multiple of 2 is advantageous.
