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

Question

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

- Eric

Answer 1

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

假设我们有一些观测数据x，为了简单起见，我们假设这些数据是正态分布的：

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

我们进行密度估计（带有一些自定义，参见?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 ...

我们想要计算曲线右侧 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: 黎曼和

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

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：梯形公式

对 (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的答案

这个答案是合法的，但可能是在作弊。OP的问题始于密度估计，但是该答案完全绕过了它。如果允许这样做，为什么不简单地这样做呢？

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

- Zheyuan Li

---

我理解你在回复我的答案时的观点，我认为我们对问题的阅读方式不同。由于OP使用density()来获取密度估计值，因此他们将拥有原始数据。在这种情况下，我不认为我的答案是作弊。他们确实可以绕过对密度估计曲线进行积分。我仍在学习中，如果我错了，请纠正我。 - Robin

@李哲源 当我们必须计算phi的平方时该怎么办？我的意思是如何计算相同积分但是针对平方phi。 - Z B

Answer 2

基于 R 的经验累积分布函数 ecdf() 非常容易使用。以李哲源的例子为例...

#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

关于李哲源对我的回答的回应。 他们的答案假设你只有密度估计曲线。 而我的答案假设您拥有原始数据，这适用于OP的问题，因为他们使用density()来获得密度估计曲线。

- Robin