r - 计算密度估计曲线下的面积，即概率

Question

我对数据learningTime有一个密度估计(使用密度函数)(见下图)，我需要找到概率Pr(learningTime > c)，即密度曲线下从给定数字c(红色垂直线)到曲线末端的面积。有什么想法吗？

Answer 1

计算密度估计曲线下的面积并不是一件困难的工作。这是一个可重现的示例。

假设我们有一些观察到的数据x，为简单起见，这些数据呈正态分布:

set.seed(0)
x <- rnorm(1000)

我们执行密度估计(进行一些自定义，请参阅？密度):

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 ...

我们要计算 x = 1 右侧曲线下的面积:

plot(d); abline(v = 1, col = 2)

从数学上讲，这是[1, Inf]上估计密度曲线的数值积分。

估计的密度曲线以离散格式存储在d$x和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`

数值积分有两种方法。

方法1:Riemann Sum

估计密度曲线下的面积为:

C <- sum(yy) * dx  ## sum(yy * dx)
# [1] 1.000976

由于黎曼和只是一个近似值，因此它与 1(总概率)有一点偏差。我们将此 C 值称为“归一化常数”。

[1, Inf] 上的数值积分可以近似为

p.unscaled <- sum(yy[xx >= 1]) * dx
# [1] 0.1691366

应该通过C进一步缩放它以获得正确的概率估计:

p.scaled <- p.unscaled / C
# [1] 0.1689718

由于我们模拟的 x 的真实密度是已知的，我们可以将此估计值与真实值进行比较:

pnorm(x0, lower.tail = FALSE)
# [1] 0.1586553

相当接近。

方法2:trapezoidal rule

我们对(xx, yy)进行线性插值，并在此线性插值上应用数值积分。

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

<小时/>

关于 Robin's answer

答案是合法的，但可能是作弊。 OP的问题从密度估计开始，但答案完全绕过了它。如果允许这样做，为什么不简单地执行以下操作？

set.seed(0)
x <- rnorm(1000)
mean(x > 1)
#[1] 0.163