4.5 Métodos específicos para la generación de algunas distribuciones notables
En el pasado se ha realizado un esfuerzo considerable para desarrollar métodos eficientes para la simulación de las distribuciones de probabilidad más importantes. Estos algoritmos se describen en la mayoría de los libros clásicos de simulación (e.g. Cao, 2002, Capítulo 5)18, principalmente porque resultaba necesario implementar estos métodos durante el desarrollo de software estadístico. Hoy en día estos algoritmos están disponibles en numerosas bibliotecas y no es necesario su implementación (por ejemplo, se puede recurrir a R o emplear su librería matemática disponible en https://svn.r-project.org/R/trunk/src/nmath). Sin embargo, además de que muchos de ellos servirían como ilustración de la aplicación de los métodos generales expuestos en secciones anteriores, pueden servir como punto de partida para la generación de otras distribuciones.
Entre los distintos métodos disponibles para la generación de las distribuciones continuas más conocidas podríamos destacar:
Método de Box-Müller para la generación de normales independientes (que se puede generalizar para otras distribuciones o incorporar dependencia).
Algoritmos de Jöhnk (1963) y Cheng (1978) para la generación de la distribución beta (como ejemplo de la eficiencia de los métodos de aceptación-rechazo).
4.5.1 Método de Box-Müller
Se basa en la siguiente propiedad. Dadas dos variables aleatorias independientes \(E \sim \exp\left( 1\right)\) y \(U \sim \mathcal{U}( 0, 1 )\), las variables \(\sqrt{2E} \cos 2\pi U\) y \(\sqrt{2E}\operatorname{sen} 2\pi U\) son \(\mathcal{N}( 0, 1 )\) independientes.
Algoritmo 4.8 (de Box-Müller 1958)
Generar \(U,V\sim \mathcal{U}(0, 1)\).
Hacer \(W_1=\sqrt{-2\ln U}\) y \(W_2=2\pi V\).
Devolver \(X_1=W_1\cos W_2\), \(X_2=W_1\operatorname{sen}W_2\).
Podemos hacer que la función rnorm() de R emplee este algoritmo estableciendo el parámetro normal.kind a "Box-Muller" en una llamada previa a set.seed() o RNGkind().
Este método está relacionado con el denominado método FFT (transformada de Fourier; e.g. Davies y Harte, 1987) para la generación de una normal multidimensional con dependencia, que resulta ser equivalente al Circular embedding (Dietrich and Newsam, 1997). La idea de estos métodos es que, considerando módulos exponenciales y fases uniformes generamos normales independientes, pero cambiando la varianza de los módulos (\(W_1\)) podemos inducir dependencia. Adicionalmente, cambiando la distribución de las fases (\(W_2\)) se generan distribuciones distintas de la normal.
4.5.2 Simulación de la distribución beta
Existen multitud de algoritmos para simular la distribución \(\mathcal{Beta}(a, b)\). Probablemente, el más sencillo de todos es el que se obtiene a partir de la distribución gamma o de la chi-cuadrado, si se dispone de un algoritmo para generar estas distribuciones, empleando la definición habitual de la distribución beta:
Si \(Y \sim \mathcal{Gamma}(a, s)\) y \(Z \sim \mathcal{Gamma}(b, s)\) son independientes, entonces \[X=\frac{Y}{Y+Z} \sim \mathcal{Beta}(a, b).\]
Como la distribución resultante no depende de \(s\) y \(\chi^2_{n} \overset{d}{=} \mathcal{Gamma}\left(\tfrac{n}{2}, \tfrac{1}{2}\right)\), se podría considerar \(Y \sim \chi^2_{2a}\) y \(Z \sim \chi^2_{2b}\) independientes.
También se podrían emplear resultado conocidos relacionados con esta distribución, como por ejemplo que la distribución del estadístico de orden \(k\) de una muestra de tamaño \(n\) de una distribución uniforme tiene una distribución beta: \[U_{(k)} \sim \mathcal{Beta}(k,n+1-k).\] El resultado es el algoritmo de Fox (1963), que podría ser adecuado para simular esta distribución cuando \(a, b \in \mathbb{N}\) y son valores pequeños.
Algoritmo 4.9 (Fox, 1963)
Generar \(U_1, U_2, \ldots, U_{a+b-1} \sim \mathcal{U}(0, 1)\).
Ordenar: \(U_{(1)}\leq U_{(2)}\leq\cdots\leq U_{(a+b-1)}\).
Devolver \(X=U_{(a)}\).
Es obvio que este algoritmo puede resultar muy lento si alguno de los dos parámetros es elevado (pues habrá que simular muchas uniformes para conseguir un valor simulado de la beta). Además, en función de cuál de los dos parámetros, \(a\) ó \(b\), sea mayor, resultará más eficiente, en el paso 2, comenzar a ordenar por el mayor, luego el segundo mayor, y así sucesivamente, o hacerlo empezando por el menor. En cualquier caso, es obvio que no es necesario ordenar todos los valores \(U_{i}\) generados, sino tan sólo encontrar el que ocupa el lugar \(a\)-ésimo.
Un método válido aunque \(a\) ó \(b\) no sean enteros es el dado por el algoritmo de Jöhnk (1964).
Algoritmo 4.10 (Jöhnk, 1964)
Generar \(U_1, U_2\sim \mathcal{U}(0, 1)\).
Hacer \(V = U_1^{\frac1a}\), \(W = U_2^{\frac1b}\) y \(S = V+W\).
Si \(S \leq 1\) devolver \(X = \frac VS\),
en caso contrario volver al paso 1.
El método resulta extremadamente ineficiente para \(a\) ó \(b\) mayores que 1. Esto es debido a que la condición \(S\leq1\) del paso 3 puede tardar muchísimo en verificarse. Por este motivo, el algoritmo de Jöhnk sólo es recomendable para \(a<1\) y \(b<1\). Como remedio a esto puede usarse el algoritmo de Cheng (1978) que es algo más complicado de implementar19 pero mucho más eficiente.
Algoritmo 4.11 (Cheng, 1978)
Inicialización:
Hacer \(\alpha = a + b\).
Si \(\min(a,b) \leq1\) entonces hacer \(\beta=\frac1{\min( a,b)}\), en otro caso hacer \(\beta=\sqrt{\frac{\alpha-2}{2pq-\alpha}}\).
Hacer \(\gamma=a+\frac1\beta\).
Simulación:
Generar \(U_1, U_2\sim \mathcal{U}(0, 1)\).
Hacer \(V=\beta\cdot\ln\left( \frac{U_1}{1-U_1}\right)\) y \(W=a\cdot e^{V}\).
Si \(\alpha\cdot\ln\left( \frac\alpha{b+W}\right) +\gamma V-\ln4 \ge \ln\left( U_1^{2}U_2\right)\) devolver \(X=\frac W{b+W}\),
en caso contrario volver al paso 1.