我正在寻找一个函数来计算正态分布的顺序统计量。不是排名或指数,而是给定分布和样本大小的预期排名或指数(例如,预期最小值、中位数、最大值)。我熟悉解析解,但无法为正态分布解决/近似积分。是否有人知道在R中实现顺序统计量的软件包?
dnorm(...)
和pnorm(...)
内置函数使用。f <- function(x, mu=0, sigma=1) dnorm(x, mean=mu, sd=sigma)
F <- function(x, mu=0, sigma=1) pnorm(x, mean=mu, sd=sigma, lower.tail=FALSE)
integrand <- function(x,r,n,mu=0, sigma=1) {
x * (1 - F(x, mu, sigma))^(r-1) * F(x, mu, sigma)^(n-r) * f(x, mu, sigma)
}
E <- function(r,n, mu=0, sigma=1) {
(1/beta(r,n-r+1)) * integrate(integrand,-Inf,Inf, r, n, mu, sigma)$value
}
E(1,1000) # expected value of the minimum
# [1] -3.241436
E(1000,1000) # expected value of the maximum
# [1] 3.241436
E(500.5,1000) # expected value of the median
# [1] -6.499977e-18
编辑 对评论的回应。
是的,从大量随机抽取的样本最大值平均将近似于E(n,n)
。然而,有两个不同之处。首先,答案将是近似值,而上述结果是精确的(至少对于数值积分的精度而言)。其次,第一种方法运行速度约快30倍。
E.max <- function(n) mean(sapply(1:100,function(i)max(rnorm(n))))
E.max(1000)
# [1] 3.267614
library(microbenchmark)
microbenchmark(E(1000,1000),E.max(1000))
# Unit: milliseconds
# expr min lq median uq max neval
# E(1000, 1000) 1.027536 1.169674 1.333428 1.50429 1.905828 100
# E.max(1000) 23.889773 28.882058 32.642485 37.37952 39.830501 100
虽然不是直接回答,但您可以使用计算机代数系统轻松计算闭合形式密度。然后对该密度进行采样以获得最小值/最大值/中位数的估计。
请参见:http://www4.stat.ncsu.edu/~hzhang/st522/08Chapter5_order.pdf