Order statistics in R?

I am looking for a function to calculate order statistics from regular distributions. Not a ranking or index, but the expectation of a given rank or index, given the distribution and size of the sample (i.e., expected minimum, average, max). I am familiar with the analytical solution, but could not solve / approximate the integrals for the normal distribution. Does anyone know of a package that implements order statistics in R?

+4
source share
2 answers

You are requesting a package. I don’t know anything, but I don’t understand why you cannot “solve / approximate the integrals for the normal distribution” in R. This is actually quite simple.

(1) :

φ - PDF N [μ, σ], Φ - CDF 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

OP.

, max 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
+3

Source: https://habr.com/ru/post/1544523/


All Articles