Why is the average of the highest value from 100 draws from a normal distribution different from the 98th percentile of the normal distribution?

Question

Why is the average of the highest value from 100 draws from a normal distribution different from the 98% percentile of the normal distribution? It seems that by definition that they should be the same. But...

Code in R:

NSIM <- 10000
x <- rep(NA,NSIM)
for (i in 1:NSIM)
{
    x[i] <- max(rnorm(100))
}
qnorm(.98)
qnorm(.99)
mean(x)
median(x)
hist(x)

I imagine that I'm misunderstanding something about what the maximum of a 100 draws from the normal distribution should be. As is demonstrated by an unexpectedly asymetrical distribution of maximum values.

Answer 1

The maximum does not have a normal distribution. Its cdf is $\Phi(x)^{100}$ where $\Phi(x)$ is the standard normal cdf. In general the moments of this distribution are tricky to obtain analytically. There is an ancient paper on this by Tippett (Biometrika, 1925).

Answer 2

I asked about why there was a difference between the average of the maximum of 100 draws from a random normal distribution and the 98th percentile of the normal distribution. The answer I received from Rob Hyndman was mostly acceptable, but too technically dense to accept without revision. I was left wondering whether it was possible to provide an answer that explains in intuitively understandable plain language why these two values are not equal.

Ultimately, my answer may be unsatisfyingly circular; but conceptually, the reason max(rnorm(100)) tends to be higher than qnorm(.98) is, in short, because on average the highest of 100 random normally distributed scores will on occasion exceed its expected value. However this distortion is non-symmetrical, since when low scores are drawn they are unlikely to end up being the highest out of the 100 scores. Each independent draw is a new chance to exceed the expected value, or to be ignored because the obtained value isn't the maximum of the 100 drawn values. For a visual demonstration compare the histogram of the maximum of 20 values to the histogram of the maximum of 100 values, the difference in skew, especially in the tails, is stark.

I arrived at this answer indirectly while working through a related problem/question I had asked in the comments. Specifically, if I found that someone's test scores were ranked in the 95th percentile, I'd expect that on average if I put them in a room with 99 other test takers that their rank would average out to be 95. This turns out to be more or less the case (R code)...

for (i in 1:NSIM)
{
    rank[i] <- seq(1,100)[order(c(qnorm(.95),rnorm(99)))==1]
}
summary(rank)

As an extension of that logic, I had likewise been expecting that if I took 100 people in a room and selected the person with 95th highest score, then took another 99 people and had them take the same test, that on average the selected person would be ranked 95th in the new group. But this is not the case (R code)...

for (i in 1:NSIM)
{
    testtakers <- rnorm(100)
    testtakers <- testtakers[order(testtakers)]
    testtakers <- testtakers[order(testtakers)]
    ranked95 <- testtakers[95]
    rank[i] <- seq(1,100)[order(c(ranked95,rnorm(99)))==1]
}
summary(rank)

What makes the first case different from the second case is that in the first case the individual's score places them at exactly the 95th percentile. In the second case their score may turn out to be somewhat higher or lower than the true 95th percentile. Since they can not possibly rank higher than 100, groups that produce a rank 95 score that is actually at the 99th percentile or higher can not offset (in terms of average rank) those cases where the rank 95 score is much lower than the true 90th percentile. If you look at the histograms for the two rank vectors provided in this answer it is easy to see that there is a restriction of range in the upper ends that is a consequence of this process I have been describing.

Answer 3

There are two issues: one is the skewness in the distribution of the top value which you have identified; the other is that you should not be looking at 98th percentile.

Instead of the mean of the highest value, consider the median. This is easier as it is an order statistic. The probability that all 100 values are less than than the quantile $q$ is $q^{100}$ so the median quantile for the maximum will be when $q^{100}=\frac12$, i.e. $q=\dfrac{1}{2^{1/100}}\approx 0.99309$, rather more than $0.98$. But because of the skewness, you would expect the mean to be higher still.

As an illustration in R

require(matrixStats)
NSIM <- 100001
cases <- 100
set.seed(1)
simmat <- matrix(rnorm(cases*NSIM), ncol=cases)
tops <- rowMaxs(simmat)
c(mean(tops), median(tops), qnorm(1/2^(1/cases)))
c(pnorm(mean(tops)), pnorm(median(tops)), 1/2^(1/cases))

which gives

> c(mean(tops), median(tops), qnorm(1/2^(1/cases)))
[1] 2.508940 2.464794 2.462038
> c(pnorm(mean(tops)), pnorm(median(tops)), 1/2^(1/cases))
[1] 0.9939453 0.9931454 0.9930925

Answer 4

Just to expand on Rob's answer a bit, suppose that we want to know the cumulative distribution function (CDF) of the highest value of $N$ independent draws from a standard normal distribution, $X_1, ..., X_N$. Call this highest value $Y_1$, the first order statistic. Then the CDF is:

$$ \begin{align*}P(Y_1 < x) &= P(\max(X_1, ..., X_N) < x) \\ &= P(X_1 < x, ..., X_N < x) \\ &= P(X_1 < x) \cdot \cdot \cdot P(X_N < x) \\ &= P(X < x)^{100}, \end{align*} $$ where the second line follows by independence of the draws. We can also write this as $$F_{Y_1}(x) = F_X(x)^{100},$$ where $F$ represents the CDF and $f$ represents the PDF of the random variable given as a subscript to this function.

Rob uses the standard notation that $\Phi(x)$ is defined as $P(X < x)$ for a standard normal---i.e., $\Phi(x)$ is the standard normal CDF.

The probability density function (PDF) of the first order statistic is just the derivative of the CDF with respect to $X$: $$f_{Y_1}(x) = 100 \cdot F_X(x)^{99} f_X(x)$$ the CDF at $x$ raised to 99 (that is, $N-1$) times the PDF at $x$ times 100 (that is, $N$).