I have a density estimate (using density function) for my data learningTime (see figure below), and I need to find probability Pr(learningTime > c), i.e., the the area under density curve from a given number c (the red vertical line) to the end of curve. Any idea?
Computing areas under a density estimation curve is not a difficult job. Here is a reproducible example.
Suppose we have some observed data x that are, for simplicity, normally distributed:
set.seed(0)
x <- rnorm(1000)
We perform a density estimation (with some customization, see ?density):
d <- density.default(x, n = 512, cut = 3)
str(d)
# List of 7
# $ x : num [1:512] -3.91 -3.9 -3.88 -3.87 -3.85 ...
# $ y : num [1:512] 2.23e-05 2.74e-05 3.35e-05 4.07e-05 4.93e-05 ...
# ... truncated ...
We want to compute the area under the curve to the right of x = 1:
plot(d); abline(v = 1, col = 2)
Mathematically this is an numerical integration of the estimated density curve on [1, Inf].
The estimated density curve is stored in discrete format in d$x and d$y:
xx <- d$x ## 512 evenly spaced points on [min(x) - 3 * d$bw, max(x) + 3 * d$bw]
dx <- xx[2L] - xx[1L] ## spacing / bin size
yy <- d$y ## 512 density values for `xx`
There are two methods for the numerical integration.
method 1: Riemann Sum
The area under the estimated density curve is:
C <- sum(yy) * dx ## sum(yy * dx)
# [1] 1.000976
Since Riemann Sum is only an approximation, this deviates from 1 (total probability) a little bit. We call this C value a "normalizing constant".
Numerical integration on [1, Inf] can be approximated by
p.unscaled <- sum(yy[xx >= 1]) * dx
# [1] 0.1691366
which should be further scaled it by C for a proper probability estimation:
p.scaled <- p.unscaled / C
# [1] 0.1689718
Since the true density of our simulated x is know, we can compare this estimate with the true value:
pnorm(x0, lower.tail = FALSE)
# [1] 0.1586553
which is fairly close.
method 2: trapezoidal rule
We do a linear interpolation of (xx, yy) and apply numerical integration on this linear interpolant.
f <- approxfun(xx, yy)
C <- integrate(f, min(xx), max(xx))$value
p.unscaled <- integrate(f, 1, max(xx))$value
p.scaled <- p.unscaled / C
#[1] 0.1687369
Regarding Robin's answer
The answer is legitimate but probably cheating. OP's question starts with a density estimation but the answer bypasses it altogether. If this is allowed, why not simply do the following?
set.seed(0)
x <- rnorm(1000)
mean(x > 1)
#[1] 0.163
The Empirical Cumulative Distribution Function ecdf() in base R makes it very easy. Using 李哲源's example...
#Reproducible sample data
set.seed(0)
x <- rnorm(1000)
#Create empirical cumulative distribution function from sample data
d_fun <- ecdf (x)
#Assume a value for the "red vertical line"
x0 <- 1
#Area under curve less than, equal to x0
d_fun(x0)
# [1] 0.837
#Area under curve greater than x0
1 - d_fun(x0)
# [1] 0.163
Regarding 李哲源's response to my answer. Their answer assumes you only have the density estimate curve. My answer assumes you have the original data, which is applicable to the the OP's question since they used density() to get the density estimate curve.
