¿Qué distribución sigue el CDF normal inverso de una variable aleatoria beta?

Supongamos que define:

$$X\sim\mbox{Beta}(\alpha,\beta)$$

$$Y\sim \Phi^{-1}(X)$$

donde $\Phi^{-1}$ es el inverso del CDF de la distribución normal estándar .

Mi pregunta es: ¿hay una distribución simple que siga $Y$ , o que pueda aproximarse a $Y$ ? Pregunto porque tengo una fuerte sospecha basada en los resultados de la simulación (que se muestra a continuación) de que $Y$ converge a una distribución normal cuando $\alpha$ y $\beta$ son altos, pero no sé por qué matemáticamente. (Por supuesto, cuando $\alpha=1;\beta=1$ , $X$ sería uniforme e $Y$ sería la normal estándar, pero ¿por qué sería cierto para valores más altos?).

Si esto converge a una normal, ¿cuáles serían los parámetros de esa normal, en términos de $\alpha$ y $\beta$ ? (Espero que la media sea $\Phi^{-1}(\frac{\alpha}{\alpha+\beta})$ya que esa es la transformación del modo, pero no sé la desviación estándar).

(Dicho de otra manera, esto podría estar preguntando "¿ converge $\Phi(\mbox{Norm}(\mu, \sigma))$ a una distribución beta, para alguna dirección de $\mu$ y $\sigma$ "? No estoy seguro de si eso es más fácil de responder).

Resultados de la simulación

Aquí muestro por qué tengo la sospecha de que el resultado es normal (ya que no puedo respaldarlo con las matemáticas). La simulación de $Y$ se puede hacer en R con qnormy rnorm. Por ejemplo, eligiendo los parámetros altos $\alpha=3000$ y $\beta=7000$ :

hist(qnorm(rbeta(5000, 3000, 7000)))

Esto parece normal, y qqnormla prueba de Shapiro-Wilk (en la que la normalidad es la hipótesis nula) también sugiere lo siguiente:

qqnorm(qnorm(rbeta(5000, 3000, 7000)))

shapiro.test(qnorm(rbeta(5000, 3000, 7000)))
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  qnorm(rbeta(5000, 3000, 7000))
#> W = 0.99954, p-value = 0.2838

Para explorar la normalidad un poco más profundo, realizo 2,000 simulaciones, simulando cada vez 5,000 valores de $Y$ , luego realizo la prueba para compararla con la normal. (Elegí valores de 5K porque ese es el máximo que shapiro.testpuede manejar, y maximiza el poder para detectar desviaciones de la norma).

Si la distribución fuera realmente normal, esperaríamos que los valores p fueran uniformes (ya que el valor nulo es verdadero). De hecho, están cerca del uniforme, lo que sugiere que la distribución es muy cercana a la normal:

hist(replicate(2000, shapiro.test(qnorm(rbeta(5000, 3000, 7000)))$p.value))

Algunos experimentos demuestran que cuanto más altas son y β , más cercana es la distribución a la normalidad (por ejemplo, está bastante lejos de lo normal, pero intente y parece estar en algún punto intermedio).$\alpha$$\beta$rbeta(5000, 3, 7)hist(replicate(2000, shapiro.test(qnorm(rbeta(5000, 30, 70)))$p.value))

Sinopsis

Ha redescubierto parte de la construcción descrita en el Teorema del límite central para medianas de muestra , que ilustra un análisis de la mediana de una muestra. (El análisis obviamente se aplica, mutatis mutandis , a cualquier cuantil, no solo a la mediana). Por lo tanto, no es sorprendente que para grandes parámetros Beta (correspondientes a muestras grandes) surja una distribución Normal bajo la transformación descrita en la pregunta. Lo que es interesante es qué tan cerca de Normal es la distribución, incluso para pequeños parámetros Beta. Eso merece una explicación.

Dibujaré un análisis a continuación. Para mantener esta publicación a una longitud razonable, implica una gran cantidad de sugerentes movimientos de mano: solo pretendo señalar las ideas clave. Permítanme resumir los resultados aquí:

1. Cuando está cerca de$\alpha$ , todo es simétrico. Esto hace que la distribución transformada ya parezca Normal.$\beta$

2. Las funciones de la forma parecen bastante normales en primer lugar, incluso para valores pequeños de α y β (siempre que ambos excedan 1 y su relación no sea demasiado cerca de 0 o $\Phi^{\alpha-1}(x)\left(1-\Phi(x)\right)^{\beta-1}$$\alpha$$\beta$$1$$0$$1$ ).

3. La aparente normalidad de la distribución transformada se debe al hecho de que su densidad consiste en una densidad normal multiplicada por una función en (2).

4. A medida que aumentan y β , la desviación de la Normalidad se puede medir en los términos restantes en una serie de Taylor para la densidad logarítmica. El término de orden n disminuye en proporción a las potencias ( n - 2 ) / 2 de α y β . Esto implica que eventualmente, para α y β suficientemente grandes , todos los términos de potencia n = $\alpha$$\beta$$n$$(n-2)/2$$\alpha$$\beta$$\alpha$$\beta$$n=3$ o mayores se han vuelto relativamente pequeños, dejando solo una cuadrática: que es precisamente la densidad logarítmica de una distribución Normal.

Colectivamente, estos comportamientos explican muy bien por qué incluso para y β pequeños$\alpha$$\beta$ los cuantiles no extremos de una muestra iid Normal se ven aproximadamente normales.

---

Análisis

Debido a que puede ser útil generalizar, dejemos que sea cualquier función de distribución, aunque tenemos en mente F = $F$$F=\Phi$ .

La función de densidad de una variable Beta ( α , β ) es, por definición, proporcional a$g(y)$$(\alpha,\beta)$

$$y^{\alpha-1}(1-y)^{\beta-1}dy.$$

Si es la transformación integral de probabilidad de x y se escribe f para la derivada de F , es inmediato que x tiene una densidad proporcional a$y=F(x)$$x$$f$$F$$x$

$$G(x;\alpha,\beta)=F(x)^{\alpha-1}(1-F(x))^{\beta-1}f(x)dx.$$

Debido a que esta es una transformación monotónica de una distribución fuertemente unimodal (una Beta), a menos que sea ​​bastante extraña, la distribución transformada también será unimodal. Para estudiar qué tan cerca de lo normal podría ser, examinemos el logaritmo de su densidad,$F$

$$\log G(x;\alpha,\beta) = (\alpha-1)\log F(x) + (\beta-1)\log(1-F(x)) + \log f(x) + C\tag{1}$$

donde es una constante irrelevante de normalización.$C$

Expanda los componentes de en la serie Taylor para ordenar tres alrededor de un valor x 0 (que estará cerca de un modo). Por ejemplo, podemos escribir la expansión de log F como$\log G(x;\alpha,\beta)$$x_0$$\log F$

$$\log F(x) = c^{F}_0 + c^{F}_1 (x-x_0) + c^{F}_2(x-x_0)^2 + c^{F}_3h^3$$

por alguna con | h | ≤ | x - x 0 | . Use una notación similar para log ( 1 - F ) y log f . $h$$|h| \le |x-x_0|$$\log(1-F)$$\log f$

Términos lineales

El término lineal en se convierte así$(1)$

$$g_1(\alpha,\beta) = (\alpha-1)c^{F}_1 + (\beta-1)c^{1-F}_1 + c^{f}_1.$$

Cuando es un modo de G $x_0$ , esta expresión es cero. Tenga en cuenta que debido a que los coeficientes son funciones continuas de x 0 , ya que α y β varían, el modo x 0 también variará continuamente. Además, una vez que α y β son suficientemente grandes, eltérmino c f 1 se vuelve relativamente intrascendente. Si nuestro objetivo es estudiar el límite como α → ∞ y β → ∞ para el cual α : β permanece en proporción constante γ$G(\,;\alpha,\beta)$$x_0$$\alpha$$\beta$$x_0$$\alpha$$\beta$$c^{f}_1$$\alpha\to\infty$$\beta\to\infty$ $\alpha:\beta$$\gamma$, por lo tanto, podemos elegir de una vez por todas un punto base para el cual$x_0$

$$\gamma c^{F}_1 + c^{1-F}_1 = 0.$$

Un buen caso es donde , donde α = β en todo momento, y F es simétrica alrededor de 0 . En ese caso, es obvio x 0 = F ( 0 ) = 1 /$\gamma=1$$\alpha=\beta$$F$$0$$x_0=F(0)=1/2$ .

Hemos logrado un método por el cual (a) en el límite, el término de primer orden en la serie de Taylor desaparece y (b) en el caso especial que se acaba de describir, el término de primer orden siempre es cero.

Términos cuadráticos

Estas son la suma

$$g_2(\alpha,\beta) = (\alpha-1)c^{F}_2 + (\beta-1)c^{1-F}_2 + c^{f}_2.$$

En comparación con una distribución normal, cuyo término cuadrático es , podemos estimar que - 1 / ( 2 g 2 ( α , β ) ) es aproximadamente la varianza de G . Vamos a estandarizar G ajustando la base de x por su raíz cuadrada. Realmente no necesitamos los detalles; es suficiente entender que este cambio de escala va a multiplicar el coeficiente de ( $-(1/2)(x-x_0)^2/\sigma^2$$-1/(2g_2(\alpha,\beta))$$G$$G$$x$ en la expansión de Taylor por ( - 1 / ( 2 g 2 ( α , β ) ) ) n / 2$(x-x_0)^n$$(-1/(2g_2(\alpha,\beta)))^{n/2}.$

Término restante

Aquí está el punto clave: el término de orden en la expansión de Taylor es, según nuestra notación,$n$

$$g_n(\alpha,\beta) = (\alpha-1)c^{F}_n + (\beta-1)c^{1-F}_n + c^{f}_n.$$

Después de la estandarización, se convierte

$$g_n^\prime(\alpha,\beta) = \frac{g_n(\alpha,\beta)}{(-2g_2(\alpha,\beta))^{n/2})}.$$

Both of the $g_i$ are affine combination of $\alpha$ and $\beta$. By raising the denominator to the $n/2$ power, the net behavior is of order $-(n-2)/2$ in each of $\alpha$ and $\beta$. As these parameters grow large, then, each term in the Taylor expansion after the second decreases to zero asymptotically. In particular, the third-order remainder term becomes arbitrarily small.

The case when $F$ is normal

The vanishing of the remainder term is particularly fast when $F$ is standard Normal, because in this case $f(x)$ is purely quadratic: it contributes nothing to the remainder terms. Consequently, the deviation of $G$ from normality depends solely on the deviation between $F^{\alpha-1}(1-F)^{\beta-1}$ and normality.

This deviation is fairly small even for small $\alpha$ and $\beta$. To illustrate, consider the case $\alpha=\beta$. $G$ is symmetric, whence the order-3 term vanishes altogether. The remainder is of order $4$ in $x-x_0=x$.

Here is a plot showing how the standardized fourth order term changes with small values of $\alpha \gt 1$:

The value starts out at $0$ for $\alpha=\beta=1$, because then the distribution obviously is Normal ($\Phi^{-1}$ applied to a uniform distribution, which is what Beta$(1,1)$ is, gives a standard Normal distribution). Although it increases rapidly, it tops off at less than $0.008$--which is practically indistinguishable from zero. After that the asymptotic reciprocal decay kicks in, making the distribution ever closer to Normal as $\alpha$ increases beyond $2$.

Convergence

Suppose that $\alpha = \beta$ and let $\alpha \to \infty$ and take any small $\varepsilon > 0$. Then $var(X) \to 0$. By Chebyshev's inequality we have $\mathbb{P} [\vert X - 0.5 \vert > \varepsilon] \to 0$ and $\mathbb{P} [\vert Y \vert > \varepsilon] \to 0$. This means that $Y$ converges in probability (not in distribution actually it converges in distribution - to singleton).

Exact distribution

Denote by $f_X$ the density of beta distribution. Then your variable $Y$ has density

$$f_Y (y) = f_X ( \Phi (y) ) \phi (y).$$ Since $\Phi$ does not have a closed form I believe that this is the furthest you can get (analytically). You can try to put it into FullSimplify function in Wolfram Mathematica to see if it finds some better form.

Here is the density in R so you can plot it instead of histogram.

f_y <- function(x, alpha, beta) {
  dbeta(pnorm(x), alpha, beta) * dnorm(x)
}

Modification

However, you are maybe interested in distribution of

$$Z = \Phi^{-1} (\sqrt{\alpha} X)$$ . (still assuming $\alpha = \beta$) This may be useful because $var(\sqrt{\alpha} X) \to 1/8$ (useful because it is not zero).

Here I present a heuristic explanation (which can be made rigorous at least asymptotically). For simplicity, take $k \in \mathbb N$, $k \geq 2$. Let $X \sim \text{Beta}(k,k)$. I want to argue that $Y = \Phi^{-1}(X)$ is approximately normal.

Now let $n=2k-1$. We start by drawing $n$ i.i.d. uniformly distributed random variables $U_1, \dotsc, U_n$. Next, form the order statistics $U_{(1)} \leq \dotsc \leq U_{(n)}$.

It is well known that $U_{(k)} \sim \text{Beta}(k, n+1-k)$, thus:

$$U_{(k)} \sim \text{Beta}(k, k)$$

In other words: The sample median of $n$ i.i.d. uniformly distributed random variables is $\text{Beta}(k,k)$ distributed.

Now let's transform by $Z_i = \Phi^{-1}(U_i)$. Then by the probability integral transform, the $Z_i$ are i.i.d. normally distributed. Also form the order statistics of the $Z_i$ ($Z_{(1)} \leq \dotsc \leq Z_{(n)}$). Since $\Phi^{-1}$ is strictly increasing, it follows that:

$$\Phi^{-1}(U_{(k)}) = Z_{(k)}$$

Therefore, to show that $Y$ is approximately normal, we just have to argue that the sample median of $n$ i.i.d. normal random variables is approximately normal.

For $k$ large, this can be made precise by a central limit theorem for sample medians. For $k$ small, say $k=2$, I will let everyone's gut feeling do the speaking.

For $a \neq b$ (but not too different) one can argue similarly by using corresponding quantiles.