Introducci´on a la estad´ıstica bayesiana, aplicaciones y m´etodos Parte 2 Ana Paula Palacios y Peter Diko Universidad Carlos III de Madrid

22 de Marzo de 2011 Instituto de Econom´ıa y Finanzas Facultad de Ciencias Econ´ omicas U.N.C.

(Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

1 / 42

Programa

1

Muestreo Monte Carlo

2

Modelo Lineal

3

WinBUGS

4

Modelos Jer´arquicos

(Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

2 / 42

Monte Carlo

Programa

1

Muestreo Monte Carlo

2

Modelo Lineal

3

WinBUGS

4

Modelos Jer´arquicos

(Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

3 / 42

Monte Carlo

Problema

El an´alisis bayesiano proporciona la distribuci´ on para θ, el par´ametro de inter´es f (θ|x) ∝ f (x|θ)f (θ) Tenemos inter´es en cuant´ıas relacionadas con la distribuci´on a posteriori: media,moda, mediana, cuantiles en general, intervalos de credibilidad.

(Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

4 / 42

Monte Carlo

Problema

El an´alisis bayesiano proporciona la distribuci´ on para θ, el par´ametro de inter´es f (θ|x) ∝ f (x|θ)f (θ) Tenemos inter´es en cuant´ıas relacionadas con la distribuci´on a posteriori: media,moda, mediana, cuantiles en general, intervalos de credibilidad. Complicaciones identificar la constante de normalizaci´ on de la a posteriori la distribuci´on a posteriori puede no ser tratable anal´ıticamente

(Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

4 / 42

Monte Carlo

Soluci´on Alternativa al tratamiento anal´ıtico construir una muestra θ1 , θ2 , . . . , θn de la distribuci´on a posteriori.

(Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

5 / 42

Monte Carlo

Soluci´on Alternativa al tratamiento anal´ıtico construir una muestra θ1 , θ2 , . . . , θn de la distribuci´on a posteriori. calcular la cuant´ıa de inter´es muestral

(Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

5 / 42

Monte Carlo

Soluci´on Alternativa al tratamiento anal´ıtico construir una muestra θ1 , θ2 , . . . , θn de la distribuci´on a posteriori. calcular la cuant´ıa de inter´es muestral por la Ley de los Grandes N´ umeros la distribuci´ on emp´ırica converge a la verdadera n

Z

1X θf (θ|x)dθ ≈ θi n Θ i=1

Z Θ

(Univ. Carlos III de Madrid)

(θ − µ)2 f (θ|x)dθ ≈

n 1X

n

Estad´ıstica bayesiana

(θi − µ)2

i=1

22-03-11

5 / 42

Monte Carlo

Media a posteriori para beta(1498, 1519)

(Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

6 / 42

Monte Carlo

Muestreo por inversi´on de F Enfoque para distribuciones univariantes. Necesitamos conocer la funci´ on de distribuci´ on F (x).

(Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

7 / 42

Monte Carlo

Muestreo por inversi´on de F Enfoque para distribuciones univariantes. Necesitamos conocer la funci´ on de distribuci´ on F (x). Algoritmo 1

Generamos un valor u de la distribuci´ on uniforme U(0, 1).

2

La cuant´ıa z = F −1 (u) es una observaci´ on aleatoria de F (x).

(Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

7 / 42

Monte Carlo

Muestreo por inversi´on de F Enfoque para distribuciones univariantes. Necesitamos conocer la funci´ on de distribuci´ on F (x). Algoritmo 1

Generamos un valor u de la distribuci´ on uniforme U(0, 1).

2

La cuant´ıa z = F −1 (u) es una observaci´ on aleatoria de F (x).

Comprobaremos que Z = F −1 (U), donde U es uniforme (0, 1) tiene distribuci´on F (x) P{Z ≤ x} = P{F −1 (U) ≤ x} = P{U ≤ F (x)} = F (x)

(Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

7 / 42

Monte Carlo

Rejection sampling La clave es encontrar una distribuci´ on g (x) f´acil de muestrear que cumpla para un m fijo f (x) ≤ m · g (x) en todo el soporte de f (x).

(Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

8 / 42

Monte Carlo

Rejection sampling La clave es encontrar una distribuci´ on g (x) f´acil de muestrear que cumpla para un m fijo f (x) ≤ m · g (x) en todo el soporte de f (x). Algoritmo 1

Generamos un valor z de una distribuci´ on g (x).

2

Calculamos el ratio R =

3

Generamos un valor u de una uniforme (0, 1). Acceptamos z como observaci´on aleatoria de f (x) si u < R.

(Univ. Carlos III de Madrid)

f (z) m·g (z) .

Estad´ıstica bayesiana

22-03-11

8 / 42

Monte Carlo

Rejection sampling

Ventajas No necesitamos conocer la constante de normalizazi´on. ∝ f (x) V´alido para el caso multidimensional. F´acil de implementar.

(Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

9 / 42

Monte Carlo

Rejection sampling

Ventajas No necesitamos conocer la constante de normalizazi´on. ∝ f (x) V´alido para el caso multidimensional. F´acil de implementar. Desventajas Encontrar la densidad g (x) puede ser dif´ıcil. Si la g (x) no es buena, el algoritmo puede ser ineficiente. Alta proporci´on de rechazos. Problemas en alta dimensi´ on.

(Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

9 / 42

Monte Carlo

Algoritmos MCMC

Soluci´on para casos de densidades complejas y de alta dimensi´on. Particionamos la densidad a muestrear en densidades multivariantes o univariantes m´as manejables. muestreo de una o varias dimensiones de la a posteriori exploraci´on de todo el soporte de la distribuci´ on paso por paso

(Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

10 / 42

Monte Carlo

Muestreo de Gibbs En estad´ıstica Bayesiana desde Gelfand and Smith (1990). Conocido en f´ısica antes de 1990. Algoritmo apropiado en casos cuando el muestreo de la distribuci´ on conjunta no es posible conocemos las distribuciones condicionadas para cada dimensi´on (o bloques de dimensiones)

f (θ1 , θ2 ) f (θ1 |θ2 ), f (θ2 |θ1 )

(Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

11 / 42

Monte Carlo

Muestreo de Gibbs Algoritmo 1 2

3

4

Empezamos con unos valores iniciales de θ10 , θ20 . j = 1 Generamos una observaci´ on θ1j de la distribuci´ on condicionada j−1 f (θ1 |θ2 ). Generamos una observaci´ on θ2j de la distribuci´ on condicionada j f (θ2 |θ1 ). Siguiente paso j := j + 1 y volvemos al paso 2.

(Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

12 / 42

Monte Carlo

Muestreo de Gibbs Algoritmo 1 2

3

4

Empezamos con unos valores iniciales de θ10 , θ20 . j = 1 Generamos una observaci´ on θ1j de la distribuci´ on condicionada j−1 f (θ1 |θ2 ). Generamos una observaci´ on θ2j de la distribuci´ on condicionada j f (θ2 |θ1 ). Siguiente paso j := j + 1 y volvemos al paso 2.

Obtenemos una cadena de Markov (θ10 , θ20 ), (θ11 , θ21 ), . . . , (θ1n , θ2n ) con distribuci´on estacionaria f (θ1 , θ2 ). Decartando las primeras observaciones generadas, nos quedamos con la muestra aleatoria de la distribuci´ on a posteriori conjunta. (Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

12 / 42

Monte Carlo

Muestreo de Gibbs

(Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

13 / 42

Monte Carlo

Muestreo de Gibbs

(Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

14 / 42

Monte Carlo

Algoritmo de Metropolis-Hastings Algoritmo basado en Metropolis et al. (1958) para explicar movimiento de part´ıculas. Generalizado por Hastings (1970) proporciona muestra del par´ametro θ conjunto no necesitamos conocer la constante de normalizaci´on contiene paso de aceptaci´ on-rechazo

(Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

15 / 42

Monte Carlo

Algoritmo de Metropolis-Hastings Algoritmo basado en Metropolis et al. (1958) para explicar movimiento de part´ıculas. Generalizado por Hastings (1970) proporciona muestra del par´ametro θ conjunto no necesitamos conocer la constante de normalizaci´on contiene paso de aceptaci´ on-rechazo Los candidatos se generan a partir de una distribuci´ on conveniente j−1 g (θ|θ ). Aceptaci´on del candidato se eval´ ua a base del ratio R=

(Univ. Carlos III de Madrid)

f (θC )g (θj−1 |θC ) f (θj−1 )g (θC |θj−1 )

Estad´ıstica bayesiana

22-03-11

15 / 42

Monte Carlo

Algoritmo de Metropolis-Hastings

Algoritmo 1

Empezamos con un valor inicial θ0 , j = 1.

2

Generamos un candidato θC de la distribuci´ on g (θ|θj−1 ).

3

Calculamos el ratio R=

f (θC )g (θj−1 |θC ) f (θj−1 )g (θC |θj−1 )

4

Generamos u de una distribuci´ on uniforme (0, 1). Si u < R j aceptamos el candidato θ := θC , en caso contrario θj := θj−1

5

Siguiente paso j := j + 1 y volvemos al paso 2.

(Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

16 / 42

Monte Carlo

Algoritmo de Metropolis-Hastings Obtenemos una cadena de Markov θ0 , θ1 , . . . , θn con distribuci´on estacionaria f (θ) pero Mala elecci´on del punto inicial puede complicar las cosas. Ratio de rechazos alto causar´a observaciones repetidas, mucha correlaci´on de la cadena y convergencia lenta. Ratio de rechazos bajo puede significar exploraci´ on lenta del espacio param´etrico.

(Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

17 / 42

Monte Carlo

Algoritmo de Metropolis-Hastings Obtenemos una cadena de Markov θ0 , θ1 , . . . , θn con distribuci´on estacionaria f (θ) pero Mala elecci´on del punto inicial puede complicar las cosas. Ratio de rechazos alto causar´a observaciones repetidas, mucha correlaci´on de la cadena y convergencia lenta. Ratio de rechazos bajo puede significar exploraci´ on lenta del espacio param´etrico. Un ejemplo de la distribuci´ on de propuesta g (θ|θj−1 ) θC ∼ N(θj−1 , C ) donde C puede adaptarse seg´ un el ratio de aceptaci´ on.

(Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

17 / 42

Monte Carlo

Otros enfoques Metropolis-within-Gibbs En caso de distribuci´on conjunta muy compleja f (θ) la distribuci´on se particiona f (θ1 |θ2 ), f (θ2 |θ1 ) para aplicar el algoritmo de Gibbs. Cada paso del algoritmo de Gibbs requiere generar observaciones de las condicionadas f (θ1 |θ2 ), f (θ2 |θ1 ) para lo que se emplea el algoritmo MH.

(Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

18 / 42

Monte Carlo

Otros enfoques Metropolis-within-Gibbs En caso de distribuci´on conjunta muy compleja f (θ) la distribuci´on se particiona f (θ1 |θ2 ), f (θ2 |θ1 ) para aplicar el algoritmo de Gibbs. Cada paso del algoritmo de Gibbs requiere generar observaciones de las condicionadas f (θ1 |θ2 ), f (θ2 |θ1 ) para lo que se emplea el algoritmo MH. Slice sampling f (θ) ∝ h(θ) U|θ ∼ uniforme(0, h(θ)) Se aplica muestreo de Gibbs a las condicionadas U|θ, θ|U para obtener la muestra de la distribuci´on conjunta f (θ, U) y de ah´ı f (θ). (Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

18 / 42

Monte Carlo

Diagn´ostico de convergencia Diagn´ostico de Gelman y Rubin repetimos el algoritmo MCMC m veces con puntos iniciales dispersos obtenemos 2N observaciones de cada cadena bas´andonos en las u ´ltimas N observaciones calculamos B N

varianza entre las m medias W la media de las varianzas dentro de las m cadenas aproximamos la densidad a posteriori con la distribuci´on t y denominamos df sus grados de libertad

(Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

19 / 42

Monte Carlo

Diagn´ostico de convergencia Diagn´ostico de Gelman y Rubin repetimos el algoritmo MCMC m veces con puntos iniciales dispersos obtenemos 2N observaciones de cada cadena bas´andonos en las u ´ltimas N observaciones calculamos B N

varianza entre las m medias W la media de las varianzas dentro de las m cadenas aproximamos la densidad a posteriori con la distribuci´on t y denominamos df sus grados de libertad

El factor de reducci´on p ˆ= R

s

N −1 m+1 B + N mN W



df df − 2

determina la posibilidad de reducir la variabilidad de la distribuci´on a posteriori al aumentar la muestra N → ∞. (Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

19 / 42

Modelo Lineal

Programa

1

Muestreo Monte Carlo

2

Modelo Lineal

3

WinBUGS

4

Modelos Jer´arquicos

(Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

20 / 42

Modelo Lineal

Modelo de regresi´on lineal Especificaci´on matricial Y = X β + e, e ∼ N(0, σe2 In ),

(Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

21 / 42

Modelo Lineal

Modelo de regresi´on lineal Especificaci´on matricial Y = X β + e, e ∼ N(0, σe2 In ), funci´on de verosimilitud L(β, σe2 ; X , Y )

=

(Univ. Carlos III de Madrid)

(2πσe2 )−n/2 exp



 1 T − 2 (Y − X β) (Y − X β) . 2σe

Estad´ıstica bayesiana

22-03-11

21 / 42

Modelo Lineal

Modelo de regresi´on lineal Especificaci´on matricial Y = X β + e, e ∼ N(0, σe2 In ), funci´on de verosimilitud L(β, σe2 ; X , Y )

=

(2πσe2 )−n/2 exp



 1 T − 2 (Y − X β) (Y − X β) . 2σe

Estimadores de m´axima verosimilitud: ˆ =σ ACOV (β) ˆe2 (X T X )−1 ,  2 1/2 2ˆ σe 2 SE (ˆ σe ) = n

βˆ = (X T X )−1 (X T Y ), 1 σ ˆe2 = e T e, n (Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

21 / 42

Modelo Lineal

Modelo de regresi´on lineal - algoritmo MH

Especificaci´on bayesiana yi ∼ N(XiT β, σe2 ) con verosimilitud igual al caso cl´asico. A priori β ∝ 1 y σe2 ∝ 1/σe2 resulta en a posteriori   1 2 T 2 −(n/2+1) f (β, σe |X , Y ) ∝ (σe ) exp − 2 (Y − X β) (Y − X β) . 2σe

(Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

22 / 42

Modelo Lineal

Modelo de regresi´on lineal - algoritmo MH

Especificaci´on bayesiana yi ∼ N(XiT β, σe2 ) con verosimilitud igual al caso cl´asico. A priori β ∝ 1 y σe2 ∝ 1/σe2 resulta en a posteriori   1 2 T 2 −(n/2+1) f (β, σe |X , Y ) ∝ (σe ) exp − 2 (Y − X β) (Y − X β) . 2σe Se puede aplicar el algoritmo MH directamente sobre el par´ametro (β, σe2 ).

(Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

22 / 42

Modelo Lineal

Modelo de regresi´on lineal - muestreo de Gibbs La distribuci´on condicional de σe2 |β f (σe2 |β, X , Y )

∝

(σe2 )(n/2)+1 exp



eT e − 2 2σe



es una gamma inversa con a = n/2 y b = e T e/2.

(Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

23 / 42

Modelo Lineal

Modelo de regresi´on lineal - muestreo de Gibbs La distribuci´on condicional de σe2 |β f (σe2 |β, X , Y )

∝

(σe2 )(n/2)+1 exp



eT e − 2 2σe



es una gamma inversa con a = n/2 y b = e T e/2. La distribuci´on condicional de β|σe2 es proporcional a   1 T exp − 2 (Y − X β) (Y − X β) 2σe y despu´es de una manipulaci´ on matricial se puede expresar como   1 T T T −1 T exp − 2 T −1 [β β − 2β (X X ) (X Y )] 2σe (X X ) completando el cuadrado en β se identifica con una normal. (Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

23 / 42

Modelo Lineal

Modelo de regresi´on lineal - muestreo de Gibbs

El muestreo de Gibbs se aplica a las condicionadas f (σe2 |β, X , Y ) ∼ IG (n/2, e T e/2) f (β|σe2 , X , Y ) ∼ N((X T X )−1 (X T Y ), σe2 (X T X )−1 ) de forma eficiente dado que las distribuciones son de familias f´acilmente muestreables.

(Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

24 / 42

WinBUGS

Programa

1

Muestreo Monte Carlo

2

Modelo Lineal

3

WinBUGS

4

Modelos Jer´arquicos

(Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

25 / 42

WinBUGS

WinBUGS WinBUGS es un software estad´ıstico desarrollado para implementar an´alisis bayesiano y que utiliza m´etodos MCMC para generar muestras de la distribuci´on a posteriori. http://www.mrc-bsu.cam.ac.uk/bugs/winbugs/contents.shtml No olviden instalar la clave de inmortalidad!! Se puede ejecutar WinBUGS desde otros softwares como R, Matlab y Excel. F´acil de usar y flexible, capaz de describir modelos altamente complejos. S´olo hay que especificar el modelo y los datos. MCMC: Metropolis-within-Gibbs −→ Rejection sampling −→ Slice sampling. (Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

26 / 42

WinBUGS

Procedimiento

1

Especificar el modelo

2

Cargar los datos

3

Compilar el modelo y los datos

4

Inicializaci´on: aleatoria o arbitraria

5

Ejecuci´on de las simulaciones y monitoreo de los par´ametros

(Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

27 / 42

WinBUGS

WinBUGS

Instalaci´on

(Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

28 / 42

WinBUGS

WinBUGS

Instalaci´on Men´ u ayuda: manual y ejemplos.

(Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

28 / 42

WinBUGS

WinBUGS

Instalaci´on Men´ u ayuda: manual y ejemplos. Estructura del c´odigo: model{ ...} par´ametros: constantes, nodos estoc´asticos y componentes l´ogicos.

(Univ. Carlos III de Madrid)

Estad´ıstica bayesiana

22-03-11

28 / 42

WinBUGS

WinBUGS

Instalaci´on Men´ u ayuda: manual y ejemplos. Estructura del c´odigo: model{ ...} par´ametros: constantes, nodos estoc´asticos y componentes l´ogicos. Ejemplo: x ∼ N(µ, σ 2 ) −→ x ∼ dnorm(mu, tau) y =x+

z 3

+

1 w

−→ y