Generación de aceptación y rechazo para v.a. continuas

En ocasiones, no es viable utilizar el método de la transformación inversa. En esos casos se puede aplicar el método de rechazo propuesto por Von Neumann (1951), el cual se basa en el siguiente teorema.

Teorema Sea \(X\) una v.a. distribuida con función de densidad \(f_X(x), x\in I\), la cual se puede representar como \[ f_X(x)=ch(x)g(x), \] donde \(c\geq 1\), \(0<h(x)\leq 1\), y \(g(x)\) es una función de densidad. Sea \(U\sim \mathcal{U}(0,1)\) e \(X\) una v.a. cuya función de densidad es \(g(x)\). Entonces se cumple que, \[ f_Y(x\mid U\leq h(X))=f_X(x) \]

Para una demostración se recomienda lectura de libro de Rubinstein (1981).

Para entender la idea de como funciona esté metodo, suponga que el objetivo es simular valores de \(f_X\) en el intervalo \([a,b]\) (ver la siguiente figura)

Notese que la función \(f_X\) tiene un máximo en \(c\). En esté caso al generar valores de \(X\) en el intervalo \([a,b]\) y al generar valores de \(Y\) en el intervalo \([0,c]\), solo aceptaremos aquellos valores tales que \(Y\leq f_X(X)\), de lo contrario se rechazan. Lo anterior se puede escribir en el siguiente pseudo-algoritmo.

1. Generar \(X\sim \mathcal{U}(a,b)\).
2. Generar \(Y\sim \mathcal{U}(0,c)\).
3. Aceptamos aquellos valores de \(X\), t. q. \(Y\leq f_X(X)\).
4. De lo contrario retorne a 1.

Notese que cada pareja \((X,Y)\) están en el rectángulo \([a,b]\times [0,c]\). Solamente se aceptan aquellas parejas \((X,Y)\) que se encuentren por debajo de la curva de la función \(f_X\).

Lo anterior se puede generalizar si encontramos una función \(g\) (la cual se puede simular de forma más sencilla que \(f_X\)) y está pueda envolver a \(f_X\), e.d., \(f_X(x)\leq c\cdot g(x)\) (ver la siguiente figura)

De lo descrito anteriormente \(\frac{f_X(x)}{g(x)}\leq c\), por lo cual, \(c=\max \frac{f_X(x)}{g(x)}\). Este hecho nos permite modificar el pseudo-algoritmo de la siguiente manera.

1. Generar \(X\) de una función de desidad \(g(x)\).
2. Generar \(Y\sim \mathcal{U}(0,c\cdot g(X))\).
3. Aceptamos aquellos valores de \(Y\), t. q. \(Y\leq f_X(X)\).
4. De lo contrario retorne a 1.

Los pasos 1 y 2 permiten generar valores para \((X,Y)\) los cuales se encuentran por debajo de la curva de la función \(c\cdot g\), de estos valores solo nos quedamos con aquellos valores de las parejas \((X,Y)\) los cuales cumplen la condición \(Y\leq f_X(X)\).

Ahora en el paso 2, observemos que se generan valores de \(Y\) los cuales siguen la distribución \(\mathcal{U}(0,c g(X))\), esto lo podemos reescribir como \(Y=U c g(X)\) para valores de \(U\sim\mathcal{U}(0,1)\). Por lo cual, en el paso 3 obtenemos lo siguiente

Lo anterior nos indica que si generamos valores de \(U\sim \mathcal{U}(0,1)\) y valores de \(X\), cuya función de densidad es \(g(x)\), podemos obtener valores de \(X\) con función de densidad \(f_X(x)\) siempre y cuando se cumpla que \[ U\leq h(X)=\frac{f_X(X)}{cg(X)} \] Por lo anterior un pseudo-código, se puede escribir de la siguiente manera:

1. Generar \(U\sim\mathcal{U}(0,1)\).
2. Generar \(X\sim g\).
3. Se aceptan aquellos valores de X, t.q., \(U\leq h(X)\).
4. De lo contrario retorne a 1.

Nota Para que el método funcione se debe tener en cuenta lo siguiente:

• La funciones \(f_X\) y \(g\) debe tener el mismo soporte.
• La constante \(c\), debe ser tal que \(f_X(x)/g(x)\leq c\), e.d., \(c=\max \frac{f_X(x)}{g(x)}\).
• Debe tenerse presente que los valores de \(g\) son fáciles de generar.
• Una forma de escoger a \(g\), deber ser t. q. \(f_X(x)\leq c\cdot g(x)\), eso significa que \(c\cdot g\) debe envolver a \(f_X\).

Ejemplo 1

Generar mediante el método de aceptación y rechazo, valores de \(X\sim \beta\text{eta}(3,4)\). Su función de densidad viene dada por: \[ f_X(x)=60x^2(1-x)^3,\quad 0<x<1 \] tomado como \(g\) la función de densidad de la distribución \(\mathcal{U}(0,1)\). Ahora se debe encontrar \[ c=\max\,{\frac{f_X(x)}{g(x)}}=\max\,{60x^2(1-x)^3} \] Haciendo el proceso, encontramos que para \(x=2/5\), la función anterior encuentra el máximo en \(c=1296/625\), por lo cual, \[ h(x)=\frac{f_X(x)}{cg(x)}=\frac{3125}{108}x^2(1-x)^3 \] Por lo cual es pseudo-código quedaría de la siguiente manera

1. Generar \(U_1\sim\mathcal{U}(0,1)\).
2. Generar \(U_2\sim\mathcal{U}(0,1)\).
3. Si \(U_1\leq\frac{3125}{108}U_2^2(1-U_2)^3\) entonces hacer \(X=U_2\).
4. De lo contrario retorne a 1.

El código en R sería el siguiente

h <- function(x){
  3125/108*(x^2)*((1-x)^3)
}

X<-numeric()
for(j in 1:10000){
  k<-0
  while(k==0){
    U1<-runif(1)
    U2 <-runif(1)
    if(U1<=h(U2)){
      k=1
      X[j]=U2
    }
  }
}

head(X)

## [1] 0.38415293 0.53398019 0.24319204 0.29332210 0.07316588 0.48616212

hist(X,probability = T)
curve(dbeta(x,shape1 = 3,shape2 = 4),add = T,col="red")

Realizando una prueba de bondad de ajuste tenemos que

ks.test(X,'pbeta',shape1 = 3,shape2 = 4)

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  X
## D = 0.013818, p-value = 0.0439
## alternative hypothesis: two-sided

Lo que observamos que el p-valor es mayor que el nivel de significancia 0.05 por lo cual no se rechaza que los valores se distribuyen beta con parametros \(\alpha=3\) y \(\beta=4\).

Ejemplo 2

Generar mediante el método de aceptación y rechazo, valores de \(X\sim \text{Gamma}(3/2,1)\). Su función de densidad viene dada por: \[ f_X(x)=\frac{2}{\sqrt{\pi}}x^{1/2}e^{-x},\quad x>0 \] tomado como \(g\) la función de densidad de la distribución \(\text{Exp}(2/3)\). Ahora se debe encontrar \[ c=\max\,{\frac{f_X(x)}{g(x)}}=\max\,\frac{\frac{2}{\sqrt{\pi}}x^{1/2}e^{-x}}{\frac{2}{3}e^{-\frac{2}{3}x}}=\max\,\frac{3}{\sqrt{\pi}}x^{1/2}e^{-x/3} \] Haciendo el proceso, encontramos que para \(x=3/2\), tenemos que \[ c=\frac{f(3/2)}{g(3/2)}=\frac{3\sqrt{3}}{\sqrt{2\pi}}e^{-1/2} \] Ahora encontrando \(h(x)=\frac{f_X(x)}{cg(x)}\), obtenemos \[ h(x)=\sqrt{\frac{2}{3}}x^{1/2}e^{1/2-x/3} \] Por último, para generar valores de \(g(x)\) utilizamos el método de la transformación inversa con el cual tenemos que \(-\frac{3}{2}\log U\), genera valores de Exponenciales con parámetro \(\lambda=\frac{2}{3}\).

Por lo cual es pseudo-código quedaría de la siguiente manera

1. Generar \(U_1\sim\mathcal{U}(0,1)\).
2. Generar \(U_2\sim\mathcal{U}(0,1)\) e \(X=-\frac{3}{2}\log(U_2)\).
3. Si \(U_1\leq \sqrt{\frac{2}{3}}Y^{1/2}e^{1/2-Y/3}\) entonces hacer \(X=Y\).
4. De lo contrario retorne a 1.

El código en R sería el siguiente

h <- function(x){
  sqrt(2/3)*x^(1/2)*exp(1/2-x/3)
}

X<-numeric()
for(j in 1:10000){
  k<-0
  while(k==0){
    U1 <- runif(1)
    U2 <- runif(1)
    Y <- -3/2*log(U2)
    if(U1<=h(Y)){
      k=1
      X[j]=Y
    }
  }
}

head(X)

## [1] 2.104846081 0.678311901 4.376956944 0.009746275 2.874845436 0.143062039

hist(X,probability = T)
curve(dgamma(x,shape = 3/2,rate = 1),add = T,col="red")

Realizando una prueba de bondad de ajuste tenemos que

ks.test(X,'pgamma',shape = 3/2, rate = 1)

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  X
## D = 0.010406, p-value = 0.229
## alternative hypothesis: two-sided

Lo que observamos que el p-valor es mayor que el nivel de significancia 0.05 por lo cual no se rechaza que los valores se distribuyen beta con parametros \(\alpha=3/2\) y \(\beta=1\).

Ejercicios

1. Mediante el método de aceptación y rechazo, generar observaciones para la distribución normal estándar.
2. Mediante el método de aceptación y rechazo, generar observaciones para la siguiente función de distribución \[ F(x)=x^n,\quad 0\leq x \leq 1 \]
3. Give a method for generating a random variable having density function \[ f_X(x) = \frac{e^x}{e-1}, 0 \leq x\leq 1 \]
4. Dé dos métodos para generar una variable aleatoria con función de densidad \[ f(x) = xe^{-x}, 0 \leq x \]
5. Use the rejection method to find an efficient way to generate a random variable having density function \[ f_X(x) = \frac{1}{2}(1 + x)e^{-x} ,\quad 0 < x \]