La primera pregunta que surge es si el número de simulaciones necesarias para obtener un valor de X será finito. Veamos que lo es con probabilidad 1. Para ello, designemos por A la sombra de f(x) y por {Pn, n ≥ 1} una sucesión de puntos elegidos al azar en [a, b]×[0, s]. La probabilidad de que uno cualquiera de ellos, Pi, no cumpla con la condición V ≤ f(U), equivale a {Pi /∈ A} y vale P (Pi /∈ A) = s(b− a)− |A| / s(b− a) = 1− 1 / s(b− a), pues |A| = ∫ b a f(x)dx = 1. Como la elección es al azar, la probabilidad de que ninguno de ellos esté en A, P (∩n≥1{Pn /∈ A}) = lim n→∞ (1− 1 / s(b− a))n = 0, y por tanto con probabilidad 1 el número necesario de simulaciones para generar un valor de X será finito. Queda ahora por comprobar que X se distribuye con densidad f(x). En efecto, P (X ≤ x) = ∑ n≥1 P (X ≤ x,X = Un), donde {X ≤ x,X = Un} = {P1 /∈ A, . . . , Pn−1 /∈ A,Pn ∈ [a, x]× [0, s]. Al tomar probabilidades, P (X ≤ x) = ∑ n≥1 (1− 1 / s(b− a))n−1 ∫ x a f(x)dx / s(b− a) = ∫ x a f(x)dx. Ejemplo 6.2. Una variable X ∼ Be(4, 3) tiene por densidad, f(x) = { 60x3(1− x)2, 0 ≤ x ≤ 1; 0, en el resto. Su función de distribución será un polinomio de grado 6, lo que dificulta la generación de una variable aleatoria con esta distribución mediante el método de la transformación inversa. Podemos, en su lugar, recurrir al método de aceptación-rechazo. Tomemos para ello s = max0≤x≤1 f(x) = f(0.6) = 2.0736. El procedimiento consistirá en generar (U, V ) ∼ U([0, 1]× [0; 2.0736]), pero como kV ∼ U(0, k) si V ∼ U(0, 1), bastará con generar sendas U(0, 1) y hacer X = U si 2.0736V ≤ 60U3(1− U)2 =⇒ V ≤ 60U3(1− U)2 / 2.0736. Generalización del método de aceptación-rechazo El método anterior puede generalizarse utilizando una función t(x) que mayore a f(x), ∀x, en lugar de una cota superior. Si t(x) es tal que 1 ≤ ∫ ∞ −∞ t(x)dx = c < +∞, la función g(x) = t(x)/c es una densidad. El método exige que podamos generar con facilidad una variable aleatoria Y con densidad g(x) y, en ese caso, los pasos a seguir son, 1. Generamos Y cuya densidad es g(y) = t(y)/c. 2. Generamos U ∼ U(0, 1). 3. Si U ≤ f(Y )/t(Y), hacemos X = Y, en caso contrario reiniciamos el proceso desde 1. La X así generada tiene por densidad f(x) porque {X ≤ x} = {YN ≤ x} = {Y ≤ x|U ≤ f(Y )/t(Y)}, donde N designa el número de la iteración en la que se ha cumplido la condición. Al tomar probabilidades P (X ≤ x) = P (Y ≤ x|U ≤ f(Y )/t(Y)) = P (Y ≤ x,U ≤ f(Y )/t(Y)) / P (U ≤ f(Y )/t(Y)). Respecto al número N de iteraciones necesarias para obtener un valor de X, se trata de una variable geométrica con p = P (U ≤ f(Y )/t(Y)). De (6.3) deducimos que P (U ≤ f(Y )/t(Y)) = 1/c y E(N) = c, además P (N = ∞) = lim n→∞ P (N > n) = lim n→∞ ∑ j≥n+1 1/c (1− 1/c)j−1 = lim n→∞ (1− 1/c)n = 0, con lo que N será finito con probabilidad 1. Ejemplo 6.3. ¿Cómo simular una N(0, 1) con el método de aceptación-rechazo. Recordemos en primer lugar la densidad de Z ∼ N(0, 1), fZ(z) = 1√ 2π e−z2/2, −∞ < z < ∞. Si definimos X = |Z|, su densidad será fX(x) = 2√ 2π e−x2/2, 0 < x < ∞. Vamos a utilizar el método general de aceptación-rechazo para simular valores de la N(0,1) utilizando como paso previo la generación de valores de X con densidad (6.4). Para ello tomaremos como función auxiliar t(x) = √ 2e/πe−x, x ≥ 0, porque mayora a fX(x) en todo su dominio, fX(x) / t(x) = √ 2/πe−x2/2√ 2e/πe−x = exp {− (x− 1)2 / 2} ≤ 1, x ≥ 0. Por otra parte, ∫ ∞ 0 t(x)dx = √ 2e π , y por tanto g(x) = t(x)√ 2e/π = e−x, que es la densidad de una Exponencial con parámetro λ = 1. De acuerdo con el procedimiento antes descrito, 1. Generaremos sendas variables aleatorias, Y y U, la primera Exp(1) y la segunda U(0, 1). 2. Si U ≤ fX(Y) t(Y) = exp {− (Y − 1)2 / 2}, haremos X = Y, en caso contrario comenzaremos de nuevo. Una vez obtenido un valor de X el valor de Z se obtiene haciendo Z = X o Z = −X con probabilidad 1/2. El procedimiento anterior puede modificarse si observamos que aceptamos Y si U ≤ exp{−(Y − 1)2/2} ⇐⇒ − lnU ≥ (Y − 1)2/2. Pero de acuerdo con el ejemplo 6.1 − lnU ∼ Exp(1) con lo que la generación de valores de una N(0, 1) se pude llevar a cabo mediante, 1. Generaremos sendas variables aleatoria Y1, Y2, ambas Exp(1). 2. Si Y2 ≥ (Y1−1)2/2, hacemos X = Y1, en caso contrario comenzaremos de nuevo. 3. Generamos U ∼ U(0, 1) y hacemos Z = { X, si U ≤ 1/2; −X, si U > 1/2. Ejemplo 6.4. Consideremos el método polar para simular una normal estándar. Si X ∼ N(0, 1) e Y ∼ N(0, 1) y son independientes, su densidad conjunta viene dada por fXY (x, y) = 1 2π exp {−x2 + y2 / 2}. Consideremos ahora las coordenadas polares del punto (X,Y), R =√ X2 + Y 2 y Θ = arctan