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

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.

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).

• Is there a way to express in plain language what the difference is between a given percentile and the maximum of N values? From a lay perspective, it is hard to see why a datapoint that comes from a given (Y) percentile wouldn't be expected to be (on average) the same as the top scorer from a group of 100/Y. For example, if I found that your answers were ranked in the 90th percentile, I'd expect that your answer would usually be the top answer among any randomly selected group of 10 answers. Jul 20, 2010 at 8:45
• @drknexus Your intuition is correct. However, extremes (and near-extremes) of sampling distributions are somewhat special: their values are constrained by the mass of data on one side of them, whereas--for parent distributions with infinite tails--there is no constraint at all on their values on the other side. Thus, for instance, the distribution of a maximum (from a distribution with no upper bound) is positively skewed. This increases its expectation relative to the corresponding percentile.
– whuber
Oct 29, 2010 at 13:32

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.

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


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$).

• There seems to be something wrong with your final equation (I fixed the typesetting but translated it exactly as you had it before). What is $X_1$? Also, the pdf of $X$ is not equivalent to $P(X=x)$. In fact, if $X$ has a normal distribution (or any continuous distribution) then $P(X=x) = 0$ any $x$, so that can't possibly be the pdf. Mar 31, 2012 at 7:33
• @Macro, $X_1$ is the first draw from $N$ independent draws; $Y_1$ is the first order statistic (you might prefer to write $X_{(1)}$ instead). I made the notation more precise in response to your other comments. Apr 1, 2012 at 4:20