为什么密度曲线下面积的和总是大于1（R）？

Question

我找到了在R语言中计算密度曲线下面积和的代码。不幸的是，我不明白为什么总是会有一个额外的 ~"0.000976" 出现在这个面积上...

nb.data = 500000
y = rnorm(nb.data,10,2)

de = density(y)

require(zoo)
sum(diff(de$x[order(de$x)])*rollmean(de$y[order(de$x)],2))

[1] 1.000976

为什么会这样呢？

它应该等于1，对吗？

- M. Beausoleil

---

舍入误差？ - jmoon

有没有办法纠正这个问题？ - M. Beausoleil

和其他语言一样，我猜。我发现这个特别有帮助，但我不确定它在你的情况下是否适用。 - jmoon

请注意，这个值几乎正好偏移了 1/(2*length(de$y))。 - Ben Bolker

这必须与分布和积分算法有关。你总是使用正态分布，我想？并且隐含地使用相同的积分算法。 - Dirk Horsten

我希望我可以给 @BenBolker 的评论点赞不止一次。 - G5W

Answer 1

这就是微积分。使用更高的n（默认为512）可以获得更精确的结果。

set.seed(42)
de = density(rnorm(500000, 10, 2))
sum(diff(sort(de$x)) * 0.5 * (de$y[-1] + head(de$y, -1)))
#[1] 1.00098

set.seed(42)
de = density(rnorm(500000, 10, 2), n = 1000)
sum(diff(sort(de$x)) * 0.5 * (de$y[-1] + head(de$y, -1)))
#[1] 1.000491

set.seed(42)
de = density(rnorm(500000, 10, 2), n = 10000)
sum(diff(sort(de$x)) * 0.5 * (de$y[-1] + head(de$y, -1)))
#[1] 1.000031

set.seed(42)
de = density(rnorm(500000, 10, 2), n = 100000)
sum(diff(sort(de$x)) * 0.5 * (de$y[-1] + head(de$y, -1)))
#[1] 1.000004

set.seed(42)
de = density(rnorm(500000, 10, 2), n = 1000000)
sum(diff(sort(de$x)) * 0.5 * (de$y[-1] + head(de$y, -1)))
#[1] 1

- d.b

Answer 2

这种差异不仅仅是由于四舍五入误差或浮点数算术所致。实际上，您正在 density 计算的点之间进行线性插值，然后计算这个近似原始函数的面积（即使用梯形法则来积分曲线），这意味着您会高估曲线凹向上的区域的面积，并低估凹向下的区域。以下是维基百科文章中演示系统误差的示例图像：

---

^{Image by Intégration_num_trapèzes.svg: Scalerderivative work: Cdang (talk) - Intégration_num_trapèzes.svg, CC BY-SA 3.0, https://commons.wikimedia.org/w/index.php?curid=8541370}

---

由于正态分布具有更多凹向上的区域（即两端），因此总体估计值过高。如另一个答案中所提到的，使用更高的分辨率（即增加N）有助于最小化误差。您还可以使用不同的数值积分方法（例如辛普森法则）获得更好的结果。

然而，没有一种数值积分方法能给出精确的答案，即使有，density 的返回值本身也只是对真实分布的近似。（对于真实数据，真实分布是未知的。）

如果您想满足于看到已知密度函数积分为1，您可以在正态密度函数上使用integrate：

> integrate(dnorm, lower=-Inf, upper=Inf, mean=10, sd=2)
1 with absolute error < 4.9e-06

- Ryan C. Thompson

---

其实，我认为这会更具挑战性！有了积分，它甚至更好了。 - M. Beausoleil