Estimate normal distribution from small sample with rankings

Question

If I have n samples from an $N(\mu,\sigma)$ distribution, how can I estimate the distribution from a subset of m of these n samples where the order within those n is known. For example, if n = 100 and m is 3, I might know that one sample is 10 and 3rd highest of the 100, one is 5 and 42nd highest, and the last is 2 and 74th highest.

Answer 1

Assuming the $n$ values (forming a sample $X$, say) are independently and randomly obtained and the ranks $m$ are not somehow based on those values, one can always apply maximum likelihood when $n$ is "not small." Another procedure is based on the "regression on order statistics" (ROS) method popular among people dealing with censored data.

Maximum Likelihood

The elements of the subset (let's call it $Y$) corresponding to $m$ are order statistics for the sample. The $k^\text{th}$ highest has rank $n+1-k$. For instance, in a subset of a sample of $n=100$ containing the third, 42nd, and 74th highest, the ranks are $98$, $59$, and $27$, respectively. To simplify some formulas, I will also include $-\infty$ at rank $0$ and $+\infty$ at rank $n+1$.

The joint distribution for the order statistics can obtained in terms of the parameters $\mu$ and $\sigma$ (using, say, the argument at https://stats.stackexchange.com/a/78559). Briefly, letting $F_{\mu,\sigma}$ be the associated Normal cumulative distribution and $f_{\mu,\sigma}$ its density, the density of the joint distribution is a product of three kinds of terms:

1. Each value $y \in Y$ contributes $f_{\mu,\sigma}(y)$.

2. Each gap between one order statistic $y$, with rank $r \ge 0$, and the next one $z$, with rank $s$ ($r \lt s \le n+1$), contributes $(F_{\mu,\sigma}(z) - F_{\mu,\sigma}(y))^{s-r-1}$.

3. The normalizing constant is the multinomial coefficient $\binom{n}{1,1,\ldots,1;k_1, k_2, \ldots, k_m}$ where the $k_i$ are the exponents which appear in (2). This term will not matter for maximum likelihood estimation (MLE) because it does not depend on the parameters.

As usual, one maximizes the log likelihood. The full machinery is described elsewhere in many places. Working code--good enough to obtain the MLE--is provided below.

ROS

Associate with each rank $r$ its Weibull plotting point $p_r$ where

$$p_r = (r - 0.3175) / (n + 0.365)$$

for $2 \le r \le n-1$. The two end plotting points are $2^{-1/n}$ and $1-2^{-1/n}$ for $r=1$ and $r=n$, respectively.

A "Normal probability plot" for $X$ pairs $F_{0,1}^{-1}(p_r)$ with the $r^\text{th}$ smallest value $x_{[r]}\in X$. Such a plot will tend to be linear when $X$ comes from a Normal distribution. Its intercept estimates $\mu$ and its slope estimates $\sigma$ Thus, the subset of pairs drawn from $Y$ should also be approximately linear and can be used in the same way to estimate the parameters by fitting a straight line. Although various sophisticated forms of fitting could be tried, it is routine (especially when $m$ is large) just to use ordinarily least squares (even though the order statistics are necessarily increasing, their joint distribution can exhibit strong correlations, and the variances are not constant).

Results

Simulations suggest that MLE can estimate $\mu$ accurately but is biased low for $\sigma$. The bias can be pretty large for small $m$ and $n$, around $10\%$ low for $m=3$ and $n=10$ even when the ranks are well positioned.

Simultaneous simulations suggest ROS may be a little less accurate and is biased high for $\sigma$. In the following, the values in $X$ were generated randomly from a Normal distribution with $\mu=0$ and $\sigma=1$. The three points in $Y$, used for the fitting, are shown as solid red dots. The gray line shows the fitted distribution. The gray dots are the values in $X$. For a good fit, they ought to follow the line closely. They do, except for three values in the bottom left (which could not possibly have been detected by $Y$ and are merely the result of chance variation in the sampling anyway).

Discussion

For either procedure to have a chance of working well, there should be at least one rank in the lower half of $X$ and at least another in the upper half. Otherwise, all values in $Y$ reflect only one tail of the distribution of $X$. It is difficult to pin down the parameters from one tail only. On the other hand, if the objective is to obtain a fit that works very well in just one tail, this behavior could be an asset!

Both procedures can give excellent results with very small values of $m$, provided the ranks are suitably situated. I recall one study (which I cannot quote specifically) suggesting that just two ranks near the $16^\text{th}$ and $84^\text{th}$ percentiles can work well.

Both procedures are robust in the sense that they enjoy bounded sensitivity to changes in the extremes of $X$--specifically, changes in those values more extreme than the smallest and largest elements of $Y$.

MLE is probably not a good procedure to use when $n$ is small: its bias can be large. Consider using a parametric bootstrap to estimate and correct that bias.

My own experience is that ROS becomes more and more unstable as the elements of $Y$ crowd toward one end of $X$ or the other. This, unfortunately, is exactly the situation where ROS is most popular: a majority of the data may be censored by detection limits (or saturation limits) in measuring devices, leaving quantitative information only for the most extreme values at one end. I would expect MLE to suffer from the same instability.

---

Code

The following is some of the R code I wrote to test both procedures. It implements ROS and MLE (without computing standard errors of estimate, though). It contains code to draw the figure and code to run Monte-Carlo simulations of both procedures.

#
# Negative log likelihood for parameters `theta`, based on a subset `y` from
# a sample of `n` iid values, associated with ranks in `m`.
#
Lambda <- function(theta, y, m, n) {
  mu <- theta[1]; sigma <- theta[2]
  powers <- diff(c(0, m, n+1)) - 1
  gaps <- diff(c(0, pnorm(y, mu, sigma), 1))
  z <- dnorm(y, mu, sigma)
  #a <- lgamma(n+1) - sum(lgamma(powers+1)) # Log multinomial coefficient
  b <- sum(log(z)) + sum(powers*log(gaps))
  return (-b)
}
#
# Return the Weibull probability plotting points for `n` values.
#
pp <- function(n) {
  mid <- c()
  if (n >= 3) mid <- (2:(n-1) - 0.3175)/(n + 0.365)
  u <- (1/2)^(1/n)
  c(1 - u, mid, u)
}
#
# ROS fit.
#
ROS <- function(y, m, n) {
  sigma <- qnorm(pp(n)[m])
  return (coef(lm(y ~ sigma)))
}
#
# MLE fit (estimates only).
#
MLE <- function(y, m, n) {
  theta <- ROS(y, m, n)      # Use ROS for an initial estimate
  fit <- optim(theta, Lambda, y=y, m=m, n=n)
  z <- c(fit$par, theta)
  names(z) <- c("mu.MLE", "sigma.MLE", "mu.ROS", "sigma.ROS")
  return (z)                 # Return MLE *and* ROS estimates
}
#
# Generate a dataset from a known Normal distribution.
#
n <- 100
m <- sort(c(98, 59, 27)) # The ranks to select
set.seed(17)             # Allow reproducible results
x <- sort(rnorm(n))      # Sort the data once and for all
y <- x[m]                # All fitting is based only on `y`, `m`, and `n`
#
# Fit using both methods.
#
fit <- MLE(y, m, n)
mu.hat <- fit[1]; sigma.hat <- fit[2]
#
# Plot the MLE fit.
#
pp.x <- qnorm(pp(n), mu.hat, sigma.hat)
pp.y <- pp.x[m]
plot(pp.x, x, col="Gray", xlab="Fit", ylab="Values",
      main="Simulated Data, n=100")
abline(c(0,1), col="Black", lty=3)
points(pp.y, y, pch=16, col="#d03020", cex=1.25)
#
# Peform a Monte-Carlo assessment of the bias and variance of both.
#
system.time({            # Takes about 4 seconds per 1000 iterations
  sim <- replicate(1e3, {
    x <- sort(rnorm(n))
    y <- x[m]
    MLE(y, m, n)
  })
})

rowMeans(sim)            # Average estimates of the parameters
apply(sim, 1, sd)        # Standard deviations of those estimates