Simulating order statistics

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)
• The "observed empirical distribution" is not the density. One way to plot the observed density is to approximate it with a KDE and plot that. – whuber Oct 21 '16 at 19:29
• @Xi'an could you please clarify – mrb Oct 21 '16 at 20:44
• Compare your equation to the second line of fos_k: your equation omits the term you have coded as fx(x). But that's only an error of exposition; the error in your code lies in prob, which is not a density function. – whuber Oct 21 '16 at 21:40
• For what it's worth, here's the computation done in a single line in Mathematica: Show[Histogram[ParallelTable[Min[RandomVariate[NormalDistribution[0, 1], 10]], {10^6}], {0.1}, "PDF"], Plot[Evaluate[PDF[OrderDistribution[{NormalDistribution[0, 1], 10}, 1], z]], {z, -4, 1}]] On my machine, the runtime is less than 1 second. – heropup Oct 21 '16 at 22:08
• Thanks to everyone. @whuber is right. @whuber would you like to go on and reply to this question with an example of your R code? – mrb Oct 22 '16 at 1:54

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)
• Wondering should use the logarithm to define the corresponding distribution function as well? Since F(x) can be zero I would avoid the log. – mrb Oct 24 '16 at 13:03
• That's a good observation, but it happens not be be relevant because almost surely $F$ will be positive at any value you have randomly generated. Regardless, good software will deal with this correctly (provided you actually ask it to compute the log of $F$ directly rather than asking it to take the log of $0$!). For instance, the R procedures will return -Inf for the logarithm in that case; and that value, when exponentiated, will be exactly $0$. Here is an example: exp(punif(-2, log.p=TRUE)) – whuber Oct 24 '16 at 13:41