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)
fos_k
: your equation omits the term you have coded asfx(x)
. But that's only an error of exposition; the error in your code lies inprob
, which is not a density function. $\endgroup$ – whuber♦ Oct 21 '16 at 21:40Show[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. $\endgroup$ – heropup Oct 21 '16 at 22:08R
code? $\endgroup$ – mrb Oct 22 '16 at 1:54