Estimating a normal distribution from three order statistics

I am interested in predicting a normal distribution, but not sure if this is possible.

I do not have information on the mean or standard deviation. However, I know the range of values, let's say from 0 to 10, and I know the sample size, let's say 1000, and I also know that value 8.1 is the 50th highest value.

Is it possible in any way from this information to infer more about what the distribution looks like?

• "However, I know the range of values" ... Do you mean by this that you know the highest and lowest values in the sample? – Silverfish Apr 14 '16 at 21:54
• Yes, sorry that is what I meant, they range from 0 to 10. So 0 is lowest and 10 is higher. – AdrianP. Apr 14 '16 at 21:55
• It sounds to me like you know three "order statistics" for your sample (the 1st, 50th and 1000th/last) and want to estimate (think that's more appropriate than "predict") the parameters of your distribution, is this correct? If not then feel free to revert my edit, better (clearer) titles generally attract answers! – Silverfish Apr 14 '16 at 21:55
• Yes! that makes sense. – AdrianP. Apr 14 '16 at 21:58
• Do you really mean 50th, or do you mean 50th%, i.e., 500th (or 501)? – Mark L. Stone Apr 14 '16 at 22:00

Use maximum likelihood.

Generally, suppose you have sorted array of order statistics $x_{[i]}$ for $i$ in a subset $\mathcal{I}$ of $\{1,2,\ldots, n\}$. Augment this vector to include $x_{[0]}=-\infty$ and $x_{[n+1]}=\infty$. You are supposing the underlying distribution is $F_{\theta}$ with corresponding density $f_{\theta}$ and you wish to estimate $\theta$. The likelihood of these order statistics, up to a factor that will not vary with $\theta$, is a product of the $f_\theta{x_{[i]}}$ (excluding $x_{[0]}$ and $x_{[n+1]}$, where necessarily $f=0$) times all the powers

$$\left(F_\theta(x_{[i_{j+1}]}) - F_\theta(x_{[i_j]})\right)^{i_{j+1}-i_j-1}$$

for $j$ from $1$ onward. In the example of the question this would be

$$f_\theta(0)f_\theta(8.1)f_\theta(10)\ \left(F_\theta(8.1) - F_\theta(0)\right)^{949}\left(F_\theta(10) - F_\theta(8.1)\right)^{48}.$$

The rest is completely routine, provided you have at least as many order statistics as there are components of $\theta$: there typically will be a unique value of $\theta$ that minimizes this expression. In the question, $\theta$ has two components and there are three order statistics, so all is fine. The machinery of ML produces, in a standard way, estimated standard errors of the parameters, too.

To illustrate, the R code below estimates $\hat\mu=5.39794$ and $\hat\sigma=1.62553$ from the data in the question (namely, $(x_{[1]}, x_{[951]}, x_{[1000]}) = (0, 8.1, 10)$, with $n=1000$). Then, as a quick visual check that the estimates are reasonable, it generates $1000$ datasets of size $n=1000$ from this Normal distribution, records the $1,951,$ and $1000$ order statistics (as given in the question), plots their histograms, and superimposes the observed order statistics on those histograms. The fit is beautiful for orders $1$ and $951$ and reasonable for order $1000$ (the value of $10$ is around the twelfth percentile of this distribution--not too extreme).

#
# Negative log likelihood.
#
Lambda <- function(theta) {
mu <- theta[1]
log.sigma <- theta[2]
sigma <- exp(log.sigma)
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)))
}
#
# The data.
#
n <- 1000
orders <- c(1, 951, 1000)
values <- c(0, 8.1, 10)
#
# Compute the estimate.
#
theta.start <- c(mu=mean(range(values)), log.sigma=log(diff(range(values))/6))
fit <- nlm(Lambda, theta.start) # It converged.
theta.hat <- fit$estimate mu.hat <- theta.hat[1] sigma.hat <- exp(theta.hat[2]) # # Check the quality of the estimate visually. # n.sim <- 1e3 n <- 1000 sim <- apply(matrix(rnorm(n.sim*n, mu.hat, sigma.hat), n), 2, sort)[orders, ] par(mfrow=c(1, length(orders))) invisible( sapply(1:length(orders), function(i) { hist(sim[i,], freq=FALSE, xlab="Value", main=paste("Order", orders[i])) abline(v=values[i], lwd=2, col="Red") }) )  You have a sample of random numbers$x_1,x_2,\dots,x_{1000}$. What's know is$\min[x_i],\max[x_i]$and 50th highest number A, which is 5th percentile in a sample 1000 long. Given the information I'd go with $$\hat\mu=\frac{\min+\max}{2}=5$$, and $$\hat\sigma=\frac{A-\hat\mu}{1.65}=\frac{8.1-5}{1.65}\approx 1.9$$ The t-statistics 1.65 is derived from 5% one tail z-value. I bet that you can come up with a rigorous MLE procedure which would produce the same result UPDATE Here's another way of estimating the dispersion using @jth idea. The max and min upper bounds are: $$\max=\mu+\sigma(\sqrt{2\log n} +1)$$ Hence, $$\hat\sigma=\frac{\max-\min}{2(\sqrt{2\log n} +1)}=\frac{10}{2(\sqrt{2\log 1000}+1)}\approx 1.4$$ • The max of$n$Gaussians scales like$\sqrt{2 \log n}$. So you'd want to renormalize$\hat{\mu}$. – nth Apr 27 '17 at 18:20 • @jth Gaussian's a symmetric distribution, I can't see how max could have different dynamics than min. The center must be right in between min and max – Aksakal Apr 27 '17 at 18:31 • This isn't a great solution. For the example data, we may compare this estimate to the Maximum Likelihood estimate by subtracting the log likelihoods and doubling them to produce their deviance. The deviance is around$1300$, amounting to a huge error. You have substantially overestimated the standard deviation. – whuber Apr 27 '17 at 18:42 • @jth - no, you wouldn't, because the min and max are, in expectation, symmetric around the true mean. – jbowman Apr 27 '17 at 19:11 • My mistake, I thought the order statistics were of$|X_i|$. In which case it should be$\sqrt{2\log(2n)}\$ also. – nth Apr 27 '17 at 19:13