Simulating order statistics

Question

I am having a problem with the density of the first order statistic of a series of n random variables iid with common distribution (standard normal).

I am using Arnold's book as a reference for such a density function for the k'th order stat: $$ f_{X_{k:n}}(x) = \frac{n!}{k-1! (n-k)!} F(x)^{k-1} (1-F(x))^{n-k} f(x). $$

The I simulate data in R store the min and compare the obtained empirical distribution against the theoretical function. These curves diverge substantially.

What am I doing wrong? or missing?

# Distribution
Fx <- function(x) pnorm(x)
fx <- function(x) dnorm(x)


n <- 10
# ATTENTION: 6M simulation may take some time 
z <- replicate(1e6, min(rnorm(n)))

# Probabiltiy z > y

prob <- function(x) {
  sapply(x, function(u) mean(z>u))
}


# Density 
fos_k <- function(x, n, k) {
  fact <- factorial(n) / (factorial(k-1) * factorial(n-k)) 
  fact * fx(x) * Fx(x)^(k-1) *(1-Fx(x))^(n-k)
}

curve(fos_k(x, n=n, k=1), from=-1, to=1)
curve(prob, add=T, col=2, from=-1, to=1)

Answer 1

Although the problem is primarily with R code, it raises issues we would have to confront in any statistical computing environment. This reply focuses on those general issues.

---

A correct formula for the density of the $k^\text{th}$ smallest of $n$ independent identically distributed (iid) values from a distribution $F$ with density $f$ is presented in my answer at https://stats.stackexchange.com/a/225990/919. It is

$$ f_{[k]}(x) = \frac{n!}{(k-1)!(1)!(n-k)!} F(x)^{k-1} (1-F(x))^{n-k} f(x).\tag{1} $$

Because factorials and potentially high powers are combined, it is numerically better to compute its logarithm as

$$\eqalign{ \log f_{[k]}(x) = &\log n! - \log(k-1)! - \log(n-k)! + \\&(k-1)\log F(x) + (n-k) \log(1-F(x)) + \log f(x).\tag{2} }$$

Furthermore, because order statistics provide a way to peer far out into the tails of the distribution, where values grow close to $1$ and $0$, it is best to compute $1-F(x)$ directly rather than subtracting $F(x)$ from $1$.

With these caveats in mind, the other issue in this thread concerns graphing an empirical density. The commonest way to do so is with a histogram: it looks like a barplot in which the bar areas (not heights!) are proportional to the relative frequencies of the data.

To illustrate these two main points, I wrote a quick solution in R that simulates many iid samples, extracts specified order statistics from each, plots the histogram of each order statistic, and overplots the density function $(1)$ computed by exponentiating the logarithm $(2)$. Here is an example of its output applied to samples of size $10$ from a standard Normal distribution and of size $25$ from a Gamma$(3/2,5)$ distribution.

Clearly the agreement between the empirical results (the histogram bars) and the theoretical formula (the colored curves) is good.

This code applies the preceding suggestions about numerical computation by exploiting the log.p, log, and lower.tail arguments that are standard in the families of R distribution functions.

f <- function(n.sim, n, k=1:n, p=pnorm, d=dnorm, r=rnorm, name, ...) {
  if (missing(name)) name <- ""
  k <- sort(unique(k))[1 <= k & k <= n]

  # Perform the simulation.
  sim <- apply(matrix(r(n.sim*n, ...), nrow=n), 2, sort)[k, , drop=FALSE]

  # Plot the requested order statistics.
  for (i in 1:length(k)) {

    # Define a function to plot the density of an order statistic.
    dord <- function(x, k, n, ...) {
      z <- lfactorial(n) - lfactorial(k-1) - lfactorial(n-k) +
        (k-1)*p(x, log.p=TRUE, ...) + 
        (n-k)*p(x, log.p=TRUE, lower.tail=FALSE, ...) + 
        d(x, log=TRUE, ...)
      return(exp(z))
    }

    # Plot the empirical distribution.
    hist(sim[i, ], freq=FALSE, 
         xlab="Value",
         sub=paste(n.sim, "iterations with sample size", n),
         main=paste(name, "order statistic", k[i]))

    # Overplot the theoretical density.
    curve(dord(x, k[i], n, ...), add=TRUE, col=hsv(runif(1), 0.8, 0.7), lwd=2)
  }
}
#
# Set up to simulate and display the results.
#
par(mfrow=c(2,4))
n.sim <- 1e4
set.seed(17)

# Study Normal order statistics.
f(n.sim, 10, c(1,3,5,9), name="Normal(0,1)")

# Study Gamma order statistics.
# This illustrates how to use `f` for general distributions.
f(n.sim, 25, c(2, 4, 16, 24), name="Gamma(1.5,5)", 
  p=pgamma, d=dgamma, r=rgamma, shape=1.5, scale=5)