r - R中的订单统计？

Question

我正在寻找一个函数来计算正态分布的订单统计信息。不是排名或指数，而是给定排名或指数的预期，给定分布和样本大小(即预期最小值、中值、最大值)。我熟悉解析解，但无法求解/近似正态分布的积分。有谁知道在 R 中实现订单统计的包？

Answer 1

你要一个包裹。我不知道，但我不明白为什么你不能在 R 中“求解/近似正态分布的积分”。它实际上很简单。

相关表达式为this paper中的方程(1) :



其中 ϕ 是 N[μ, σ] 的 PDF，Φ 是 N[μ, σ] 的 CDF。这些函数作为 dnorm(...) 内置在 R 中和 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

编辑 对 OP 评论的回应。

是的，从大量随机抽取中平均样本最大值将近似于 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