Siguiente: Método de composición Nivel anterior: Métodos generales de generación de vv.aa. Previa: Método de la transformación inversa

Método de aceptación y rechazo

Sea X una variable aleatoria con función de densidad $f:I\rightarrow \mathbb{R}$ . Supóngase que f(x)=Cg(x)h(x) con C una constante, $C\geq 1$ , $0\leq g(x) \leq 1$ y h(x) una función de densidad en I. Si $U\equiv U(0,1)$ e Y (con función de densidad h) son independientes, entonces

$\begin{displaymath}h(x\vert U\leq g(Y))=f(x).\end{displaymath}$

Demostración:

$\begin{displaymath}h(x\vert U\leq g(Y)) = \frac {P(U\leq g(Y)\vert Y=x) h(x)} {\... ...\frac {P(U\leq g(x)) h(x)} {\int_I {P(U\leq g(x)) h(x)} dx} = \end{displaymath}$

$\begin{displaymath}=\frac {g(x) h(x)} {\int_I {g(x) h(x)} dx} = \frac { \frac{f(x)}{C}}{\int_I {\frac{f(x)}{C}} dx} = f(x)\end{displaymath}$

Q.E.D.

Algoritmo de aceptación y rechazo:

1.: Se genera $u\equiv U(0,1)$ y un valor y de la variable Y (de forma independiente).
2.: Si u>g(y), se va al paso 1.
3.: Si $u\leq g(y)$ , entonces se hace x=y.

Ahora interesa conocer cuál es la probabilidad de que en una iteración concreta se rechace el valor generado:

$\begin{displaymath}P(\text{aceptar}\ y_i) = P(U\leq g(Y)) = \text{denominador de la demostración} = \frac{1}{C}.\end{displaymath}$

Como el número de iteraciones del método hasta aceptar un valor sigue una distribución geométrica de parámetro $\frac{1}{C}$ (cuya esperanza es C), entonces si queremos optimizar el método habrá que intentaar que C sea próxima a uno y que h sea sencilla de generar.

$\begin{displaymath}g(x)=\frac{f(x)}{C h(x)} \leq 1 \rightarrow f(x)\leq C h(x)\end{displaymath}$

Si hacemos $\Phi(x)=C h(x)$ , entonces buscamos que sea $\Phi(x)\geq f(x)$ . Se tendrá

$\begin{displaymath}\int_I{\Phi(x)}dx = C \int_I {h(x)} dx = C \cdot 1 = C \quad \rightarrow \quad h(x)=\frac{\Phi(x)}{C}.\end{displaymath}$

Como $\Phi$ coincide con h salvo la constante y h es una función de densidad en I, entonces también se tiene que cumplir que $\int_I {\Phi(x)} dx$ sea finita.

Caso particular: f acotada, definida en un intervalo [a,b] acotado

Sea $f\leq M$ en I=[a,b].

En este caso tomamos $\Phi = M$ . Así

$\begin{displaymath}C=\int_a^b {M} dx =M(b-a),\end{displaymath}$

$\begin{displaymath}h(x)= \frac{\Phi(x)}{C} = \frac {M}{M(b-a)} = \frac {1}{b-a} \ \text{en} \ [a,b].\end{displaymath}$

$\begin{displaymath}H(x)=\frac{x-a}{b-a};\end{displaymath}$

H^-1(u)= a + (b-a)u.

Además,

$\begin{displaymath}g(x)=\frac{f(x)}{M}.\end{displaymath}$

El algortimo resultante es:

1.: Genérese u₁, $u_2 \equiv U(0,1)$ de forma independiente.
2. Si Mu₁ > f(a+(b-a)u₂), entonces se vuelve al paso 1. En caso contrario, se hace x=a+(b-a)u₂.

Ejemplo: $f(x)=\frac{2}{\pi R^2} \sqrt{R^2-x^2}, \ -R\leq x \leq R.$

En este caso [a,b]=[-R,R], por lo que $M=\frac{2}{\pi R}$ , pues f alcanza su máximo en x=0. Tratemos de simplificar el test de aceptación $M u_1\leq f(a+(b-a)u_2)$ .

$\begin{displaymath}\frac{2}{\pi R} u_1 \leq \frac{2}{\pi R^2} \sqrt { R^2-(-R+2Ru_2)^2} \end{displaymath}$

si y sólo si

$\begin{displaymath}u_1^2 \leq \frac{1}{R^2} R^2(1-(-1+2u_2)^2)\end{displaymath}$

si y sólo si

$\begin{displaymath}u_1^2 \leq (1+1-2u_2)(1-1+2u_2)=(2-2u_2)2u_2=4u_2(1-u_2)\end{displaymath}$

si y sólo si

$\begin{displaymath}\frac{u_2(1-u_2)}{u_1^2} \geq \frac{1}{4}.\end{displaymath}$

La probabilidad de aceptación es

$\begin{displaymath}\frac{1}{C}=\frac{1}{M(b-a)}=\frac {\pi R}{2} \frac{1}{2R} = \frac {\pi}{4}.\end{displaymath}$

Ejemplo: $X\equiv Be(\alpha,\beta) \ (\alpha,\beta \geq 1).\\ f(x)=\frac{1}{Be(\alpha,\beta)} x^{\alpha-1} (1-x)^{\beta-1}, \ 0<x<1$ . La función $h(x)=\alpha x^{\alpha -1}$ es función de densidad en (0,1). Así

$\begin{displaymath}H(x)=x^{\alpha} \Rightarrow \end{displaymath}$

$\begin{displaymath}H^{-1}(u)=u^{\frac{1}{\alpha}}, \ 0<u<1.\end{displaymath}$

$\begin{displaymath}f(x)=\alpha x^{\alpha -1} (1-x)^{\beta -1} \frac{1}{\alpha} ... ...)(1-x)^{\beta -1} \frac{1}{\alpha} \frac{1}{Be(\alpha,\beta)}.\end{displaymath}$

Tomamos $g(x)=(1-x)^{\beta-1}$ . De este modo se tiene que

$\begin{displaymath}C=\frac {1}{\alpha Be(\alpha,\beta)} \geq 1\end{displaymath}$

y el criterio de aceptación es

$\begin{displaymath}u_1 \leq (1-u_2^{\frac{1}{\alpha}})^{\beta-1},\end{displaymath}$

tomándose

$\begin{displaymath}x=u_2^{\frac{1}{\alpha}}\end{displaymath}$

Ejemplo:

$X\equiv \Gamma(1,\alpha), \ 0<\alpha<1.$

$\begin{displaymath}f(x)=\frac{x^{\alpha-1}e^{-x}}{\Gamma(\alpha)}, \ x\geq 0.\end{displaymath}$

Para que se cumplan todas las condiciones que necesitamos, tomamos

$\begin{displaymath}g(x)=\left \{ \begin{array}{ccc} e^{-x} & si & 0<x\leq 1, \\ x^{\alpha -1} & si & x> 1. \end{array} \right.\end{displaymath}$

$\begin{displaymath}\Phi(x)= \left \{ \begin{array}{ccc} \frac{x^{\alpha-1}}{\G... ...frac{e^{-x}}{\Gamma(\alpha)} & si & x>1. \end{array} \right. \end{displaymath}$

Es inmediato que $0\leq g \leq 1$ . Veamos que $\Phi$ integra un valor finito:

$\begin{displaymath}C=\int_0^{+\infty} {\Phi(x)} dx = \frac {1}{\Gamma(\alpha)} ... ...ha} \right]_0^1 + \left [ -e^{-x} \right]_1^{+\infty} \right)=\end{displaymath}$

$\begin{displaymath}= \frac{1}{\Gamma(\alpha)} \left( \frac{1}{\alpha} + \frac{1}{e} \right) = \frac{\alpha+e}{\alpha e \Gamma(\alpha)}.\end{displaymath}$

Ahora

$\begin{displaymath}h(x)=\left\{\begin{array}{ccc} \frac{\alpha e}{\alpha+e} x... ...{\alpha e}{\alpha + e} e^{-x} & si & x>1. \end{array} \right.\end{displaymath}$

$\begin{displaymath}H(x)=\left\{ \begin{array}{ccc} \frac{e}{e+\alpha} x^{\alp... ... \frac{\alpha}{\alpha+e} e^{1-x} & si & x>1 \end{array} \right.\end{displaymath}$

Como $H(1)=\frac{e}{e+ \alpha}$ , entonces

$\begin{displaymath}H^{-1}(u)=\left\{ \begin{array}{ccc} \left( \frac{e+\alpha... ...ght) & si & \frac{e}{e+\alpha} \leq u <1. \end{array} \right.\end{displaymath}$

El algoritmo para simular la variable X es:

1.: Se genera u₁, u₂ según una distribución U(0,1).
2.: Si $u_2\leq g(H^{-1}(u_1))$ , entonces se toma x=H^-1(u₁). En caso contrario, se repite el paso 1.