Generación de variables aleatorias discretas

Método de la transformación inversa

Sea X una variable aleatoria discreta con función de distribución F y probabilidades puntuales

$$P(X=x_i)=p_i, \ i=1,2,3\ldots$$

Considerando la función F, que es escalonada, se tiene el siguiente algoritmo:

1.

Se hace s=p₁, i=1.

2.

Se genera $u\equiv U(0,1)$.

Si $u\leq s$, entonces x=x_i es el valor que se genera y finaliza el algoritmo. Si u>s, entonces se hace i=i+1, s=s+p_i y se repite el paso 2.

Método del alias

(Este método sólo es válido para variables cuya probabilidad está concentrada en un número finito de puntos) Sea X tal que P(X=x_i)=p_i, $i=1,2,\ldots ,k$. Tras una fase de preprocesamiento que se detalla más adelante, se tiene el siguiente algoritmo:

1.

Se genera $u\equiv U(0,1)$. Sean y=1+[ku], z=frac(ku).

2.

Si $z\leq Q(y)$, entonces k=y. Si z>Q(y), entonces k=A(y). Se toma x=x_k.

Falta determinar los valores Q(i) y los alias A(i) de modo que se tenga

$$p_i = \frac{Q(i)}{k} + \sum_{j/A(j)=i} \frac{1-Q(j)}{k}.$$

Fase de preprocesamiento:

1.

Para cada $i=1,\ldots, k$ se define

$$Q(i)=1, \, a_i=p_i, \, I_i=\text{true}.$$

2.

Se repiten las siguientes operaciones (a lo sumo) k-1 veces:

• Selecciónese i tal que $a_i<\frac{1}{k}$, I_i=true. Si esto no es posible, finaliza el preprocesamiento.
• Selecciónese j tal que $a_j>\frac{1}{k}$, I_j=true.
• Hágase I_i=false, A(i)=j, Q(i)=ka_i, $a_j=a_j-\frac{1-ka_i}{k}$.

Ejemplo: Sea $X\equiv B(3,\frac{1}{3}), \, P(X=x)= \small {\left( \begin{array}{cc} 3 \\ x ... ...array} \right)} \left( \frac{1}{3}\right)^x \left( \frac{2}{3} \right)^{3-x}$ con x=0,1,2,3. El vector de probabilidades es ( $\frac{8}{27},\frac{12}{27},\frac{6}{27},\frac{1}{27}$)

$$\begin{tabular}{\vert c\vert c\vert c\vert c\vert c\vert} \h... ... 4 & 1 & True & $\frac{1}{27}$ & - \\ \hline \end{tabular}$$

i=4, j=1, I(4)=False, A(4)=1, $Q(4)=4\frac{1}{27}$, $a(1)=\frac{8}{27}-\frac{1-Q(4)}{4}=0.08\widehat{3}$

$$\begin{tabular}{\vert c\vert c\vert c\vert c\vert c\vert} \h... ...}{27}$ & False & $\frac{1}{27}$ & 1 \\ \hline \end{tabular}$$

i=3, j=2, I(3)=False, A(3)=2, $Q(3)=4·\frac{6}{27}=\frac{24}{27}$, $a(2)=\frac{12}{27}-(\frac{1}{4}-\frac{6}{27})=0,41\widehat{6}$

$$\begin{tabular}{\vert c\vert c\vert c\vert c\vert c\vert} \h... ...}{27}$ & False & $\frac{1}{27}$ & 1 \\ \hline \end{tabular}$$

i=1, j=2, $Q(1)=4\cdot 0,08\widehat{3}=\frac{9}{27}$, A(1)=2, $a(2)=0,41\widehat{6}-(\frac{1-\frac{9}{27}}{4})=\frac{1}{4}$

$$\begin{tabular}{\vert c\vert c\vert c\vert c\vert c\vert} \h... ...27}$ & False & $\frac{1}{27}$ & 1. \\ \hline \end{tabular}$$

Ahora asignamos el 1 al valor 0, el 2 al 1, el 3 al 2 y el 4 al 3, con lo cual cuando apliquemos la segunda parte del algoritmo del alias nos dará los valores que buscamos.

Distribuciones concretas

Distribución geométrica

$$X\equiv G(p), \, P(X=x)=p(1-p)^x, \, x=0,1,2\dots$$

Representa el número de fracasos hasta que se produce el primer éxito en un experimento de Bernouilli de parámetro p.

La variable geométrica se puede relacionar fácilmente con la variable exponencial:

Sea $Y\equiv Exp(\lambda), \, F_Y(y)=1-e^{-\lambda y}, \, y\geq 0$.

Sea x>0.

$$P(x<Y\leq x+1)=F_Y(x+1)-F_Y(x)=1-e^{-\lambda(x+1)} - 1 + e^{-\lambda x} =$$

$$=-e^{-\lambda x - \lambda} + e^{-\lambda x} = \left( e^{-\lambda} \right)^x \left( 1- e^{-\lambda} \right)$$

Como $e^{-\lambda} \in (0,1)$, tomemos $\lambda$ tal que $1 - e^{-\lambda} = p$ para conseguir la expresión de probabilidad puntual de una distribución G(p). Basta tomar $\lambda=-\ln (1-p)$. Tras esto se toma un valor y generado según una $Exp(-\ln(1-p))$ y se toma x=[y]. Ya se vio que para ello hay que hacer $y=-\frac{\ln u}{-\ln (1-p)}$, con $u\equiv U(0,1)$, por lo que se concluye que

$$x=\left[ \frac{\ln u}{\ln(1-p)}\right].$$

Distribución binomial negativa

Representa el número de fracasos antes del r-ésimo éxito en un experimento de Bernouilli de parámetro p.

$$X\equiv BN(p,r), \quad P(X=x)=\left( \begin{array}{cc} x+r-1 \\ x \end{array} \right) p^r (1-p)^r, \, x=0,1,2\ldots$$

Existen varias formas de simular variables con distribución binomial negativa: la más inmediata es ir simulando experimentos de Bernouilli de parámetro p y contabilizar los fracasos hasta que se obtiene el éxito r-ésimo. Otro modo consiste en observar que una variable BN(p,r) es suma de r variables geométricas G(p) independientes, por lo que basta simular r veces una geométrica con parámetro p y sumar los valores obtenidos.

Veamos también un algoritmo directo basado en el método de composición.

Sean $X\equiv P(\lambda)$ (Poisson) y $\lambda \equiv \Gamma(\beta,\alpha)$.

$$P(X=x)=\int_0^{+\infty} {P(X=x\vert\lambda=\tilde{\lambda}) ... ...ilde{\lambda}} {\tilde{\lambda}}^{\alpha-1}} d\tilde{\lambda}=$$

$$=\frac{\beta^{\alpha}}{x!(\alpha-1)!} \int_0^{+\infty} {{\ti... ...{x! (\alpha-1)!} \frac {(x+\alpha-1)!}{(\beta+1)^{x+\alpha}} =$$

$$= \left( \begin{array}{cc} x+\alpha-1 \\ x \end{array} \right... ...{\beta+1} \right)^{\alpha} \left( \frac{1}{\beta+1} \right)^x.$$

Luego esta composición es $BN(\frac{\beta}{\beta+1}, \alpha)$. Como nosotros queremos simular una BN(p,r), entonces debemos tomar $\alpha=r$, $\beta=\frac{p}{1-p}$, resultando el siguiente algoritmo:

1.

Se genera $\lambda \equiv \Gamma(\frac{p}{1-p},r)$.

2.

Se genera $x\equiv P(\lambda)$.

Este valor x está distribuido según una distribución binomial negativa BN(p,r).

Distribución de Poisson

En la primera parte de la asignatura se vio que si los tiempos entre sucesos consecutivos son $Exp(\lambda)$, entonces el número de sucesos por unidad de tiempo sigue una distribución de Poisson de parámetro $\lambda$. En consecuencia, para generar una variable $P(\lambda)$, podemos proceder del siguiente modo: se generan valores $T_j \equiv Exp(\lambda)$ tales que $\sum_{j=1}^x T_j \leq 1 < \sum_{j=1}^{x+1} T_j$. Ese valor x corresponderá a una variable $P(\lambda)$.

Como los valores T_j son exponenciales de parámetro $\lambda$, entonces se tiene que $T_j=-\frac{1}{\lambda} \ln u_j$, con $u_j\equiv U(0,1)$, por lo que

$$\sum_{j=1}^x T_j \leq 1 < \sum_{j=1}^{x+1} T_j\Rightarrow$$

$$\sum_{j=1}^x -\frac{1}{\lambda} \ln u_j \leq 1 < \sum_{j=1}^{x+1} -\frac{1}{\lambda} \ln u_j\Rightarrow$$

$$\ln \prod_{j=1}^x u_j \geq -\lambda > \ln \prod_{j=1}^{x+1} u_j\Rightarrow$$

$$\prod_{j=1}^x u_j \geq e^{-\lambda} > \prod_{j=1}^{x+1} u_j.$$

En resumen, el algoritmo es el siguiente:

1.

Sean k=1, x=1.

2.

Se genera $u_i\equiv U(0,1)$ y se hace k=ku_i.

3.

Si $k\leq e^{-\lambda}$, entonces x-1 es el valor buscado. En caso contrario se hace x=x+1 y se va al paso 2.