# How do you find the population size N based on the highest n values?

For example assume $N$ people performed a selection test like GMAT. Assume the distribution of the scores is a normal distribution (but parameters are not known). If you have a list of the $n$ highest scores (approved people), how do you estimate $N$? Intuitively, it is totally possible to make such an estimation, if $n$ is equal to $N$ then you would expect the distribution of the n scores to be a normal (a few very high scores, a lot of average scores and a few very low scores) if $n\ll N$ then you would expect to see an almost linear increase in the number of scores as you look at the lowest scores (just a few high scores and a lot of low scores, where low here means low in comparison to the other available scores). The motivation is that in my country (Brazil), many times there is a test to get a public job and sometimes they do not publish how many people participated but publish a list of the $n$ people that passed and their scores. I would like to find a way to determine the $N/n$ (candidates per job position).

• Are the parameters of this "normal distribution" known or not? – whuber Apr 2 '15 at 14:45
• They are not known. Ideally I would not even assume it is normal. But since we are using only the extreme values of the real distribution, I don't think we have a better option. – Mandrill Apr 2 '15 at 15:17
• I think you are likely to have some good options (I don't know which one you are comparing to when you write "better," though). At a minimum you should do sensitivity analyses to see how much your estimates might depend on making a specific assumption about the distributional shape and parameters. – whuber Apr 2 '15 at 15:23

This is a nice question. I’ll give it a try... Denote by $\Phi$ the cdf of the normal distribution and by $\phi$ its density. The joint distribution of the $n$ highest scores $y_{(N-n+1)}, \dots, y_{(N)}$ is $${N ! \over (N-n) !} \Phi(y_{(N-n+1)})^{N-n} \phi(y_{(N-n+1)}) \cdots \phi( y_{(N)} ),$$ (that is, $N-n$ values among $N$ below $y_{(N-n+1)}$, and the others where observed).

This is a bit puzzling at first sight: the usual setting would be that you know $N$ and $n$ and you want to infer the parameters of the distribution. Here you are interested only in $N$ so the MLE is obtained by maximizing $$(N - n) \log \Phi(y_{(N-n+1)}) + \log(N!) - \log\left((N-n) !\right).$$

This is kind of intuitive: the grade of the last admitted $y_{(N-n+1)}$ and the number of admitted are all that matters.

A quick numerical experiment:

> N <- 1000
> set.seed(1)
> x <- sort(rnorm(N), decreasing=TRUE)[1:10]
> x
[1] 3.810277 3.055742 2.675741 2.649167 2.497662 2.446531 2.401618 2.350554
[9] 2.349493 2.321334
> f <- function(N, n = 10, xn = x[n])
+      (N-n)*log(pnorm(xn)) + lfactorial(N) - lfactorial(N-n)
> plot( 800:1200, sapply(800:1200, f), type="l")


This looks promising. Let's have a look on the properties of this estimator, again for $n = 10$:

> MLE <- replicate( 1e4, {x <- sort(rnorm(N), decreasing=TRUE)[1:10];
+        optimize(f, c(100,20000), maximum=TRUE, xn = x[10])$maximum} ) > mean(MLE) [1] 1112.798 > sd(MLE) [1] 393.086 > hist(MLE, breaks=40)  However from you edits to your question I think that you want to estimate both$n$and the parameters of the normal distribution. This could be done by maximizing the log of above the joint density. However for your concrete application, the underlying distribution is unlikely to be normal. It is surely a mixture between the grades of diversely prepared candidates, and I am not optimistic about the possibility of obtaining good estimates. So let’s try that again, with unknown parameters for the normal distribution: g <- function(N, mu, sd, X) { n <- length(x); (N-n)*pnorm(X[n], mean = mu, sd = sd, log.p = TRUE) + sum(dnorm(X, mean=mu, sd=sd, log=TRUE)) + lfactorial(N) - lfactorial(N-n) } > optim( c(1000,0,1), function(theta) -g(theta[1], theta[2], theta[3], x) )$par
[1] 292951.707061     -3.650634      1.498264

$value [1] -13.78503  So the MLE says here that the best guess is a mean of -3.65 (instead of 0), a standard deviation of 1.49 (instead of 1, well, ok), and a size$N = 293 000$... hu... Well let’s try again with$n = 100$: 10 observations is not much! > set.seed(17) > x <- sort(rnorm(N), decreasing=TRUE)[1:100] > optim( c(1000,0,1), function(theta) -g(theta[1], theta[2], theta[3], x) )$par
[1] 1031.1320174   -0.2112833    1.1694436

$value [1] -321.6677  Not so bad...? but if I try again with set.seed(18), the estimate for$N$is 5000...! I might be missing something, but for the moment I continue to be pessimistic. Moreover, in the real world the grades are not normal. It is not rare to have a frankly bi-modal distribution, and the right tail is often quite special. The best students/candidates are far off the distribution, I have had multiple occasions to check this. So it is wrong to rely on the right tail for making these estimates: for example, if the$n = 20$best candidates are all from a (relatively) homogeneous group of$100$very smart and well prepared candidates, you will estimate only the size of this group; you won’t have any information on the (much more numerous) less prepared candidates. • This is a nice idea, but it seems you are assuming the distribution is always standard normal, which almost surely is far from the truth. Your simulation seems to work only because it simulates standard normal variates, of course. – whuber Apr 2 '15 at 14:37 • Well yes when I started this, from the formulation of the question I understood that the underlying distribution was known. Hence my last paragraph... – Elvis Apr 2 '15 at 14:39 • Also since we are estimating a distribution only using extreme data (highest values) I think it is best to just assume normal than something more sophisticated. Once I have an implementation I can test on real data where N is know. My intuition (I trust my math intuition) is that if n is not too small (n>20) then I think it is possible to have a decent estimate of N (at least a good estimate of order of magnitude 10,100,1000,10000 etc). Thank you so far. – Mandrill Apr 2 '15 at 15:00 • I simulated data from a standard normal with$N = 1000$. Assuming that the parameters of the normal are known, the maximum likelihood is obtained by maximizing the second large expression in the text. If not, you need to maximize the (log of the) first large expression.$\Phi$and$\phi\$ are the cumulative distribution function and the density of the normal distribution, respectively. – Elvis Apr 2 '15 at 15:02
• I used the (free) programming language R! – Elvis Apr 2 '15 at 15:12