Posterior a consultar con expertos en la materia, programadores de deferentes ramas y maestros hemos dado con la respuesta al problema y la dejamos plasmada en esta publicación.
Solución:
Calcular áreas bajo una curva de estimación de densidad no es un trabajo difícil. Aquí hay un ejemplo reproducible.
Supongamos que tenemos algunos datos observados x que están, por simplicidad, distribuidas normalmente:
set.seed(0)
x <- rnorm(1000)
Realizamos una estimación de densidad (con alguna personalización, ver ?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 ...
Queremos calcular el área bajo la curva a la derecha de x = 1:
plot(d); abline(v = 1, col = 2)
Matemáticamente, esta es una integración numérica de la curva de densidad estimada en [1, Inf].
La curva de densidad estimada se almacena en formato discreto en d$x y 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`
Hay dos métodos para la integración numérica.
método 1: Suma de Riemann
El área bajo la curva de densidad estimada es:
C <- sum(yy) * dx ## sum(yy * dx)
# [1] 1.000976
Ya que Suma de Riemann es solo una aproximación, esto se desvía un poco de 1 (probabilidad total). Llamamos a esto C valor una "constante de normalización".
Integración numérica en [1, Inf] se puede aproximar por
p.unscaled <- sum(yy[xx >= 1]) * dx
# [1] 0.1691366
que debería ser escalado aún más por C para una estimación de probabilidad adecuada:
p.scaled <- p.unscaled / C
# [1] 0.1689718
Desde el true densidad de nuestra simulación x es saber, podemos comparar esta estimación con la true valor:
pnorm(x0, lower.tail = FALSE)
# [1] 0.1586553
que está bastante cerca.
método 2: regla trapezoidal
Hacemos una interpolación lineal de (xx, yy) y aplique la integración numérica en este interpolante lineal.
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
Con respecto a la respuesta de Robin
La respuesta es legítima pero probablemente engañosa. La pregunta de OP comienza con una estimación de densidad, pero la respuesta la omite por completo. Si esto está permitido, ¿por qué no simplemente hacer lo siguiente?
set.seed(0)
x <- rnorm(1000)
mean(x > 1)
#[1] 0.163
La función de distribución acumulativa empírica ecdf() en base R lo hace muy fácil. Usando el ejemplo de 李哲源...
#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
Con respecto a la respuesta de 李哲源 a mi respuesta. Su respuesta asume que solo tiene la curva de estimación de densidad. Mi respuesta supone que tiene los datos originales, que son aplicables a la pregunta del OP ya que usaron density() para obtener la curva de estimación de densidad.
Comentarios y puntuaciones del artículo
Si piensas que ha sido de ayuda este artículo, te agradeceríamos que lo compartas con otros desarrolladores y nos ayudes a dar difusión a nuestra información.