我找到了用来计算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
Run Code Online (Sandbox Code Playgroud)
为什么会这样?
它应该等于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
Run Code Online (Sandbox Code Playgroud)
这种差异不仅仅是由于舍入误差或浮点运算引起的.您有效地在计算的点之间进行线性插值density,然后在此近似下计算原始函数的面积(即,您使用陷阱规则对曲线进行积分),这意味着您高估了曲线区域中的面积.凹陷并在向下凹陷的区域中低估它.以下是维基百科文章中展示系统错误的示例图片:
![]()
图片来自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)有助于最小化误差.使用不同的数值积分方法(例如Simpson规则),您也可能获得更好的结果.
然而,没有数值积分方法可以给你一个确切的答案,即使有,但返回值density只是实际分布的近似值.(对于真实数据,真正的分布是未知的.)
如果您想要的是将已知密度函数积分为1的满足感,则可以使用integrate常规密度函数:
> integrate(dnorm, lower=-Inf, upper=Inf, mean=10, sd=2)
1 with absolute error < 4.9e-06
Run Code Online (Sandbox Code Playgroud)