Cálculos Computacionales Bayesianos

Métodos de aproximación bayesianos
Técnicas y métodos
M.Sc. Christian Amao Suxo
Universidad Nacional de Ingenierı́a

Escuela Profesional de Ingenierı́a Estadı́stica
Semestre I, 2020
M.Sc. Christian Amao Suxo (UNI) Métodos de aproximación bayesianos Semestre I, 2020 1 / 54
Temario
1 Introducción
2 Simulación de variables aleatorias
Método de la transformada inversa
Métodos generales de transformación
Métodos de aceptación y rechazo
3 Método de Monte Carlo
Muestreo de Monte Carlo
Muestreo por importancia
4 Métodos de Markov Chain - Monte Carlo (MCMC)
¿Qué son los métodos MCMC?
Cadenas de Markov
Muestreo de Gibbs
Algoritmo Metropolis - Hastings
Diagnosis de convergencia
Introducción
¿Cuál es el contexto del problema?

En un problema de análisis bayesiano, toda la información requerida
es descrita en la distribución posteriori:
1 Probabilidades.
2 Valores esperados.
3 Optimización de funciones.
Sin embargo, existen casos en los cuales estas cantidades no pueden
ser obtenidas analı́ticamente. Por ejemplo, supóngase que
X|θ ∼ N (θ, σ02 ) y θ ∼ Cauchy(µ, τ ), entonces se obtiene que
2 2
π(θ|x) ∝ (τ 2 + (θ − µ)2 )−1 e−(θ−x) /(2σ ) y por tanto:
R∞ 2 2 −1 −(θ−x)2 /(2σ02 ) dθ
−∞ θ(τ + (θ − µ) ) e
E(θ|x) = R∞ (1)
2 2 −1 −(θ−x)2 /(2σ02 ) dθ
−∞ (τ + (θ − µ) ) e
Interrogante: ¿Qué métodos computacionales nos permiten

aproximar efectivamente las cantidades bayesianas de interés?
Simulación de variables aleatorias
¿Cómo se simulan variables aleatorias?
Los métodos de aproximación basados en muestreo requieren de la

generación de variables aleatorias de diversa ı́ndole.
La distribución U(0, 1) brinda la representación probabilı́stica básica
de aleatoriedad en un computador (Ver aquı́ para más detalle.).
Las simulaciones de variables aleatorias usan como base la
distribución de la variable aleatoria: U(0, 1). La generación de
cualquier distribución requerirá de la generación de variables
uniformes simuladas.
Observación: No es posible la generación de números aleatorios en el

sentido estricto, por lo que se habla de números pseudo-aleatorios.
Simulación de variables aleatorias Método de la transformada inversa
1 Introducción
Cadenas de Markov
Muestreo de Gibbs
¿En qué consiste el método de la transformada inversa?

Para una variable aleatoria continua X con densidad pX positiva en todo
su dominio, la transformada inversa usa el hecho de que
U = FX (X) ∼ U(0, 1). (2)
Ası́, si se desea simular X, entonces se simula primero U para obtener
X = FX−1 (U ). (3)
Observación: Para una variable aleatoria arbitraria X con fda FX ,

defı́nase la operación FX− (u) = inf {x; F (x) ≥ u}. Si U ∼ U(0, 1),
entonces FX− (U ) se distribuye como X.
Ejemplo
1 Si se desea simular X ∼ Exp(1), entonces se calcula
u = FX (x) = 1 − e−x y despejando x con respecto a u se obtiene
que x = −log(1 − u). Por tanto, si U ∼ U(0, 1), entonces
X = −log(1 − U ) ∼ Exp(1).
Compare los resultados obtenidos con el generador de R para variables

aleatorias exponenciales. Para ello genere 10000 valores usando rexp.
Ejercicios
1 Genere 10000 valores de la distribución logı́stica Logis(µ, β), con
µ ∈ R y β > 0, cuya densidad viene dada por
1 e−(x−µ)/β
p(x) = .
β (1 + e−(x−µ)/β )2
2 Genere 10000 valores de la distribución Cauchy C(µ, σ), con µ ∈ R y
σ > 0, cuya densidad viene dada por
1 1
p(x) = .
πσ 1 + ( x−µ
σ )
2
3 Compare los valores obtenidos con los generados por R usando la

función rlogis y rcauchy.
¿Cómo se simulan variables aleatorias discretas?
Transformación inversa para variables aleatorias discretas

Para generar valores de la variable aleatoria discreta X ∼ pX , se calculan
primero:
p0 = pX (X ≤ 0), p1 = pX (X ≤ 1), . . . , (4)
y luego se genera U ∼ U(0, 1) y se toma
X = k si pk−1 ≤ U ≤ pk . (5)
Ejemplo
Haciendo uso de la transformación inversa, genere 10000 valores de una
distribución Bin(10, 3) y compare estos valores con los generados por la
función rbinom.
Simulación de variables aleatorias Métodos generales de transformación
1 Introducción
Cadenas de Markov
Muestreo de Gibbs
¿En qué consisten los métodos generales de transformación?
Estos métodos se abren paso cuando una variable aleatoria X está

relacionada de manera simple con otra variable que es sencilla de
simular.
Ejemplos
i.i.d.
En el ejemplo anterior, se pudieron simular Xi ∼ Exp(1). Entonces, se
pueden simular las variables:
1. Y = 2 nj=1 Xj ∼ χ2 (2n), n ∈ N.
P
2. Y = β aj=1 Xj ∼ G(a, β), a ∈ N.

P
Pa
j=1 Xj
3. Y = Pa+b ∼ Be(a, b), a, b ∈ N.
j=1 Xj
Ejercicio. Simule 10000 valores de la distribución χ2 (10), G(5, 1.5) y

Be(8, 15) y compárelos con los obtenidos mediante las funciones de R
rchisq, rgamma y rbeta.
¿En qué consisten los métodos generales de transformación?
Algoritmo de Box-Muller
Si U1 y U2 son dos variables U(0, 1) independientes, entonces
p p
X1 = −2log(U1 )cos(2πU2 ) y X2 = −2log(U1 )sen(2πU2 ) (6)
son variables N (0, 1) independientes.
Ejemplo
Genere 10000 valores de una distribución normal estándar usando el
algoritmo de Box-Muller y compárelos con los generados por R con la
función rnorm.
Simulación de variables aleatorias Métodos de aceptación y rechazo
1 Introducción
Cadenas de Markov
Muestreo de Gibbs
¿En qué consisten los métodos de aceptación y rechazo?
Método de aceptación-rechazo
Supóngase que se desea simular los valores de una variable aleatoria X ∼ p
y se tiene una densidad q cuyos valores son fáciles de simular y además:
a. p y q tienen soportes compatibles (i.e. q(x) > 0 cuando p(x) > 0).
b. ∃M > 0 tal que supx {p(x)/q(x)} ≤ M .
El método consiste en:
1. Generar una variable Y proveniente de la densidad q.
2. Generar U ∼ U(0, 1) (independiente de Y ).
3. Si U ≤ Mp(Y )
q(Y ) , entonces X = Y . Caso contrario, rechazar y retornar
al paso 1.
Este método es conocido como método indirecto de simulación

porque simula valores de una densidad pX basada en los valores de
una densidad qX más fácil de simular.
Este método asegura la simulación de valores de la densidad X porque
P(Y ≤ x|U ≤ p(Y )/M q(Y )) = P(X ≤ x). (7)
Sea N : La cantidad de iteraciones hasta obtener un acierto. N es

una variable aleatoria geométrica de parámetro p = P(U ≤ Mp(Y )
q(Y ) ).
Con esto se demuestra que E(N ) = 1/p.
Se demuestra que p = 1/M y por tanto E(N ) = M . Con esto, se
debe buscar un M que sea lo más pequeño posible, i.e.
M = supx {p(x)/q(x)}
Ejercicios
1 Simule 10000 valores de una distribución Be(1.5, 4.3) usando el
método de aceptación y rechazo tomando como densidad candidata a
la distribución U(0, 1).
2 Simule 10000 valores de una distribución N (0, 1) usando el método
de aceptación y rechazo tomando como densidad candidata a la
distribución de Laplace cuya densidad está dada por
1
p(x) = e−|x| , x ∈ R
2
3 Compare los resultados obtenidos con los brindados por las funciones
de R rbeta y rnorm.
Método de Monte Carlo Muestreo de Monte Carlo
1 Introducción
Cadenas de Markov
Muestreo de Gibbs
¿En qué consiste la aproximación mediante muestreo?
Supóngase que la cantidad de interés es:

Z
E(h(X)) = h(x)pX (x)dx (8)
X
y que esta cantidad no se puede calcular analı́ticamente.

El principio del método de Monte Carlo consiste en generar una
m.a.s. (x1 , x2 , . . . , xn ) provenientes de una densidad pX y proponer
como aproximación empı́rica de (8)
n
1X
hn = h(xi ). (9)
n
i=1
De la ley fuerte de los grandes números, se tiene que

a.s.
hn → E(h(X))
Ası́ para n suficientemente grande, hn es una buena aproximación

de E(h(X)).
Si E(h2 (X)) < ∞, entonces se tiene que
Z
V ar(h(X)) 1
V (hn ) = = (h(x) − E(h(X)))2 pX (x)dx (10)
n n X
V (hn ) puede ser estimada mediante

n
1 X
V\
(hn ) = (h(xj ) − hn )2 , (11)
n(n − 1)
j=1
donde (x1 , x2 , . . . , xn ) es una m.a.s. proveniente de la densidad pX .

Bajo ciertas condiciones y por el teorema del lı́mite central:
hn − E(h(X)) D
q → N (0, 1) (12)
\
V (hn )
Con la ecuación (12) y para un n suf. grande, se pueden construir los

lı́mites para la aproximación de E(h(X)).
Ejercicios
1 Use el método de aproximación de Monte Carlo para estimar la
siguiente integral:
Z 1
(cos(50x) + sen(20x))2 dx
0
2 Replique el ejercicio anterior pero considerando como lı́mites de

integración a, b con a < b.
3 Use el método de aproximación de Monte Carlo para estimar la
integral Z ∞
xα−1 e−x dx
0
para α = 1, 2, 3, . . . , N .
Ejercicios
1 Supóngase que se desea estimar la función
R x de1 distribución de una
−t 2 /2
distribución normal estándar: Φ(x) = −∞ 2π e
√ dt
a. Exprese el problema para solucionarlo con un muestreo de Monte Carlo.
b. Calcule la varianza del estimador de Monte Carlo de Φ(x) y dé una
cota superior para la misma.
c. ¿Cuál es el número de iteraciones mı́nimo para obtener estimaciones de
Φ(x) con una precisión de al menos 4 decimales?
d. Obtenga estimaciones de Φ(x) para x = −2, 0 y 2 con la cantidad de
iteraciones hallada en c.
2 Supóngase que se desea estimar la probabilidad P(Z > 4.5) donde
Z ∼ N (0, 1). Observe lo siguiente:
a. ¿Cuántas simulaciones son necesarias para obtener un valor simulado
proveniente de una distribución Z que sea mayor a 4.5?
b. Realice una estimación por MC con n = 10000 simulaciones. ¿Las
estimaciones son buenas? ¿Por qué?
Método de Monte Carlo Muestreo por importancia
1 Introducción
Cadenas de Markov
Muestreo de Gibbs
¿En qué consiste el muestreo por importancia?
Muestreo de Monte Carlo por importancia

R
Supóngase que se desea estimar Ep (h(X)) = X h(x)p(x)dx, pero existe
una densidad alternativa q(x) (función de importancia) que tiene el
mismo soporte que h × p, entonces:
Z Z
Ep (h(X)) = h(x)p(x)dx = {h(x)w(x)}q(x)dx = Eq (h(X)w(X)) (13)
X X
donde w(x) = p(x)/q(x). La aproximación se da simulando observaciones

i.i.d X1 , . . . , Xn con densidad q con n suf. grande:
n
1X
hwn = h(xi )w(xi ) (14)
n
i=1
¿En qué consiste el muestreo por importancia?
Ejemplo
1 Mediante el método de Monte Carlo, aproxime el área bajo la curva
de la distribución normal estándar en el intervalo [−3, 3]. Para ello
use como funciones de importancia la densidad de una variable
U[−3, 3] y una distribución Cauchy truncada en [−3, 3]. ¿Qué
densidad hace que el algoritmo sea más eficiente? ¿Por qué?
2 Supóngase que se desea estimar la probabilidad P(Z > 4.5) donde
Z ∼ N (0, 1). Use un muestreo MC por importancia para estimar
dicha probabilidad usando como densidad alternativa la distribución
Exp(1) truncada en [4.5, ∞].
OJO: Un estimador de Monte Carlo es un “mejor” estimador si tiene

menor varianza (por el mismo costo computacional) que sus competidores.
Selección óptima de la función de importancia

Una función de importancia q(x) adecuada deberı́a tener las siguientes
propiedades:
1 q(x) > 0 cuando h(x)p(x) > 0.
2 q(x) deberı́a ser aproximadamente proporcional a |h(x)p(x)|, para
alcanzar la mı́nima varianza.
3 Deberı́a ser fácil simular valores provenientes de la densidad q.
4 Deberı́a ser fácil calcular la densidad q(x) para cualquier valor x que
pueda tomar.
5 Los valores de q(x) en las colas de la distribución no deben ser muy
pequeños a comparación de h(x)p(x), i.e. los valores de
h(x)p(x)/q(x) deberı́an ser “uniformes”.
Muestreo MC para estimar cantidades bayesianas

En el caso de estimaciones de cantidades bayesianes estas se expresarán:
R
g(θ)p(X|θ)p(θ)dθ
Z
E(g(θ)|X) = g(θ)p(θ|X)dθ = ΘR (15)
Θ Θ p(X|θ)p(θ)dθ
Luego, suponiendo que u es una función de importancia adecuada, la

ecuación (12) queda como:
R
g(θ)w(θ)u(θ)dθ Eu (g(θ)w(θ))
E(g(θ)|X) = ΘR = (16)
Θ w(θ)u(θ)dθ Eu (w(θ))
Donde w(θ) = p(X|θ)p(θ)/u(θ). Finalmente aproximamos la esperanza

generando θ1 , . . . , θn variables i.i.d. con densidad u con n suf. grande:
n
X n
X
π
E (g(θ)|X) ≈ w(θi )g(θi )/ w(θi ) (17)
i=1 i=1
Ejemplo
Para el estimador de bayes Normal - Cauchy:
R ∞ θ −(x−θ)2 /2
1+θ2
e dθ
δ(x) = R−∞
∞ 1 2
−(x−θ) /2 dθ
,
−∞ 1+θ2 e
resuelva los siguientes items para x = 0, 2, 4:

a. Grafique los integrandos y use el método de Monte Carlo para estimar
δ(x). Use la simulación de variables provenientes de una distribución
Cauchy estándar mediante el método de transformada inversa.
b. Repita el experimento usando el método de Monte Carlo pero esta
vez basado en la simulación de valores de una distribución normal
obtenidas mediante el algoritmo de Box-Muller.
c. Compare los resultados obtenidos en a y b. ¿Qué experimento
asegura una convergencia más rápida? ¿Por qué?
Métodos de Markov Chain - Monte Carlo (MCMC) ¿Qué son los métodos MCMC?
1 Introducción
Cadenas de Markov
Muestreo de Gibbs
Métodos de Markov Chain - Monte Carlo (MCMC) ¿Qué son los métodos MCMC?
¿En qué consisten los métodos Markov Chain - Monte Carlo?
Los métodos Markov Chain - Monte Carlo (MCMC) permiten simular

valores de una distribución posteriori de la cual no es fácil generar
muestras.
En el proceso de simulación MCMC, cada valor generado depende
solo del anterior valor simulado, de ahı́ la noción de cadena de
Markov.
Los métodos MCMC son muy usados en la computación moderna
Bayesiana ya que en la mayorı́a de los casos la forma analı́tica de
π(θ|x) es desconocida y la dimensión de θ elevada.
Métodos de Markov Chain - Monte Carlo (MCMC) Cadenas de Markov
1 Introducción
Cadenas de Markov
Muestreo de Gibbs
¿Qué es una Cadena de Markov?
Cadena de Markov
Una cadena de Markov es una secuencia de variables aleatorias discretas,
{Θi , i = 0, 1, 2, . . . } tal que la distribución de Θt dados los valores previos
Θ0 , Θ1 , . . . , Θt−1 solo depende de Θt−1 , i.e.
P (Θt ∈ A|Θ0 = θ0 , . . . , Θt−1 = θt−1 ) = P (Θt ∈ A|Θt−1 = θt−1 ), ∀A ∈ 2Ω
Nos restringiremos al caso temporalmente homogéneo, en el cual la

probabilidad P (Θt+1 = j|Θt = i) = pij no depende del tiempo t.
La matriz P = (pij )i,j∈S es llamada la matriz de transición de la
cadena {Θt }.
¿Qué propiedades tiene una Cadena de Markov?
Se dice que la distribución π es la distribución estacionaria asociada

a la cadena de Markov {Θt }, si la distribución de probabilidad de Θ0
es π y la distribución de Θt es π, ∀t ≥ 1. Ası́, diremos que {Θt }t≥0 es
una cadena de Markov con distribución estacionaria π = (πi : i ∈ S) si
π = πP , (18)
donde P es la matriz de transición de la cadena.

Una cadena de Markov {Θn }n≥0 con espacio de estados contable S y
matriz de transición P es llamada irreducible si para cualquier par de
estados i, j ∈ S, existe n ≥ 1 tal que la probabilidad de llegar al
estado j comenzando en el estado i es positiva, es decir
∃n ≥ 1 tal que P (Θn = j|Θ0 = i) > 0 (19)
Una cadena de Markov es recurrente si para cualquier estado i, si la

cadena comienza en i, esta eventualmente retornará a i con
probabilidad 1.
Una cadena de Markov es recurrente positiva si el tiempo esperado
medio de retorno al estado i es finito. Caso contrario, se dice que la
cadena es recurrente nula
Una cadena de Markov es aperiódica si para algún i ∈ S se cumple
que
m.c.d.{n : pnii > 0} = 1 (20)
Teorema
Sea {Θn }n≥0 una cadena de Markov con espacio de estados contable S y
matriz de transición P . Además, supóngase que la cadena es irreducible y
aperiódica. Si π es una distribución estacionaria de {Θn }n≥0 , entonces
para cualquier distribución inicial de Θ0
X
lim |P (Θn = j) − πj | = 0 (21)
n→∞
j
Observación: El teorema presentado postula que para n suf. grande, sea

cual fuese la distribución inicial, la distribución de probabilidad de Θn será
muy cercana a la de π.
¿Cómo ayudan las Cadenas de Markov?
1 Supongamos que queremos simular valores de una distribución

posteriori π(θ|x). La idea de los métodos MCMC consiste en simular
una cadena de Markov {Θn }n≥0 cuya distribución estacionaria sea
π(θ|x).
2 Cada valor simulado, θn , depende únicamente de su predecesor, θn−1 .
3 Es necesario simular la cadena para un número elevado de iteraciones
para aproximarse a la distribución estacionaria. Si el algoritmo se
implementa correctamente, la convergencia de la cadena está
garantizada independientemente de cuáles sean los valores iniciales.
¿Cómo ayudan las Cadenas de Markov?
1 Los primeros valores simulados, (iteraciones de burn-in), se eliminan

porque no están en el estado estacionario.
2 Una vez que se tiene que la cadena de Markov ha convergido a la
distribución estacionaria p(θ|y), entonces se entenderı́a que los valores
simulados de esta cadena pueden ser usados para hallar cantidades
bayesianas de interés.
Problema: Las realizaciones no son independientes.
Solución: El Teorema Ergódico.
Teorema Ergódico
Sean Θ1 , Θ2 , . . . , ΘN valores de una cadena de Markov {Θn } que es
aperiódica, irreducible y positiva recurrente. Si E(g(θ)) < ∞, entonces
N Z
1 X a.s.
g(Θi ) → g(θ)π(θ)dθ, (22)
N Θ
i=1
donde π es la distribución estacionaria de la cadena {Θn }.
Interrogante: ¿Qué métodos generan cadenas de Markov ergódicas que

convergan a una distribución candidata π cualquiera?
Métodos de Markov Chain - Monte Carlo (MCMC) Muestreo de Gibbs
1 Introducción
Cadenas de Markov
Muestreo de Gibbs
¿En qué consiste el muestreo de Gibbs?
Supongamos que se desea simular la densidad conjunta
p(θ) = p(θ1 , θ2 , . . . , θk )
y que para todo i = 1, 2, . . . , k; es posible simular valores de la distribución

condicional completa,
p(θi |θ −i ) = p(θi |θ1 , . . . , θi−1 , θi+1 , . . . , θk )
¿En qué consiste el muestreo de Gibbs?
Algoritmo de Gibbs
1. Fijar un valor inicial, θ 0 = (θ1,0 , θ2,0 , . . . , θk,0 )0 .
2. Para t ≥ 0, generar:
Simular θ1,t+1 ∼ p(θ1 |θ2,t , . . . , θk,t )
Simular θ2,t+1 ∼ p(θ2 |θ1,t+1 , θ3,t . . . , θk,t )
Simular θ3,t+1 ∼ p(θ3 |θ1,t+1 , θ2,t+1 , θ4,t . . . , θk,t )
..
.
Simular θk,t+1 ∼ p(θk |θ1,t+1 , . . . , θk−1,t+1 )
3. Hacer t = t + 1 y retornar al paso 2 hasta lograr convergencia.
Aplicación del algoritmo de Gibbs
Ejemplo 1: Simulación de una normal bivariada

Aplique el algoritmo de Gibbs para generar pares aleatorios (x, y)
correspondientes a una distribución normal bivariada de la forma

0 1 ρ
(X, Y ) ∼ N ,
0 ρ 1
. Para esto se hallan primero las condicionales completas:

Y | X = x ∼ N (ρx, 1 − ρ2 )
X | Y = y ∼ N (ρy, 1 − ρ2 )
Aplicación del Algoritmo de Gibbs
Ejemplo: Simulación de una normal bivariada

1. Inicializar los valores X (0) e Y (0)
2. Para j ≥ 1, simular
Y (j) ∼ N (ρX (j−1) , 1 − ρ2 )
X (j) ∼ N (ρY (j) , 1 − ρ2 )
3. Hacer j = j + 1 e iterar hasta alcanzar convergencia.
¿Cómo se aplica el muestreo de Gibbs?
Ejemplo 2
Suponga que se tiene una muestra Xi |λi ∼ P(λi ) para i = 1, 2, . . . , 10.
Además, se establecen las jerarquı́as:
iid
λi |β ∼ G(1.8, β)
β ∼ G(0.01, 1)
Si los valores observados de la muestra son
X = (5, 1, 5, 14, 3, 19, 1, 1, 4, 22),
entonces mediante el muestreo de Gibbs simule valores de la densidad

posteriori de λ = (λ1 , λ2 , . . . , λ10 )0 .
Ejemplo 2
Para aplicar el muestreo de Gibbs se requiere definir la densidad
condicional completa. Para ello se halla primero la densidad posteriori de
todos los parámetros involucrados como sigue:
p(λ, β|X) = p(X|λ, β) × p(λ|β) × p(β)

10 −λi xi
β α α−1 −βλi δ γ γ−1 −δβ

Y e λi
p(λ, β|X) = × λi e × β e
xi ! Γ(α) Γ(γ)
i=1
10
Y
p(λ, β|X) ∝ λixi +α−1 e−(1+β)λi × β 10α+γ−1 e−δβ
i=1
Ejemplo 2
Hallando las distribuciones condicionales completas se obtiene que:
p(λi |λ−i , β, X) ∝ λxi i +α−1 e−(1+β)λi

P10
p(β|λ, X) ∝ β 10α+γ−1 e−β(δ+ i=1 λi )
Con esto se obtiene que:

p(λi |λ−i , β, X) es una densidad proveniente de una distribución
G(xi + α, β + 1).
p(β|λ, X) es una
Pdensidad proveniente de una distribución
G(10α + γ, δ + 10i=1 λi ).
Métodos de Markov Chain - Monte Carlo (MCMC) Algoritmo Metropolis - Hastings
1 Introducción
Cadenas de Markov
Muestreo de Gibbs
¿En qué consiste el algoritmo de Metropolis - Hastings?

Supóngase que existe una función q(x|θ) la cual es sencilla de simular y se
desean simular los valores de una distribución π(θ):
1. Hacer i = 0 e iniciamos en θ(0) .
2. Para j ≥ 1. Sea, x = θ(j−1) .
3. Se genera y ∼ q(y | x) y u ∼ U(0, 1).
4. Si u < α(x, y) se acepta y, θ(j) = y; en otro caso, θ(j) = x, donde

π(y) q(x | y)
α(x, y) = min ,1
π(x) q(y | x)
5. Hacer j = j + 1 y retornar a 2 hasta alcanzar convergencia.
Observaciones:
En caso la distribución a querer simular se la distribución posteriori
π(θ|x), la probabilidad de aceptación α, no depende de la constante
de integración de la distribución posteriori, de modo que:
π(θ∗ )f (x|θ∗ )q(θt |θ∗ )

α = min 1,
π(θt )f (x|θt )q(θ∗ |θt )
El algoritmo MH no es muy eficiente si α es muy elevada. Esto

porque la autocorrelación entre los valores simulados serı́a alta y la
cadena convergerı́a lentamente hacia la distribución posteriori de θ.
En el caso univariado una elección habitual es usar
q(θ∗ |θt ) = N (θt , σ 2 )
donde el valor σ se ajusta para obtener un valor aceptable de α.
Métodos de Markov Chain - Monte Carlo (MCMC) Diagnosis de convergencia
1 Introducción
Cadenas de Markov
Muestreo de Gibbs
¿Cómo saber si los valores simulados son los correctos?
Cuando se ejecuta un algoritmo MCMC, es importante examinar si los

valores simulados, θt , han convergido aproximadamente a la
distribución estacionaria, π(θ|x). Para ello es recomendable:
1 Examinar empı́ricamente cómo de bien está explorando el algoritmo
MCMC el espacio de estados.
2 Verificar la convergencia
PN de las medias de los valores simulados en el
MCMC, e.g. si N1 t=1 θt → E(θ|x).
3 Analizar si los valores simulados son aproximadamente una muestra de
valores independientes e idénticamente distribuidos.
También se podrı́a ejecutar el algoritmo varias veces comenzando en
distintos valores iniciales y comprobar si el comportamiento de la
distribución estacionaria es la misma.

Para examinar si la cadena está explorando bien el espacio de estados,
se pueden graficar los valores simulados de θt frente a t y analizar si
se producen desviaciones de la estacionaridad.
(a) (b) (c)
Figura: Diferentes tipos de simulaciones para examinar si el algoritmo MCMC

explora bien el espacio de estados. (a) Buena exploración, (b) mala exploración y
(c) mala exploración.

Para examinar la convergencia de las medias posteriori de los
parámetros estimados, se pueden graficar el valor de la media muestral
de los valores simulados en función de t y verificar si han convergido.
Figura: Simulaciones de los valores estimados de E(θ|x) para examinar si dos

distintos algoritmos MCMC convergen a la media posteriori.

Finalmente, se pueden graficar las funciones de autocorrelación de los
valores generados. En general, como se simulan valores de una
cadena de Markov, los valores de θt estarán correlados positivamente.
(a) (b)
Figura: Simulación del correlograma de las cadenas de Markov mediante el
algoritmo MCMC: (a) Paseo Aleatorio MH (q(θ∗ |θ) = q(|θ∗ − θ|)), (b) Muestreo
independiente MH (q(θ∗ |θt ) = q(θ∗ )).

Cálculos Computacionales Bayesianos

Cargado por

Copyright:

Formatos disponibles

Cálculos Computacionales Bayesianos

Cargado por

Información del documento

Descripción original:

Título original

Derechos de autor

Formatos disponibles

Compartir este documento

Compartir o incrustar documentos

Opciones para compartir

¿Le pareció útil este documento?

¿Este contenido es inapropiado?

Copyright:

Formatos disponibles

Cálculos Computacionales Bayesianos

Cargado por

Copyright:

Formatos disponibles

Métodos de aproximación bayesianos

M.Sc. Christian Amao Suxo

Universidad Nacional de Ingenierı́a

¿Cuál es el contexto del problema?

Interrogante: ¿Qué métodos computacionales nos permiten

¿Cómo se simulan variables aleatorias?

Los métodos de aproximación basados en muestreo requieren de la

Observación: No es posible la generación de números aleatorios en el

¿En qué consiste el método de la transformada inversa?

Método de la transformada inversa

U = FX (X) ∼ U(0, 1). (2)

Ası́, si se desea simular X, entonces se simula primero U para obtener

Observación: Para una variable aleatoria arbitraria X con fda FX ,

¿En qué consiste el método de la transformada inversa?

Compare los resultados obtenidos con el generador de R para variables

¿En qué consiste el método de la transformada inversa?

3 Compare los valores obtenidos con los generados por R usando la

¿Cómo se simulan variables aleatorias discretas?

Transformación inversa para variables aleatorias discretas

¿En qué consisten los métodos generales de transformación?

Estos métodos se abren paso cuando una variable aleatoria X está

2. Y = β aj=1 Xj ∼ G(a, β), a ∈ N.

Ejercicio. Simule 10000 valores de la distribución χ2 (10), G(5, 1.5) y

¿En qué consisten los métodos generales de transformación?

son variables N (0, 1) independientes.

¿En qué consisten los métodos de aceptación y rechazo?

¿En qué consisten los métodos de aceptación y rechazo?

Este método es conocido como método indirecto de simulación

P(Y ≤ x|U ≤ p(Y )/M q(Y )) = P(X ≤ x). (7)

Sea N : La cantidad de iteraciones hasta obtener un acierto. N es

¿En qué consisten los métodos de aceptación y rechazo?

¿En qué consiste la aproximación mediante muestreo?

Supóngase que la cantidad de interés es:

y que esta cantidad no se puede calcular analı́ticamente.

¿En qué consiste la aproximación mediante muestreo?

De la ley fuerte de los grandes números, se tiene que

Ası́ para n suficientemente grande, hn es una buena aproximación

¿En qué consiste la aproximación mediante muestreo?

V (hn ) puede ser estimada mediante

donde (x1 , x2 , . . . , xn ) es una m.a.s. proveniente de la densidad pX .

Con la ecuación (12) y para un n suf. grande, se pueden construir los

¿En qué consiste la aproximación mediante muestreo?

2 Replique el ejercicio anterior pero considerando como lı́mites de

¿En qué consiste la aproximación mediante muestreo?

¿En qué consiste el muestreo por importancia?

Muestreo de Monte Carlo por importancia

donde w(x) = p(x)/q(x). La aproximación se da simulando observaciones

¿En qué consiste el muestreo por importancia?

OJO: Un estimador de Monte Carlo es un “mejor” estimador si tiene

Selección óptima de la función de importancia

Muestreo MC para estimar cantidades bayesianas

Luego, suponiendo que u es una función de importancia adecuada, la

Donde w(θ) = p(X|θ)p(θ)/u(θ). Finalmente aproximamos la esperanza

resuelva los siguientes items para x = 0, 2, 4:

¿En qué consisten los métodos Markov Chain - Monte Carlo?

Los métodos Markov Chain - Monte Carlo (MCMC) permiten simular

¿Qué es una Cadena de Markov?

P (Θt ∈ A|Θ0 = θ0 , . . . , Θt−1 = θt−1 ) = P (Θt ∈ A|Θt−1 = θt−1 ), ∀A ∈ 2Ω