Método de Composición
Caso Especial
Método Alias (Walter 1977)
Permite generar de manera eficiente v.a.d. Con soporte
finito. Supongamos que se desea generar la v.a.d. X con
función de cuantía P = { pi : i = 1,2,...,n }
P
1
n 1
Q
(k )
n  1 k 1
donde Q(k) es una distribución concentrada en a lo sumo
dos puntos {1,2,...,n}.
La demostración de esta
descomposición se basa en:
Transformaciones
Lema1: Sea P = { pi : i=1,2,...,n} función de cuantía
Entonces:
a) Existe i {1,2,...,n} tal que pi <  n11 
b) Para tal i, existe j con i  j tal que pi + pj   n11 
Dem : Reducción al absurdo
Métodos Específicos
Distribución Binomial
Para generar una v.a.d. X ~ B(n,p)
n
X
Z
i 1
i
;
Z i ~ B (1, p )
independientes
Algoritmo
P1 : Hacer X = 0
P2 : Efectuar n réplicas
- Generar U ~ U(0,1)
Si U < p , Hacer X = X + 1
Si U  p , Hacer X = X + 0
P3 : Generar salida X
Observación: Método propuesto requiere de generar “n”
números aleatorios y n comparaciones
Métodos Específicos
Un método alternativo es el método de inversión basado
en la relación recursiva siguiente
[Fórmula recursiva]
Sea
P( X  i  1) 
(n  i ) p
(i  1)(1  p)
P( X  i )
P  P( X  i ) ; F  P( X  i )
Métodos Específicos
Algoritmo
P1 : Genera U ~ U(0,1)
P2 : Hacer i = 0 ,
P = F = (1-p)n
Hasta que U < F
Hacer P =
(n  i ) p
(i  1)(1  p )
i=i+1
P3 : Generar salida X = i
P ,F=F+P
Métodos Específicos
Distribución Poisson
Para generar la distribución de Poisson P() con 
pequeño, utilizando el método de inversión.
P(X = i + 1) =

(i  1)
P( X  i )
usando P = P(X = i) , F = P(X  i)
Métodos Específicos
Algoritmo
P1 : Genera U ~ U(0,1)
P2 : Hacer i = 0
F = P = Exp(-)
Hasta que U < F

Hacer P = (i  1) P , F = F + P
i=i+1
P3 : Generar salida X = i
Métodos Específicos
Distribución Geométrica
Para generar una v.a.d. X ~ G(p), es posible discretizar
Y ~ exp().
Sea X = [y]
Entonces P[x = r] =P(r  Y < r +1), r=0,1,2,..
r 1
=   exp( s)ds  exp( r )  exp(  (r  1));
r
es la función de cuantía de una Geo(p=1-exp(-))
lnU
[
Tomando  = -ln(1-p)  X = ln(1 p ) ] ~ Geo(p)
Métodos Específicos
Distribución Hipergeométrica
Para generar una distribución Hipergeométrica H(m,n,p)
se efectúan n extracciones sin reposición de un conjunto
de m elementos de dos clases {p m  C1 y m(1-p)  C2 }
Algoritmo
P1 : Hacer X = 0, C1 = mp C2 = m-C1
P2 : Repetir n veces
Generar U ~ U(0,1)
Si U  C1/m hacer X = X+1 , C1 = C1 - 1
sino , C2 = C2 - 1
Hacer m = m - 1
P3 : Generar salida X
Métodos Específicos
Distribuciones Multivariadas
Distribuciones Independientes
El caso más simple lo constituye el de distribuciones
marginales independientes
p
F ( x )   Fxi ( xi )
i 1
con x = (x1, x2,...,xp) Basta con generar cada componente
Xi, como distribución univarianda y salir con
X = (X1, X2, ..., Xp)
Métodos Específicos
Distribuciones Dependientes
Distribuciones Dependientes con condicionadas
disponibles. Utilizando la descomposición
F(x) = F1(x1) • F2(x2 / x1)...• F(xp / x1,x2,...,xp-1)
Si disponemos de las distribuciones
Xi / X1, ..., Xi-1
i = 1,2,...,p
Algoritmo
P1 : Desde i=1,2,...,p
Generar Xi ~ Xi / x1, ..., xi-1
P2 : Generar salida x = (x1,x2,...,xp)
Métodos Específicos
Estadísticos de Orden
Para muestrear (X(1), X(2),...,X(p)), el estadístico de orden
asociado a m.a.s. X1,X2,...,Xp de X. La forma obvia de
muestrear es hacerlo de (X1,X2,...,Xp). Alternativamente,
podemos generar la muestra de orden. Por ejemplo, si
conocemos la inversa generalizada F, podemos generar
números aleatorios (U(1), U(2),...,U(p)) y salir X(i) = F(U(i)).
Para ello es necesario generar una muestra ordenada de
números aleatorios (U(1), U(2),...,U(p)) .
Métodos Específicos
Algoritmo
P1 : Generar U(1), U(2),...,U(p) ~ U(0,1)
P2 : Hacer U(p) = (Up)1/p
U(k) = U(k+1) Uk1/k
Métodos Específicos
Distribuciones Discretas
Las distribuciones discretas multivariadas no difieren de las
univariadas. El soporte puede ser grande, pero los
métodos, inversión, alias, etc. funcionan bien.
Ejemplo : Distribución bivariada (X,Y) con soporte
{1,2,...,L}x{1,2,...,M} tenemos
Pxy = P(X  x) + P(X=x, Y=y)
indexado en x.
Métodos Específicos
Métodos Específicos
Para generar X = (X1, X2,...,Xp) ~ N(, ) se usa el método
de descomposición de Cholesky.
Sea  = L Lt, para alguna matriz L.
Entonces si Z = (Z1, Z2,...,Zp) ~ N(0, Ip)
la variable X = (, LZ) ~ N(, )
Métodos Específicos
Distribución de Wishart
Para generar una v.a.c. W ~ W(n,,) para  = 0, si  = LLt
y V =  Zi Zit ; Zi normales p-variantes N(0, Ip) , i = 1,2,...,n
n
i 1
Entonces:
W = L V Lt ~ W (n,,0)
Métodos Específicos
Algoritmo
P1 : Generar Zij ~ N(0,1)
n
P2 : Hacer V = Zi Zit
i 1
P3 : Hacer W = L V Lt
P4 : Salida W
i = 1,2,...,n j=1,2,...,n
Métodos Específicos
El algoritmo implica generar “np” normales estándar. Una
reducción del esfuerzo de cálculo se obtiene utilizando la
descomposición de Bartlett.
En el caso no centrado (  0),  es una matriz simétrica
definida no negativa. Sea  =  t su descomposición de
Cholesky y u1, u2, ..., up las filas de .
Entonces, podemos escribir :
p
W   (  k  LZ k )(  k  LZ k ) 
t
k 1
n
 LZ k Z k L
t
t
k  p 1
donde se genera W, similar al caso  = 0 usando np
normales estándares.
Métodos Específicos
Distribución Multinomial (p-dimensional).
Para generar la Distribución Multinomial de parámetros
q1, q2, ..., qp X = (X1, X2, ..., Xp) ~ M(n, q1,...,qp) con :
p
q
i
p
 1,
qi  0
i 1

X i  n,
i 1
Como Xi ~ B(n, qi)
i = 1,2,...,p
Xi / X1=x1,..., Xi-1=xi-1, ~ B(n-x1...-xi-1, wi)
qi
i = 2,3,...,p
con wi = 1  q1  q2 ......  qi 1
i  1,2,..., p
Métodos Específicos
Entonces resulta el Algoritmo
P1 :
Hacer mientras m=n
i=1, w=1, Xi = 0, i=1,...,p
Mientras m  0
Generar Xi ~ B(m, qi/w)
Hacer m = m-Xi , w =1 - qi , i = i+1
Métodos Específicos
Generación de Procesos Estocásticos
Generación de Familias de v.a. {Xt}t T
Comenzaremos con las cadenas de Markov homogéneas.
Cadena de Markov en Tiempo Discreto
Para generar una cadena de Markov con espacio de estado S
y matriz de transición P = [pij] donde pij = P(Xn+1=j / X = i). La
forma más simple de simular la transición (n+1)-ésima,
conocida Xn, es generar Xn+1~{pxnj : j  S}
Métodos Específicos
Alternativamente se puede simular Tn, el tiempo hasta el
siguiente cambio de estado y, después el nuevo estado
Xn+Tn. Si Xn = s, Tn ~ G(pss) y Xn+Tn tiene una distribución
discreta con cuantía {psj / (1 - pss) : j  S \ {s}}.
Para muestrear N transiciones de la cadena suponiendo
Xo = io
Algoritmo
Hacer t=0, Xo = io
Mientras t < N
Generar h ~ G(pxtxt)
Generar Xt+h ~ {pxtj / (1 - pxtxt) : j  S \ {s}}.
Hacer t=t+h
Métodos Específicos
OBS. 1) La primera forma de simular una cadena de
Markov, que corresponde a una estrategia
sincrónica, es decir en la que el tiempo de
simulación avanza a instantes iguales.
2) La estrategia asincrónica es más complicada de
simular [Ver. B. Ripley 1996]
Métodos Específicos
Cadenas de Markov en Tiempo Continuo
La simulación asincrónica de cadenas de Markov en
tiempo continuo es sencilla de implantar.
- Las cadenas de Markov de Tiempo Continuo vienen
caracterizadas por los parámetros vi de las distribuciones
exponenciales de tiempo de permanencia en el estado i y
la matriz de transición P; con pii = 0;  pij = 1
i j
- Sea Pi la distribución de la fila i-ésima. Entonces si Xo= io,
para simular hasta T se tiene :
Métodos Específicos
Algoritmo
Hacer t = 0, Xo = io , j = 0
Mientras t < N
Generar tj ~ exp(vxj)
Hacer t = t + tj
Hacer j = j + 1
Generar Xj ~ Pxj-1
Métodos Específicos
Proceso de Poisson
En el Proceso de Poisson P(), el número de eventos NT
en un intervalo (0,T) es P(T) y los NT ~ U(0,T)
Algoritmo
- Generar NT ~ P(T)
- Generar U1, ..., UT ~ U(0,T)
Métodos Específicos
OBS :
1) Para procesos de Poisson no homogéneos, con
t
intensidad (t) y u(t) =  (s) ds . Entonces
0
- Generar NT ~ P(u(t))
- Generar T1, T2 ,..., TNT ~

 (t )
I[ 0,T ]
2) Los procesos de Poisson son un caso particular
de los procesos de renovación. La forma de generar
los primeros se extiende a los procesos de
renovación.
Métodos Específicos
- Sean S0 = 0,
S1, S2, ...
Los tiempos de ocurrencia
- Ti = Si - Si-1 los tiempos entre sucesos.
- Para un proceso de renovación, los Ti son v.a.i.i.d. según
cierta distribución .
- Simular hasta el instante T.
Hacer S0 = 0
Mientras Si < T
Generar Ti ~ 
Hacer Si = Ti + Si-1
Hacer i = i + 1
Métodos Específicos
Procesos no Puntuales (Movimiento Browniano)
- La simulación de procesos (no puntuales) en tiempo
continuo es más complicada que la simulación de procesos
puntuales.0
- Una solución es generar procesos en suficientes
instantes discretos y aproximar la trayectoria por
interpolación.
Métodos Específicos
Como ejemplo, consideremos el movimiento Browniano
con parámetro 2
- X0 = 0
- Para s1  t1 s2  t2 .....  sn  tn las v.a.
Xt1 - Xs1, ..., Xtn - Xsn son independientes
- Para s < t, Xt - Xs ~ N(0, (t-s) 2)
- Las trayectorias son continuas
Métodos Específicos
Entonces para t fijo,
Hacer X0 = 0
Desde i = 1 hasta n
Generar Yi ~ N(0, (t-s) 2)
Hacer Xit = X(i-1)t + Yi
Interpolar la trayectoria en {(it, Xit)}
Otros ejemplos de Simulación de Procesos continuos [Ver
B. Ripley 1987]
Métodos Específicos
El Proceso de Gibbs
El creciente interés en los métodos de cadenas de Markov,
se debe al uso en Inferencia Bayesiana del Muestrador de
Gibbs. [Geman (1984)]
Ejemplo: Sean (X,Y) v.a.d. Bernoulli con distribución
x
0
1
0
1
y P(X,Y)
0
p1
0
p2
1
p3
1
p4
 pi = 1
pi > 0
Métodos Específicos
P(X=1)
= p2 + p4 (Marginal)
P(X/Y=1)
 p3 x  0
=
 p4 x  1
P(X=1/Y=1) =
p4
p3  p 4
Las Distribuciones condicionales
 P( y  0 / x  0) p( y  1 / x  0) 
Ayx  

 P( y  0 / x  1) P( y  1 / x  1) 
p3
 p p1p

p1  p3
1
3


Ayx 
p4
 p2

p2  p4 
 p2  p4
Métodos Específicos

Axy  


p1
p2
p1  p2
p1  p2
p3
p4
p3  p4
p3  p4




Algoritmo
Escoger Y0 = y0 , j =1
Repetir
Generar Xj ~ X/Y = yj-1
Generar Yj ~ Y/X = xj
j=j+1
Entonces {Xn} define una cadena de Markov con matriz de
transición
A = Ayx Axy
Métodos Específicos
Como las probabilidades pi > 0, la cadena es ergódica y
tiene distribución límite, que es la marginal de X
Xn
X ; Yn
Y ; (Xn, Yn)
(X,Y)
OBS: 1) El procedimiento descrito se llama muestrador de
Gibbs [Gibbs Sampler] y nos proporciona una
cadena de Markov, con distribución límite deseada y
se puede generalizar.
Para muestrear un vector aleatorio p-variante
X = (X1, X2, ..., Xp) con distribución , conociendo
las distribuciones condicionadas Xs/Xr, r  s
Métodos Específicos
Sea  (xs/xr, r  s) Dist. Condicionada
El [Gibbs Sampler] en este caso es
- Escoger X10, X20,..., Xp0 ; j = 1
Repetir
Generar X1j ~ X1/ X2j-1,..., Xpj-1
Generar X2j ~ X2/ X1j, X3j-1,..., Xpj-1
....
Generar Xpj ~ Xp/ X1j, X2j,..., Xp-1j
j = j+1
Métodos Específicos
Se puede verificar que Xn = (X1n, X2n,..., Xpn) define una
cadena de Markov con Matriz de transición
p
n 1
n
n 1

(
x
/
X
;
j

i
;
X
, j  i)
Pg(Xn, Xn+1) = 
i
j
j
j 1
Bajo condiciones suficientemente generales [Ver Roberts
Smith (1994)]
Métodos Específicos
Ejemplo : Muestrear la densidad
 (x1/x2) =
siendo D =
R+ 
1

exp[  x1 (1  x2 )] I ( x1 , x2 )
2
D
R
 (x1/x2) = ( x(1x, x)2 )  exp[  x1 (1  x22 )]
2
 (x2/x1) = exp[  x1 x2 ]
2
x1/x2 ~ exp[1  x22 ]
x2/x1 ~ N(0, 2=(1/2x1))
Métodos Específicos
El muestreador Gibbs
Escoger x20 ; j = 1
Repetir
Generar X1j ~ exp[1+(x2j-1)2]
Generar X2j ~ N(0, 1/2x1j)
OBS: Las secuencias podrían efectuarse en forma
aleatoria en lugar de usar la secuenciación natural
Estudiar el Algoritmo de Metropolis-Hastings.
Descargar

Generación de Variables Aleatorias