R中的顺序统计量?

5

我正在寻找一个函数来计算正态分布的顺序统计量。不是排名或指数,而是给定分布和样本大小的预期排名或指数(例如,预期最小值、中位数、最大值)。我熟悉解析解,但无法为正态分布解决/近似积分。是否有人知道在R中实现顺序统计量的软件包?


1
这是关于交叉验证的讨论:http://stats.stackexchange.com/questions/9001/approximate-order-statistics-for-normal-random-variables - MrFlick
2个回答

3
您正在寻求一个软件包。我不知道有没有这样的软件包,但我认为您可以在R中“解决/近似正态分布的积分”。实际上这很简单。
相关表达式是这篇论文中的方程(1):
其中ϕ是N [μ,σ]的概率密度函数,Φ是N [μ,σ]的累积分布函数。这些函数在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

谢谢!这当然是很有道理的。我还在考虑将被积函数简化为单一分布。直接抽样并没有想到。我忘记了计算机实际上是多么有用。您知道这种方法与从正态分布中进行大量重复抽样并平均样本最大值是否有所不同吗?我的直觉告诉我它们应该是相同的。但我不确定。 - user3738669

1

我缺少计算机代数系统的访问权限(作为一名生物学家,它们超出了我的日常工具箱),所以如果我无法手动解决被积函数,那么我只能自认倒霉。但我应该考虑研究一下开源计算机代数系统。谢谢。 - user3738669

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接