Best estimator of the mean of a normal distribution based only on box-plot statistics

Question

Suppose $X_1,\ldots,X_n\sim\operatorname N(\mu,\sigma^2)$ and you can observe only the sample size $n,$ the two extreme values, and the first, second, and third quantiles of the sample. Among unbiased estimators of $\mu$ based only on those five statistics, which has the smallest variance?

Answer 1

An exact answer will be difficult, so first I will look at asymptotic theory. Answers from that could be tested by simulation, comparing it to a maximum likelihood estimator computed by maximizing an exact likelihood based on the joint distribution of these order statistics.

One problem is that the extreme order statistics maximum and minimum do not have normal asymptotic distributions, so the form of the joint asymptotic distribution of the extrema and the central order statistics is not clear. So for the moment I will look at the simpler problem of using only the median and the first and third quantiles.

From the book "A first course in order statistics" by Arnold at.al, we have the following theorem 8.5.2: Let $i_r =[n p_r]+1, 1 \le r \le k$, $0 < p_1 < \dotsm < p_k < 1$ such that $f(F^{-1}(p_r)$ is finite and positive. Then the joint distribution of $\sqrt{n} \left( X_{i_r:n} - F^{-1}(p_r) \right)$ is asymptoticaly $k$-dim normal with mean zero and covariance matrix $\Sigma^A = (\sigma_{rs}^A)$ with $$ \sigma_{rs}^A = \frac{p_r (1-p_s)}{f(F^{-1}(p_r)) f(F^{-1}(p_s))} $$ Let the vector of the three observations used be $Q=(Q_1, M, Q_3)$. An linear unbiased estimator of the mean based on $Q$ must have the form (based on symmetry) $$ a M + \frac{1-a}{2} (Q_1 + Q_3) $$ The solution cannot depend on $(\mu, \sigma)$ so we can for the calculations consider the standard normal. Then, with a little help from R we find the covariance matrix as in the following:

sigmaA <- function(p) {
  p <- sort(p)
  p[1]*(1-p[2]) / ( dnorm(qnorm(p[1]))*dnorm(qnorm(p[2]))  )
}

SIGMA <- matrix(0, 3, 3)

ps <- 1:3/4

for (i in 1:3) for (j in 1:3) {
  SIGMA[i, j] <- sigmaA( c(ps[i], ps[j]))
}

giving

SIGMA
          [,1]      [,2]      [,3]
[1,] 1.8567675 0.9860026 0.6189225
[2,] 0.9860026 1.5707963 0.9860026
[3,] 0.6189225 0.9860026 1.8567675

Now, we can minimize (jumping some details ...) the variance by minimizing the quadratic form $A^T \Sigma A$ with $A=((1-a)/2, a, (1-a)/2)$ which is a simple problem.

(it is late night now, so I leave it there to continue after some sleep)

Answer 2

My thoughts so far:

1. We can work with the standard normal; the general case will follow immediately from that.

2. (Please note that my original argument here was wrong; I neglected some terms that I failed to recognize were functions of $n$) I think the median and the midrange should be asymptotically uncorrelated as should the midhinge and the midrange, but at small sample sizes are clearly dependent. The midhinge and the median will still be substantively correlated asymptotically.

   I think the asymptotic correlation matrix of the three quartiles $Q = (Q_1,Q_2,Q_3)^\top$ should be

   $$C= \left[ \begin{array}{ccc} 1 & \frac{1}{\sqrt{3}} & \frac13 \\ \frac{1}{\sqrt{3}} & 1 & \frac{1}{\sqrt{3}} \\ \frac13 & \frac{1}{\sqrt{3}} & 1 \end{array} \right]$$

3. By symmetry, at a fixed sample size we should be looking at a function of the median, the midhinge and the midrange $-$ here I am ignoring, as Wikipedia does in its article on the midhinge, distinction between hinges and sample quartiles. Asymptotically the median and midhinge will combine linearly (due to their joint normality)^[1]. I'll denote these by $M$ (rather than the $Q_2$ used above), $\bar{H}=\frac{Q_1+Q_3}{2}$ and $\bar{R}=\frac{X_{(1)}+X_{(n)}}{2}$ respectively. $\bar{R}$ will not be normally distributed, but will be symmetric. (Indeed, asymptotically, for samples from normal populations, the midrange should converge to a logistic with mean $\mu$, but it does seem to come in fairly slowly. Simulations suggest that around $n=2001$ the approximation looks fairly reasonable and by $n=10001$ it seems to be an excellent approximation - the ecdf is extremely close to a logistic - though a logistic Q-Q plot does still have a small degree of wiggle even then, though that's not surprising)

   The asymptotic variance of $(1-q)M+q\bar{H}=u^\top Q$ (where $u=(\frac{q}{2}, (1-q), \frac{q}{2})^\top$) should be $u^\top Vu$, where $V$ is the variance-covariance matrix of $Q$. Hence we should be able to write down the $q$ that will minimize that asymptotic variance.

4. As sample size increases the contribution of the midrange should decrease, since the variance of the midrange decreases more slowly as a function of n than the other order statistics. Indeed, from the data in Tippett (1925)[2], and via simulation, the expected value of both the sample maximum and the range at large sample sizes seems to be progressing close to linearly in $\log(n)$, and similarly with the precision of the maximum (and by symmetry, the minimum); asymptotically the variance of the midrange should be proportional to the variance of the maximum.

   Even though it likely won't be quite optimal I'll focus on considering linear combinations of all three statistics, since the contribution of the midrange will be less less and less critical in large samples. In short, the estimator $(1-q-r)M+q\bar{H}+r\bar{R}$ for suitably chosen $q$ and $r$ (each is constant for a fixed $n$, but each changes as a function of $n$) should be quite close to optimal at all but perhaps quite small $n$; this simple linear functional of the statistics should perform quite well in general.

   At $n=5$, the contribution of all five order statistics should be the same -- it's just the average of the five observations. So $q=r=0.4$ there.

5. As $n$ increases, we should expect that $q$ asymptotes to a constant, but given the variance behavior of the midrange, we should expect that $r$ asymptotes down to $0$.

   We could compute asymptotic weights on $M$ and $\bar{H}$ using the asymptotic variance-covariance matrix of the corresponding order statistics. This should provide a useful reasonableness check on any estimate for $q(n)$ as $n$ grows large.

6. One approach to getting the coefficients for this not-quite-efficient estimator would be to use simulation to choose approximate $q$ and $r$ at a range of finite $n$ above $n=5$ $-$ given the estimator will be unbiased, we're presumably just trying to minimize variance, so we're looking for reasonably a smooth additive functions $q(n)$ and $r(n)$ that minimize variance. We know that $q$ should increase from $0.4$ to the asymptotic value and $r$ should (slowly) decrease from $0.4$ downward toward $0$.

   However, simulation shouldn't be required; since we're minimizing the variance, and the estimator is linear in a set of order statistics, we should be able to proceed to get coefficients that minimize variance directly from the variance-covariance matrix of the relevant order statistics.

edit: It should be noted that the functions of $n$ will not be smooth unless one focuses separately on the $n$ being $0,1,2,3$ $\mod 4$.

I'll have to come back to this, so I'll post what I have for now but I believe this is a good outline of a reasonable way to proceed. I'll be surprised if someone hasn't already performed this exercise or one very close to it, to be honest. I'll see what I can turn up.

[Edit: Jacob mentions Luo's paper in a comment under kjetil's answer, which suggests an estimator of exactly the form suggested here -- a linear-weighted average of median, midhinge and midrange, where the weights are functions of $n$. As Henry points out in comments, that paper is for estimating a sample mean rather than the population mean, but it's interesting to see the same form of estimator come up.]

---

[1]: To clarify; at a given $n$, the $q$ is purely a constant. We can write down its form explicitly; if $u=(q/2,1-q,q/2)^\top$ then we seek to minimize $\operatorname{var}(u'V)=u'Vu$ where $V=\operatorname{var}(Q)$, the covariance matrix of the "quartiles" (or we can write it directly in terms of $(q,1-q)$ and the covariance matrix of $(M,\bar{H})$; it makes one step easier and another step harder). Either way, it's a quadratic in $q$ bounded to $[0,1]$ (the more general case is a simple quadratic programming problem with linear constraints). We can work out the elements of $V$ at the given $n$ and find the $q$ to minimize the quadratic, and then check it's between the endpoints (checking the variance at the $q=0$ and $q=1$ values themselves is simple, since we already compute those along the way). The elements of $V$ can be found from the joint distribution of the involved order statistics. Asymptotically in the normal case, when the quartiles are all individual order statistics this involves products of terms like $\phi(\Phi^{-1}(p))$ for specific $p$ with $0<p<1$. In small samples there's integrals to compute (or tables, or approximations might be used). If you want to write $q$ as a function of $n$, the exact form will be complicated, but suitable approximations should be possible.

[2]: L. H. C. Tippett (1925). "On the Extreme Individuals and the Range of Samples Taken from a Normal Population". Biometrika 17 (3/4): 364–387

NB, this paper is some years out of copyright. It can be located legally for free online, don't let rapacious publishers* charge you large amounts (like 39 USD for 24 hours access, I kid you not) to look at something they know perfectly well is in the public domain. $\,$ $\,$ * I really don't see how tamer words do justice to their behavior on this.

Answer 3

The maximum likelihood solution I posted long ago at Estimating a normal distribution from three order statistics works well. The only questions concern (A) how well and (B) whether its estimate of the mean is unbiased.

There are some surprises.

I won't prove (B), but intuitively the symmetric positions of the five order statistics used for the 5-letter summary (the boxplot statistics) will produce an estimate symmetrically distributed around the correct parameter, which therefore is unbiased.

This leaves (A). Let's compare the maximum likelihood estimates of the parameters $\mu$ and $\sigma$ to those that would be obtained using all the data.

A quick simulation shows us how things are shaking out. I chose a sample size of just $n=8,$ because maximum likelihood is most likely to perform poorly for the smallest samples. The 5-letter summary I chose uses the order statistics at positions $1,3,4,6,8,$ which isn't even symmetric. (The median is usually at position $4.5,$ interpreted as the arithmetic mean of the order statistics at positions $4$ and $5$. For simplicity, I don't attempt to use quantiles that interpolate among order statistics, so this solution isn't even going to be optimal.)

Here are the results of applying this 5-letter summary estimator to $500$ independent samples of a standard Normal distribution (for parameters $\mu=0$ and $\sigma=1,$ with no loss of generality because the estimator is location and scale covariant), yielding estimates $\hat\mu$ and $\hat\sigma,$ along with the results of applying the usual maximum likelihood estimates $\bar x$ (the mean) and $\operatorname{SD}(x)$ (computed as the root mean squared deviation from the mean, not with any bias correction).

The surprises are

1. The distribution of $\hat \mu$ scarcely differs from that of $\bar x.$ Compare the histograms in the left column and see the summary below.

2. The distribution of $\hat\sigma$ is narrower, less skewed, and less biased than the distribution of $\operatorname{SD}(x)$! (Compare the histograms in the right column and notice the slightly wider range of the x-axis at the bottom.)

3. Because I used asymmetric order statistics, the estimate of the mean is (slightly) biased (high in this case). (To obtain an unbiased estimate, one could take two estimates where the median is considered as either at position $4$ or position $5$ and average the two estimates.)

Here are the averages and standard errors from these $500$ samples:

             Mean     SE
mu.hat     0.0046 0.0130
x.bar     -0.0200 0.0130
sigma.hat  0.9500 0.0098
SD         0.8900 0.0180

To demonstrate the lack of bias in $\hat\mu$ when the order statistics are symmetric (that is, they are at positions $1,$ $k,$ $(n+1)/2$, $n+1-k,$ and $n$ -- and therefore $n$ must be odd), I re-ran the simulation for a sample size of $n=9$ (for order statistics $1,3,5,7,$ and $9$) and created $5000$ samples for greater precision. The histograms are about the same as before, but here are the summary statistics:

            Mean     SE
mu.hat    0.0062 0.0049
x.bar     0.0071 0.0048
sigma.hat 0.9100 0.0034
SD        0.8900 0.0063

Now neither of the estimates of the mean (first two rows) differs significantly from zero (both are well within two SEs of zero) and $\hat\sigma$ is still superior to $\operatorname{SD}(x),$ although not by as much as before.

Finally, here is a summary for $n=79,$ again with $5000$ samples (using order statistics $1,20,40,60,$ and $79$).

                Mean      SE
mu.hat     0.000160 0.00168
x.bar     -0.000624 0.00160
sigma.hat  0.978000 0.00135
SD         0.988000 0.00226

TANSTAAFL, and reality intrudes: $\bar x$ must be superior to $\hat\mu,$ and so it is, as evidenced by its slightly smaller standard error. But $\hat\mu$ remains unbiased. Also, there is no detectable difference now between $\hat\sigma$ and $\operatorname{SD}(x),$ both of which are closely approximating $\sigma = 1.$

Because the performance of $\hat\mu$ is practically the same as that of $\bar x,$ which is a minimum variance unbiased estimator based on all the data, in most applications it would be pointless to try to improve on $\hat\mu.$

It would still be interesting to compare these estimators to the linear-combination-of-order-statistics estimators proposed in other answers, because those are attractive approaches and can admit simpler computation.

---

Here are details of the simulations (set up to perform the first one). At the end, the four estimates are contained in the four rows of sim. Plotting their histograms is routine and omitted for brevity.

#
# Negative log likelihood.
#
Lambda <- function(theta, values, orders) {
  mu <- theta[1]
  sigma <- exp(theta[2])
  f <- dnorm(values, mu, sigma, log = TRUE)
  F <- log(diff(pnorm(c(-Inf, values, Inf), mu, sigma)))
  -(sum(f) + sum(F * (diff(c(0, orders, n+1))-1)))
}
#
# A simulation study.
#
n <- 8
q <- c(`X-` = 0, `H-` = 1/4, M = 1/2, `H+` = 3/4, `X+` = 1)
orders <- round(quantile(seq_len(n), q))
mu <- 0; sigma <- 1
set.seed(17)
sim <- replicate(5e2, {
  x <- rnorm(n, mu, sigma)  # The full dataset
  nletter <- quantile(x, q) # Its 5-letter summary
  #
  # Compute the estimate.
  #
  theta.start <- c(mu = mean(range(nletter)), 
                   log.sigma = log(diff(range(nletter)) / 6))
  fit <- nlm(Lambda, theta.start, values = nletter, orders = orders) 
  if (fit$code > 2) warning("nlm did not succeed.")
  theta.hat <- setNames(fit$estimate, names(theta.start))
  c(theta.hat, xbar = mean(x), sd = mean((x - mean(x))^2))
})
sim[2, ] <- exp(sim[2, ]) # (The log of sigma was estimated.)
rownames(sim) <- c("mu.hat", "sigma.hat", "x.bar", "SD")