¿Por qué la media del valor más alto de 100 extracciones de una distribución normal es diferente del percentil 98 de la distribución normal?

Question

¿Por qué la media del valor más alto de 100 extracciones de una distribución normal es diferente del percentil del 98% de la distribución normal? Parece que por definición deberían ser iguales. Pero...

Código en 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)

Imagino que estoy entendiendo mal algo sobre cuál debe ser el máximo de 100 extracciones de la distribución normal. Como lo demuestra una inesperada distribución asimétrica de los valores máximos.

Answer 1

El máximo no tiene una distribución normal. Su fdc es $\Phi(x)^{100}$ donde $\Phi(x)$ es la fdc normal estándar. En general, los momentos de esta distribución son difíciles de obtener analíticamente. Hay un artículo antiguo sobre esto de Tippett ( Biometrika , 1925) .

Answer 2

Pregunté por qué había una diferencia entre la media del máximo de 100 extracciones de una distribución normal aleatoria y el percentil 98 de la distribución normal. La respuesta que recibí de Rob Hyndman fue en gran parte aceptable, pero demasiado densa desde el punto de vista técnico para aceptarla sin revisarla. Me quedé con la duda de si era posible dar una respuesta que explicara en un lenguaje sencillo e intuitivo por qué estos dos valores no son iguales.

En última instancia, mi respuesta puede ser insatisfactoriamente circular; pero conceptualmente, la razón por la que max(rnorm(100)) tiende a ser mayor que qnorm(.98) es, en resumen, porque en promedio la más alta de 100 puntuaciones distribuidas normalmente al azar superará en ocasiones su valor esperado. Sin embargo, esta distorsión no es simétrica, ya que cuando se extraen puntuaciones bajas es poco probable que acaben siendo las más altas de las 100 puntuaciones. Cada sorteo independiente es una nueva oportunidad de superar el valor esperado, o de ser ignorado porque el valor obtenido no es el máximo de los 100 valores sorteados. Para una demostración visual, compare el histograma del máximo de 20 valores con el histograma del máximo de 100 valores, la diferencia en la inclinación, especialmente en las colas, es notable.

Llegué a esta respuesta de forma indirecta mientras trabajaba en un problema/pregunta relacionado que había planteado en los comentarios. Específicamente, si descubro que los resultados de una persona en un examen están en el percentil 95, esperaría que, en promedio, si los pongo en una habitación con otros 99 examinados, su rango sería en promedio de 95. Esto resulta ser más o menos así (código R)...

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

Como extensión de esa lógica, también esperaba que si tomaba a 100 personas en una sala y seleccionaba a la persona con la 95ª puntuación más alta, luego tomaba a otras 99 personas y les hacía hacer la misma prueba, eso de media la persona seleccionada se situaría en el puesto 95 del nuevo grupo. Pero este no es el caso (código R)...

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)

Lo que diferencia el primer caso del segundo es que en el primero la puntuación del individuo lo sitúa exactamente en el percentil 95. En el segundo caso, su puntuación puede resultar algo superior o inferior al verdadero percentil 95. Como no es posible que se clasifiquen por encima de 100, los grupos que producen una puntuación de rango 95 que en realidad está en el percentil 99 o superior no pueden compensar (en términos de rango medio) los casos en los que la puntuación de rango 95 es mucho más baja que el verdadero percentil 90. Si se observan los histogramas de los dos vectores de rango proporcionados en esta respuesta, es fácil ver que hay una restricción de rango en los extremos superiores que es consecuencia de este proceso que he estado describiendo.

Answer 3

Para ampliar un poco la respuesta de Rob, supongamos que queremos conocer la función de distribución acumulativa (FDA) del valor más alto de $N$ extracciones independientes de una distribución normal estándar, $X_1, ..., X_N$ . Llama a este valor más alto $Y_1$ , la estadística de primer orden. Entonces la FCD es:

$$ \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*} $$ donde la segunda línea sigue por independencia de los empates. También podemos escribir esto como $$F_{Y_1}(x) = F_X(x)^{100},$$ donde $F$ representa la FCD y $f$ representa la PDF de la variable aleatoria dada como subíndice de esta función.

Rob utiliza la notación estándar que $\Phi(x)$ se define como $P(X < x)$ para una norma es decir, $\Phi(x)$ es la FCD normal estándar.

La función de densidad de probabilidad (FDP) del estadístico de primer orden no es más que la derivada de la FCD con respecto a $X$ : $$f_{Y_1}(x) = 100 \cdot F_X(x)^{99} f_X(x)$$ el CDF en $x$ elevado a 99 (es decir, $N-1$ ) por el PDF en $x$ por 100 (es decir, $N$ ).

Answer 4

Hay dos cuestiones: una es la asimetría en la distribución del valor máximo que has identificado; la otra es que no deberías estar mirando el percentil 98.

En lugar de la media del valor más alto, considera la mediana. Esto es más fácil, ya que se trata de una estadística de orden. La probabilidad de que los 100 valores sean menores que el cuantil $q$ es $q^{100}$ por lo que el cuantil medio para el máximo será cuando $q^{100}=\frac12$ es decir $q=\dfrac{1}{2^{1/100}}\approx 0.99309$ más que $0.98$ . Pero debido a la asimetría, se espera que la media sea aún mayor.

Como ejemplo en 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))

que da

> 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