225
Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers [email protected] Departamento de Matem´ atica Aplicada e Estat´ ıstica Universidade de S˜ ao Paulo Cap. 6 a 10 em Robert & Casella Cap. 5 a 7 em Gamerman & Lopes

Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers [email protected] Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

  • Upload
    others

  • View
    19

  • Download
    0

Embed Size (px)

Citation preview

Page 1: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Monte Carlo via Cadeias de Markov (MCMC)

Ricardo [email protected]

Departamento de Matematica Aplicada e EstatısticaUniversidade de Sao Paulo

Cap. 6 a 10 em Robert & CasellaCap. 5 a 7 em Gamerman & Lopes

Page 2: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Introducao

Page 3: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

The top 10 algorithms with the greatest influence on the

development and practice of science and engineering in the 20th

century (in chronological order):

1. Metropolis Algorithm for Monte Carlo

2. Simplex Method for Linear Programming

3. Krylov Subspace Iteration Methods

4. The Decompositional Approach to Matrix Computations

5. The Fortran Optimizing Compiler

6. QR Algorithm for Computing Eigenvalues

7. Quicksort Algorithm for Sorting

8. Fast Fourier Transform

9. Integer Relation Detection

10. Fast Multipole Method

Dongarra and Sullivan (2000) Guest Editors’ Introduction:

The Top 10 Algorithms, Computing in Science and Engineering,

2, 22-23.

1

Page 4: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Cadeias de Markov

Definicao

Seja X0,X1, . . . um processo estocastico com espaco de estadosfinito ou infinito enumeravel.

Propriedade de Markov: a probabilidade de que a cadeia assumaum certo valor futuro, quando o seu estado atual e conhecido, naose altera se conhecemos o seu comportamento passado.

Em termos probabilisticos, temos que,

P(Xt+1 = y |Xt = x , xt−1, . . . , x0) = P(Xt+1 = y |Xt = x).

2

Page 5: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Objetivo: Gerar valores de uma distribuicao π(·) simulando umacadeia de Markov.

• A cadeia deve ser homogenea, irredutivel e ergodica cujadistribuicao estacionaria seja π(·).

• Dadas as realizacoes {X(t), t = 0, 1, . . . } de uma cadeia deMarkov que tenha π como distribuicao de equilibrio entao, sobas condicoes acima,

X (t) t→∞−→ π(x)

1

n

n∑

t=1

g(X(t)i )

n→∞−→ Eπ(g(Xi ))

3

Page 6: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

• A cadeia e por definicao dependente.

• Media aritmetica dos valores da cadeia e um estimadorconsistente da media teorica.

• Cadeia Irredutivel: probabilidade positiva de atingir qualquerponto a partir de qualquer outro ponto em numero finito deiteracoes.

• Cadeia Aperiodica: nao atinge o mesmo ponto comregularidade fixa.

4

Page 7: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Nucleo de transicao de uma cadeia de Markov.

• No caso discreto,

P(Xt+1 = y |Xt = x).

• No caso continuo, se K (·|x) e o nucleo de transicao entao,

P(X ∈ A|x) =∫

A

K (y |x)dy .

• Iremos considerar cadeias de Markov cujo nucleo de transicaoe a funcao de densidade condicional de Xt+1|Xt .

5

Page 8: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Exemplo. Uma sequencia de variaveis aleatorias {Xt , t = 0, 1, . . . }e dita um passeio aleatorio se,

Xt+1 = Xt + ǫt ,

com ǫt independente de Xt ,Xt−1, . . .

Propriedade Markoviana,

P(Xt+1 ∈ A|x0, x1, . . . , xt) = P(Xt+1 ∈ A|xt)

6

Page 9: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Recorrencia

Em uma cadeia de Markov com espaco de estados finito, umestado e dito recorrente se o numero medio de visitas e infinito.Caso contrario, o estado e dito transiente.

• A cadeia e dita recorrente se o numero medio de visitas a umconjunto qualquer A e infinito.

• Em MCMC estamos interessados em cadeias recorrentes, queexplorem todo o espaco de estados.

Recorrencia de Harris

Se alem de recorrente a cadeia tem Recorrencia de Harris entao,

P(Xt ∈ A|x0) = 1

Ou seja, comecando a cadeia de qualquer ponto inicial (x0) adistribuicao estacionaria e a mesma.

7

Page 10: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Na pratica os metodos MCMC requerem,

• Especificar o nucleo de transicao da cadeia.

• Escolha do valor inicial da cadeia.

• Monitorar a convergencia. Como decidir se a cadeia atingiu oequilibrio?

• A cadeia resultante precisa ser homogenea, irredutivel eaperiodica.

8

Page 11: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Amostrador de Gibbs

O amostrador de Gibbs e uma cadeia de Markov onde as transicoesde um estado para outro sao feitas de acordo com as distribuicoescondicionais completas,

π(xi |x−i ) =π(x)

π(x)dxi

∝ π(x).

sendox−i = (x1, . . . , xi−1, xi+1, . . . , xd)

′.

9

Page 12: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

• Ao inves de gerar valores da distribuicao conjunta π(x), seraogerados valores das distribuicoes condicionais.

• π(xi |x−i ) e obtida considerando os termos da distribuicaoconjunta que nao dependem de xi .

• Gerar valores de π(xi |x−i ), i = 1, . . . , d e equivalente a gerarvalores de π(x) ?

10

Page 13: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

1. inicialize o contador de iteracoes da cadeia t = 0;

2. especifique valores iniciais x(0) = (x(0)1 , . . . , x

(0)d )′;

3. obtenha um novo valor de x(t) a partir de x(t−1) atraves dageracao sucessiva dos valores

x(t)1 ∼ π(x1|x (t−1)

2 , x(t−1)3 , . . . , x

(t−1)d )

x(t)2 ∼ π(x2|x (t)1 , x

(t−1)3 , . . . , x

(t−1)d )

...

x(t)d ∼ π(xd |x (t)1 , x

(t)2 , . . . , x

(t)d−1)

4. Incremente o contador de t para t + 1 e retorne ao passo 3ate obter convergencia.

11

Page 14: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

• Gerar valores das distribuicoes condicionais completas esuficiente para recuperar a distribuicao conjunta.

• A sequencia x(t), t = 1, 2, . . . e uma cadeia de Markov comdistribuicao invariante π(x).

• Cada sequencia x(t)j , t = 1, 2, . . . , j = 1, . . . , d e uma cadeia

de Markov com distribuicao invariante πj(xj).

• As componentes do vetor x podem ser escalares, vetores oumatrizes.

12

Page 15: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Exemplo. Seja um vetor aleatorio X = (X1,X2,X3) com funcaode probabilidade p(x) que admite a seguinte simplificacao,

p(x) ∝ p(x1|x2)p(x2)p(x3).

sendo p(x1|x2), p(x2) e p(x3) conhecidas.

As probabilidades condicionais completas ficam,

p(x1|x2, x3) ∝ p(x1|x2)p(x2|x1, x3) ∝ p(x1|x2)p(x2)p(x3|x1, x2) ∝ p(x3)

13

Page 16: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Exemplo. Sejam X e Y variaveis aleatorias continuas comdensidade de probabilidade conjunta,

f (x , y) = kx4 exp(−x(2 + y)), x > 0, y > 0.

As densidades condicionais completas sao,

f (y |x) ∝ kx4 exp(−x(2 + y)) ∝ exp(−xy), y > 0.

f (x |y) ∝ kx4 exp(−x(2 + y)) ∝ x4 exp(−x(2 + y)), x > 0.

Portanto,

Y |X = x ∼ Exponencial(x),

X |Y = y ∼ Gama(5, 2 + y)

14

Page 17: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

> gibbs <- function(x0,y0,niter=1000) {

+ x = y = rep(0, niter)

+ x[1] = x0

+ y[1] = y0

+ for (i in 2:niter) {

+ x[i] = rgamma(1,5,2+y[i-1])

+ y[i] = rexp(1,x[i])

+ }

+ return(cbind(x,y))

+ }

15

Page 18: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

10 valores gerados com ponto inicial x = 5 e y = 1.5.

x

y

0.02

0.04

0.06

0.08

0.1

0.12

0.14

0.16 0.18

0.2 0.22

0.24 0.26

0 1 2 3 4 5 6

0.0

0.5

1.0

1.5

2.0

16

Page 19: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

1000 valores gerados com ponto inicial x = 5 e y = 1.5.

x

y

0.02

0.04

0.06

0.08

0.1

0.12

0.14

0.16 0.18

0.2 0.22

0.24 0.26

0 1 2 3 4 5 6

0.0

0.5

1.0

1.5

2.0

17

Page 20: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

x y

0 200 400 600 800 1000 0 200 400 600 800 1000

0.0

2.5

5.0

7.5

0

2

4

6

Traços das cadeias simuladas

18

Page 21: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

y

x

0.0 2.5 5.0 7.5 10.0

Distribuições a posteriori com medianas e intervalos de 95%

19

Page 22: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Exemplo. Seja o vetor aleatorio (X ,Y ) com distribuicao normalbivariada,

(X ,Y ) ∼ N

([

00

]

,

[

1 ρρ 1

])

Gerar valores de (X ,Y ) usando as distribuicoes condicionaiscompletas, i.e.

X |Y = y ∼ N(ρy , 1− ρ2)

Y |X = x ∼ N(ρx , 1− ρ2)

Amostrador de Gibbs. Dado um valor yt na iteracao t,

1. gere xt+1|yt ∼ N(ρyt , 1− ρ2),

2. gere yt+1|xt+1 ∼ N(ρxt+1, 1− ρ2).

20

Page 23: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

> gibbs <- function(x0,y0,niter=1000) {

+ x = y = rep(0, niter)

+ x[1] = x0

+ y[1] = y0

+ for (i in 2:niter) {

+ x[i] = rnorm(1,mean=rho*y[i-1],sd=sqrt(1-rho^2))

+ y[i] = rnorm(1,mean=rho*x[i],sd=sqrt(1-rho^2))

+ }

+ return(cbind(x,y))

+ }

21

Page 24: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

x

y

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

2.2

2.4

2.6

2.8

3

3.2

3.4 3

.6

3.8 4

−3 −2 −1 0 1 2 3

−3

−2

−1

01

23

ρ = 0.75

22

Page 25: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

x

y

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8 0.9

1

1.1

1.2

1.3

1.4

1.5

1.6

1.7

−3 −2 −1 0 1 2 3

−3

−2

−1

01

23

ρ = 0.95

23

Page 26: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

x

y

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

2.2

2.4

2.6

2.8

3

3.2

3.4 3

.6

3.8 4

−3 −2 −1 0 1 2 3

−3

−2

−1

01

23

ρ = 0.75

24

Page 27: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

x y

0 200 400 600 800 1000 0 200 400 600 800 1000

−2

0

2

−3

−2

−1

0

1

2

3

Traços das cadeias simuladas

25

Page 28: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

y

x

−2 0 2

Distribuições a posteriori com medianas e intervalos de 95%

26

Page 29: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Exemplo. Modelo autoexponencial (Besag, 1974). Seja o vetoraleatorio y = (y1, y2, y3) cuja funcao de densidade e,

f (y) ∝ exp[−(y1 + y2 + y3 + θ12y1y2 + θ23y2y3 + θ31y3y1)],

com y1, y2, y3 > 0 e θij conhecidos.

Verifique que,

f (y2|y1) ∝ exp[−(y1 + y2 + θ12y1y2)]

1 + θ23y2 + θ31y1,

f (y1) ∝ exp(−y1)

0

exp[−(y2 + θ12y1y2)]

1 + θ23y2 + θ31y1dy2

que sao dificeis de serem simuladas.

27

Page 30: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

No entanto as distribuicoes condicionais completas sao todasexponenciais, por exemplo,

Y3|y1, y2 ∝ Exponencial(1 + θ23y2 + θ31y1).

Verifique as outras!

28

Page 31: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Exemplo. Em um modelo Bayesiano para os dadosy = (y1, . . . , yn) que depende dos parametros θ, λ e δ suponha quea distribuicao conjunta e dada por,

p(y, θ, λ, δ) ∝ p(y|θ, δ) p(θ|λ) p(λ) p(δ).

Apos observar y as distribuicoes de cada parametro dados todos osoutros sao,

π(θ|y, λ, δ) ∝ p(y|θ, δ) p(θ|λ)

π(λ|y, θ, δ) ∝ p(θ|λ) p(λ)

π(δ|y, θ, λ) ∝ p(y|θ, δ) p(δ).

29

Page 32: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Exemplo. Suponha que Y1, . . . ,Yn sao independentes tais que,

Yi ∼ N(µ, σ2), i = 1, . . . , n

µ ∼ N(0, s2)

τ ∼ Gama(a, b)

sendo τ = σ−2 com s2, a e b conhecidos.

A funcao de verossimilhanca e dada por,

p(y|µ, τ) ∝ τn/2 exp

[

−τ

2

n∑

i=1

(yi − µ)2

]

.

30

Page 33: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Portanto, a densidade conjunta de (µ, τ) apos observar y fica,

p(µ, τ |y) ∝ p(y|µ, τ) p(µ) p(τ)

∝ τn/2 exp

[

−τ

2

n∑

i=1

(yi − µ)2

]

exp

[

− µ2

2s2

]

τ a−1e−bτ .

Esta distribuicao conjunta nao tem forma padrao mas asdensidades condicionais completas sao faceis de obter.

31

Page 34: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

p(µ|y, τ) ∝ exp

[

−τ

2

n∑

i=1

(yi − µ)2

]

exp

[

− µ2

2s2

]

∝ exp

[

−1

2(nτ + s−2)µ2 − 2µτ y)

]

∝ exp

[

− 1

2C(µ−m)2

]

sendo C−1 = nτ + s−2 e m = Cnτ y .

p(τ |y, µ) ∝ τ a+n/2−1 exp

[

−τ

(

b +1

2

n∑

i=1

(yi − µ)2

)]

.

32

Page 35: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Portanto, as distribuicoes condicionais completas sao,

µ|y, τ ∼ N(m,C )

τ |y, µ ∼ Gama

(

a+n

2, b +

1

2

n∑

i=1

(yi − µ)2

)

33

Page 36: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

1. inicialize o contador de iteracoes da cadeia t = 0;

2. especifique valores iniciais (µ(0), τ (0));

3. obtenha (µ(t), τ (t)) a partir de (µ(t−1), τ (t−1)) gerandovalores,

µ(t)|τ (t−1), y ∼ N(m(t−1),C (t−1)),

C (t−1) = (nτ (t−1) + s−2)−1,

m(t−1) = C (t−1)nτ (t−1)y

τ (t)|µ(t), y ∼ Gama

(

a+n

2, b +

1

2

n∑

i=1

(yi − µ(t))2

)

4. Incremente o contador de t para t + 1 e retorne ao passo 3ate obter convergencia.

34

Page 37: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

> gibbs1 <- function(y,s2,a,b,mu0,tau0,niter) {

+ n= length(y)

+ tau= rep(0,niter)

+ mu = rep(0,niter)

+ tau[1]= tau0

+ mu [1]= mu0

+ ybar= mean(y)

+ for (i in 2:niter) {

+ C= 1/(n*tau[i-1] + 1/s2)

+ m= C * ybar * n * tau[i-1]

+ mu[i] = rnorm (1, mean=m, sd= sqrt(C) )

+ tau[i]= rgamma(1,a+0.5*n,b+0.5*sum((y-mu[i])^2))

+ }

+ theta = cbind(mu,1/tau)

+ colnames(theta) = c("mu","sigma2")

+ return(theta)

+ }

35

Page 38: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Aplicando a dados simulados Y1, . . . ,Yn ∼ N(2, 4), n = 50 e 2000iteracoes.

> y= rnorm(50,mean=2, sd=2)

Hiperparametros, valores iniciais e numero de iteracoes doamostrador de Gibbs,

> s2 = 4

> a = 0.1

> b = 0.1

> mu0 = 0

> tau0 = 1

> niter= 2000

> m = gibbs1(y,s2,a,b, mu0,tau0,niter)

36

Page 39: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Alguns pacotes do R para analise de cadeias de Markov,

• coda (Convergence Diagnosis and Output Analysis)

• boa (Bayesian Output Analysis)

• ggmcmc (Tools for Analyzing MCMC Simulations fromBayesian Inference)

• bayesplot (Plotting for Bayesian Models)

37

Page 40: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

> library(coda)

> theta = as.mcmc(m)

> summary(theta)

Iterations = 1:2000

Thinning interval = 1

Number of chains = 1

Sample size per chain = 2000

1. Empirical mean and standard deviation for each variable,

plus standard error of the mean:

Mean SD Naive SE Time-series SE

mu 1.363 0.2824 0.006315 0.006315

sigma2 4.091 0.8484 0.018970 0.018581

2. Quantiles for each variable:

2.5% 25% 50% 75% 97.5%

mu 0.8268 1.176 1.362 1.540 1.943

sigma2 2.7074 3.501 3.985 4.597 6.03038

Page 41: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

> plot_title <- ggtitle("Tracos das cadeias simuladas")

> bayesplot_theme_set(theme_minimal())

> color_scheme_set("blue")

> mcmc_trace(as.matrix(m)) + plot_title

39

Page 42: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

mu sigma2

0 500 1000 1500 2000 0 500 1000 1500 2000

2.5

5.0

7.5

0.0

0.5

1.0

1.5

2.0

2.5

Traços das cadeias simuladas

40

Page 43: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

sigma2

mu

0.0 2.5 5.0 7.5

Distribuições a posteriori com medianas e intervalos de 95%

41

Page 44: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

mu sigma2

0 5 10 15 20 0 5 10 15 20

0.0

0.5

1.0

Lag

Au

toco

rre

latio

n

42

Page 45: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Evolucao dos quantis amostrais (2.5%, 50%, 97.5%).

0 500 1500

0.0

0.5

1.0

1.5

2.0

Iterations

mu

0 500 1500

12

34

56

Iterations

sigma2

43

Page 46: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Exemplo. Sejam Y1, . . . ,Yn ∼ Poisson(λ) com λ ∼ Exp(β). Aoinves de especificar um valor para β este tambem recebe umadistribuicao, β ∼ Gama(c , d).

Temos entao que,

p(y|λ) =n∏

i=1

λyi e−λ

yi !

p(λ|β) = βe−βλ

p(β) =dc

Γ(c)βc−1e−dβ

44

Page 47: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Portanto,

p(λ, β|y) ∝n∏

i=1

λyi e−λ

yi !p(λ|β) p(β)

∝ λte−nλ β e−βλ βc−1e−dβ, t =n∑

i=1

yi .

p(λ|β, y) ∝ λte−(β+n)λ

p(β|λ, y) ∝ βce−(λ+d)β

Portanto, as distribuicoes condicionais completas sao

λ|β, y ∼ Gama(t + 1, β + n)

β|λ, y) ∼ Gama(c + 1, λ+ d)

45

Page 48: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

> Gibbs2 <- function(c,d,y,niter){

+ N = length(y)

+ lambda = matrix(0, nrow=niter)

+ beta = matrix(0, nrow=niter)

+ lambda[1]= 1

+ beta [1]= 1

+ t1 = sum(y)

+ for (i in 2:niter) {

+ lambda[i]= rgamma(1, 1+t1, beta[i-1]+N)

+ beta [i]= rgamma(1, 1+c, lambda[i]+d)

+ }

+ return(theta = list(lambda=lambda,beta=beta))

+ }

46

Page 49: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Exemplo. Aplicando o algoritmo anterior com dados simuladosY1, . . . ,Yn ∼ Poisson(4), n = 100, c = d = 0.01 e 2000 iteracoes.

> y= rpois(100,lambda=4)

> niter = 2000

> q = Gibbs2(c = 0.01, d = 0.01, y, niter= niter)

> q1= mcmc(cbind(q$lambda, q$beta))

> q2= window(q1, start=201)

> varnames(q2)= c("lambda","beta")

> y

[1] 6 6 5 2 2 1 5 0 4 5 3 5 3 1 2 5 6 6 1 3

[26] 2 3 4 5 1 5 5 2 3 7 2 3 3 1 1 4 1 3 5 1

[51] 5 3 2 7 4 1 5 6 4 13 4 7 5 1 4 3 3 4 3 5

[76] 2 5 4 3 5 1 7 3 2 3 6 5 6 3 7 4 7 4 2 3

47

Page 50: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Iterations = 201:2000

Thinning interval = 1

Number of chains = 1

Sample size per chain = 1800

1. Empirical mean and standard deviation for each variable,

plus standard error of the mean:

Mean SD Naive SE Time-series SE

lambda 3.7393 0.1878 0.004426 0.003652

beta 0.2817 0.2775 0.006541 0.007670

2. Quantiles for each variable:

2.5% 25% 50% 75% 97.5%

lambda 3.390461 3.61014 3.7403 3.8625 4.113

beta 0.006326 0.08143 0.1972 0.3943 1.019

48

Page 51: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

lambda beta

0 500 1000 1500 0 500 1000 1500

0.0

0.5

1.0

1.5

2.0

2.5

3.5

4.0

Traços das cadeias simuladas

49

Page 52: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

beta

lambda

0 1 2 3 4

Distribuições a posteriori com medianas e intervalos de 95%

50

Page 53: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

lambda beta

0 5 10 15 20 0 5 10 15 20

0.0

0.5

1.0

Lag

Au

toco

rre

latio

n

51

Page 54: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Medias ergodicas ao longo das iteracoes.

3.6

03

.70

3.8

0

lam

bd

a

0.2

00

.30

0.4

0

0 500 1000 1500

be

ta

iteracoes

means

52

Page 55: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Evolucao dos quantis amostrais (2.5%, 50%, 97.5%).

500 1500

3.4

3.6

3.8

4.0

4.2

Iterations

lambda

500 1500

0.0

0.2

0.4

0.6

0.8

1.0

1.2

Iterations

beta

53

Page 56: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Exemplo. Os dados abaixo sao contagens anuais de casamentospor 1000 habitantes na Italia entre 1936 e 1951.

> y = c(7,8,9,7,7,6,6,5,5,7,9,10,8,8,8,7)

As contagens medias nao devem ser constantes ao longo dos anosdevido ao efeito do perıodo de guerra. Vamos assumir entao que,

Yi |λi ∼ Poisson(λi)

λi ∼ Exp(β), i = 1, . . . , n

β ∼ Gama(c , d)

54

Page 57: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

p(β, λ1, . . . , λn | y) ∝n∏

i=1

λyii e

−λi

n∏

i=1

βe−βλi βc−1eβd

p(β|y, λ1, . . . , λn) ∝ βc+n−1e−β(d+∑n

i=1 λi )

p(λi |β, y, λ−i ) ∝ λyii e

−λi (β+1)

Portanto,

β|y, λ1, . . . , λn ∼ Gama(c + n, d +n∑

i=1

λi )

λi |β, y, λ−i ∼ Gama(yi + 1, β + 1)

55

Page 58: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

> Gibbs3 <- function (c, d, y, niter){

+ N = length(y)

+ lambda = matrix(0, nrow=niter, ncol=N)

+ beta = matrix(0, nrow=niter)

+ lambda[1,] = rep(1,N)

+ beta[1] = 1

+ for (i in 2:niter){

+ for (j in 1:N) lambda[i,j]=rgamma(1,y[j]+1,beta[i-1]+1)

+ beta[i] = rgamma(1,c+1, d+sum(lambda[i,]))

+ }

+ return(list(lambda = lambda, beta = beta))

+ }

56

Page 59: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Iterations = 201:2000

Thinning interval = 1

Number of chains = 1

Sample size per chain = 1800

1. Empirical mean and standard deviation for each variable,

plus standard error of the mean:

Mean SD Naive SE Time-series SE

lambda 1 7.998055 2.835081 0.0668235 0.0668235

lambda 2 8.944898 3.120508 0.0735511 0.0735511

lambda 3 9.915417 3.216012 0.0758021 0.0758021

lambda 4 7.957833 2.832384 0.0667599 0.0667599

lambda 5 7.930075 2.799427 0.0659831 0.0659831

lambda 6 6.927021 2.657406 0.0626357 0.0626357

lambda 7 7.014190 2.710355 0.0638837 0.0692631

lambda 8 5.960430 2.408541 0.0567699 0.0567699

lambda 9 5.974878 2.455295 0.0578719 0.0631951

57

Page 60: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

lambda 10 8.001614 2.798360 0.0659580 0.0743087

lambda 11 10.005084 3.177227 0.0748880 0.0748880

lambda 12 10.870986 3.229228 0.0761136 0.0761136

lambda 13 8.862042 2.986458 0.0703915 0.0634856

lambda 14 8.936886 2.984165 0.0703374 0.0625564

lambda 15 8.971401 2.971617 0.0700417 0.0700417

lambda 16 7.849945 2.746720 0.0647408 0.0625234

beta 0.008244 0.007924 0.0001868 0.0001868

2. Quantiles for each variable:

2.5% 25% 50% 75% 97.5%

lambda 1 3.5246480 5.900112 7.726475 9.6493 14.39759

lambda 2 3.9049885 6.710652 8.582586 10.7658 16.06745

lambda 3 4.7946997 7.552836 9.487644 11.8334 17.30265

lambda 4 3.5766547 5.918542 7.534403 9.6779 14.49540

lambda 5 3.4459195 5.918151 7.558955 9.5475 14.33402

lambda 6 2.6414691 4.991869 6.630721 8.3837 13.17473

lambda 7 2.8075146 5.037588 6.600702 8.7135 13.10404

58

Page 61: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

lambda 8 2.2862922 4.169336 5.668539 7.4024 11.62166

lambda 9 2.2348913 4.213837 5.572444 7.3770 11.52448

lambda 10 3.4472962 6.004650 7.680049 9.6216 14.49164

lambda 11 4.8053440 7.745268 9.719038 11.8548 17.36491

lambda 12 5.5583188 8.558385 10.509884 12.8736 17.83346

lambda 13 4.1312262 6.623102 8.578227 10.6630 15.32376

lambda 14 4.0582237 6.780241 8.597522 10.7409 15.43769

lambda 15 4.1557576 6.852114 8.582123 10.7980 15.46722

lambda 16 3.3427246 5.828858 7.542555 9.5362 13.98021

beta 0.0003024 0.002647 0.005813 0.0117 0.02936

59

Page 62: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

media erropadrao

1936 7.998055 2.835081

1937 8.944898 3.120508

1938 9.915417 3.216012

1939 7.957833 2.832384

1940 7.930075 2.799427

1941 6.927021 2.657406

1942 7.014190 2.710355

1943 5.960430 2.408541

1944 5.974878 2.455295

1945 8.001614 2.798360

1946 10.005084 3.177227

1947 10.870986 3.229228

1948 8.862042 2.986458

1949 8.936886 2.984165

1950 8.971401 2.971617

1951 7.849945 2.746720

60

Page 63: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

1940 1945 1950

46

81

01

21

4

61

Page 64: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Boxplots dos valores simulados das variaveis latentes.

1936 1938 1940 1942 1944 1946 1948 1950

05

10

15

20

25

62

Page 65: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Exemplo. Em um processo de contagem no qual foram observadosY1, . . . ,Yn suspeita-se que houve um ponto de mudanca m tal que

Yi |λ ∼ Poisson(λ), i = 1, . . . ,m

Yi |φ ∼ Poisson(φ), i = m + 1, . . . , n.

Deseja-se estimar m, λ e φ.

Distribuicoes a priori independentes:

λ ∼ Gama(a, b)

φ ∼ Gama(c , d)

p(m) = 1/n.

63

Page 66: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

A densidade a posteriori fica,

p(λ, φ,m|y) ∝m∏

i=1

e−λλyi

n∏

i=m+1

e−φφyi λa−1e−bλ φc−1e−dφ 1

n

∝ λa+t1−1e−(b+m)λφc+t2−1e−(d+n−m)φ 1

n

sendo

t1 =m∑

i=1

yi e t2 =n∑

i=m+1

yi .

.

64

Page 67: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Portanto, as densidades condicionais completas ficam,

p(λ|φ,m, y) ∝ λa+t1−1e−(b+m)λ

p(φ|λ,m, y) ∝ φc+t2−1e−(d+n−m)φ

p(m|λ, φ, y) ∝ λt1e−mλφt2e−(n−m)φ, m = 1, . . . , n.

ou seja,

λ|φ,m, y ∼ Gama(a+ t1, b +m)

φ|λ,m, y ∼ Gama(c + t2, d + n −m)

65

Page 68: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

> Gibbs <- function(a,b,c,d,y,niter){

+ N = length(y)

+ lambda = phi = m = matrix(0, nrow=niter)

+ lambda[1] = phi[1] = 1; m[1] = 10

+ for (i in 2:niter) {

+ t1 = sum(y[1:m[i-1]]); t2 = 0; prob = NULL

+ if (m[i-1] < N) t2 = sum(y[(m[i-1]+1):N])

+ lambda[i] = rgamma(1,(a + t1), (b + m[i-1]))

+ phi[i] = rgamma(1,(c + t2), (d + N-m[i-1]))

+ for (j in 1:N){

+ t1 = sum(y[1:j])

+ t2 = 0

+ if (j < N) t2 = sum(y[(j+1):N])

+ aux=(lambda[i-1]^t1)*exp(-j*lambda[i-1])*(phi[i-1]^t2)*exp(

+ prob = c(prob,aux)

+ }

+ soma = sum(prob)

+ probm = prob/soma

+ m[i] = sample(x=N, size=1, prob=probm)

+ }

+ return(list(lambda=lambda, phi=phi, m=m))}66

Page 69: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Exemplo. Testando a funcao Gibbs com 40 dados simulados deprocessos com medias 2 e 5 e ponto de mudanca 23.

> set.seed(124)

> y= c(rpois(n=22, lambda=2),rpois(n=18, lambda=5))

> y

[1] 0 2 2 1 1 1 2 2 4 1 3 3 3 3 2 0 2 3 4 0

[26] 3 7 4 7 4 7 4 4 7 2 10 3 9 3 4

> x = Gibbs(a=0.1, b=0.1, c=0.1, d=0.1, y=y, niter=5000)

67

Page 70: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

> theta = mcmc(cbind(x$lambda,x$phi,x$m))

> theta = window(theta, start=1001)

> colnames(theta) = names(x)

> summary(theta)

Iterations = 1001:5000

Thinning interval = 1

Number of chains = 1

Sample size per chain = 4000

1. Empirical mean and standard deviation for each variable,

plus standard error of the mean:

Mean SD Naive SE Time-series SE

lambda 2.025 0.3036 0.004801 0.005302

phi 5.068 0.5827 0.009214 0.009872

m 23.614 1.6418 0.025959 0.029750

68

Page 71: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

2. Quantiles for each variable:

2.5% 25% 50% 75% 97.5%

lambda 1.470 1.817 2.011 2.218 2.653

phi 3.975 4.667 5.044 5.455 6.239

m 20.000 23.000 23.000 25.000 26.000

69

Page 72: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Probabilidades a posteriori do ponto de mudanca.

> tm= table(theta[,"m"])

> print(tm[1:10])

8 10 16 17 18 19 20 21 22 23

1 1 7 16 11 5 60 150 468 1444

> print(tm[11:length(tm)])

24 25 26 27 28 30

767 410 599 37 23 1

70

Page 73: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

> plot(tm)

05

00

10

00

15

00

tm

8 10 16 18 20 22 24 26 28 30

71

Page 74: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

m

lambda phi

0 1000 2000 3000 4000

0 1000 2000 3000 4000 0 1000 2000 3000 4000

3

4

5

6

7

1.5

2.0

2.5

3.0

3.5

10

15

20

25

30

Traços das cadeias simuladas

72

Page 75: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

phi

lambda

2 4 6 8

Distribuições a posteriori com medianas e intervalos de 95%

73

Page 76: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

lambda phi m

0 5 10 15 20 0 5 10 15 20 0 5 10 15 20

0.0

0.5

1.0

Lag

Au

toco

rre

latio

n

74

Page 77: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Exemplo. Os dados a seguir referem-se ao numero de acidentespor ano em minas de carvao na Inglaterra (acidentes queenvolveram pelo menos 10 mortes) entre 1851 e 1962. Parecehaver uma mudanca em torno do ano 1900. Vamos usar osmesmos passos do exemplo anterior para fazer inferencias.

> yr= 1851:1962

> counts = c(4,5,4,1,0,4,3,4,0,6,3,3,4,0,2,6,3,3,5,4,5,3,1,4,4,1,5,5,3,4

+ 5,2,2,3,4,2,1,3,2,2,1,1,1,1,3,0,0,1,0,1,1,0,0,3,1,0,3,2,2,0,1,

+ 1,1,0,1,0,1,0,0,0,2,1,0,0,0,1,1,0,2,3,3,1,1,2,1,1,1,1,2,4,2,0,

+ 0,0,1,4,0,0,0,1,0,0,0,0,0,1,0,0,1,0,1)

75

Page 78: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

> plot(yr, counts)

1860 1880 1900 1920 1940 1960

01

23

45

6

yr

co

un

ts

76

Page 79: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

> x = Gibbs(a=0.1,b=0.1,c=0.1,d=0.1,y=counts,niter=5000)

> theta = mcmc(cbind(x$lambda,x$phi,x$m))

> theta = window(theta, start=1001)

> colnames(theta) = names(x)

> summary(theta)

Iterations = 1001:5000

Thinning interval = 1

Number of chains = 1

Sample size per chain = 4000

1. Empirical mean and standard deviation for each variable,

plus standard error of the mean:

Mean SD Naive SE Time-series SE

lambda 3.1087 0.2910 0.004601 0.004896

phi 0.9223 0.1155 0.001826 0.001927

m 39.9940 2.4386 0.038558 0.043829

77

Page 80: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

2. Quantiles for each variable:

2.5% 25% 50% 75% 97.5%

lambda 2.5656 2.9009 3.1040 3.301 3.708

phi 0.7086 0.8408 0.9188 0.996 1.160

m 36.0000 39.0000 40.0000 41.000 46.000

78

Page 81: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Probabilidades a posteriori do ponto de mudanca.

1885 1890 1895

02

00

40

06

00

80

01

00

0

yrs

fre

q

79

Page 82: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

m

lambda phi

0 1000 2000 3000 4000

0 1000 2000 3000 4000 0 1000 2000 3000 4000

0.75

1.00

1.25

1.50

2.0

2.5

3.0

3.5

4.0

35

40

45

Traços das cadeias simuladas

80

Page 83: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

phi

lambda

1 2 3 4

Distribuições a posteriori com medianas e intervalos de 95%

81

Page 84: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

lambda phi m

0 5 10 15 20 0 5 10 15 20 0 5 10 15 20

0.0

0.5

1.0

Lag

Au

toco

rre

latio

n

82

Page 85: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Exemplo. Seja X1, . . . ,Xn uma amostra aleatoria da distribuicaot-Student com parametros de locacao θ, escala σ (conhecido) egraus de liberdade ν > 0.

Funcao de verossimilhanca,

p(x|ν, θ, σ2) =n∏

i=1

Γ

(

ν + 1

2

)

Γ(ν

2

)√πν σ

[

1 +(x − θ)2

νσ2

]−(ν+1)/2

.

83

Page 86: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

O modelo pode ser reescrito criando as variaveis latentes (naoobservaveis) λ1, . . . , λn,

Xi |θ, σ, λi ∼ N

(

θ,σ2

λi

)

λi |ν ∼ Gama(ν

2,ν

2

)

para i = 1, . . . , n. A funcao de verossimilhanca completa fica,

f (x,λ|θ, ν) = f (x|θ,λ)f (λ|ν).

84

Page 87: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Especificando distribuicoes a priori, θ ∼ N(µ0, τ20 ) e ν ∼ Exp(β)

tem-se a seguinte densidade a posteriori,

p(θ, ν,λ|x) ∝ f (x|θ,λ)f (λ|ν)f (θ)f (ν)

∝n∏

i=1

(2πσ2/λi )−1/2 exp

{

− λi

2σ2(xi − θ)2

}

n∏

i=1

(ν/2)ν/2

Γ(ν/2)λν/2−1i exp(−νλi/2)

(2πτ20 )−1/2 exp

{

− 1

2τ20(θ − µ0)

2

}

exp(−βν).

85

Page 88: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Densidades condicionais completas,

p(λ|θ, ν, x) ∝n∏

i=1

λ1/2i exp

{

− λi

2σ2(xi − θ)2

}

λν/2−1i exp(−νλi/2)

Portanto,

λi |λ−i , θ, ν, x ∼ Gama

(

ν + 1

2,ν

2+

(xi − θ)2

2σ2

)

86

Page 89: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

p(θ|λ, ν, x) ∝ exp

{

− 1

2σ2

n∑

i=1

λi (xi − θ)2

}

exp

{

− 1

2τ20(θ − µ0)

2

}

∝ exp

{

− 1

2σ2

(

θ2n∑

i=1

λi − 2θn∑

i=1

λixi

)

− 1

2τ20

(

θ2 − 2θµ0

)

}

∝ exp

{

−1

2

[

θ2(σ−2n∑

i=1

λi + τ−20 )− 2θ(σ−2

n∑

i=1

λixi + τ−20 µ0)

]}

Portanto,θ|λ, ν, x ∼ N(µ1, τ

21 )

τ−21 = σ−2

n∑

i=1

λi + τ−20 µ1 =

σ−2∑n

i=1 λixi + τ−20 µ0

σ−2∑n

i=1 λi + τ−20

87

Page 90: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

p(ν|θ,λ, x) ∝n∏

i=1

(ν/2)ν/2

Γ(ν/2)λν/2−1i exp(−νλi/2) exp(−βν).

que nao tem forma funcional conhecida.

88

Page 91: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Exemplo. Testando com 20 dados simulados da t com ν = 6 eσ2 = 1.

> options(width=50)

> round(x,3)

[1] -1.216 3.584 0.700 -1.358 0.850 0.339

[7] -0.034 -0.542 0.009 1.216 0.488 -1.028

[13] 0.982 -1.214 -1.755 0.243 -1.172 -2.216

[19] 2.775 1.008

89

Page 92: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

> mu0 = 0

> tau0= 10

> sigma2= 1

> niter = 5000

> nburn = 1000

> theta = array(0,niter)

> lambda= matrix(0,nrow=niter, ncol= 20)

> theta[1] = mean(x)

> lambda[1,] = 1

> for (i in 2:niter) {

+ s1 = sum(lambda[i-1,] )

+ s2 = sum(lambda[i-1,] * x)

+ tau1= 1 / ((1/tau0) + s1/sigma2)

+ mu1 = (s2/sigma2 + mu0/tau0) * tau1

+ theta[i] = rnorm(1,mean=mu1, sd=sqrt(tau1))

+ for (j in 1:20) {

+ lambda[i,j]=rgamma(1,(nu+1)/2,nu/2+(x[j]-theta[i])^2/(2*sigma2))

+ }

+ }

90

Page 93: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

0 1000 3000

−1

.0−

0.5

0.0

0.5

1.0

Iterations

Trace of var1

−1.0 0.0 1.0

0.0

0.5

1.0

1.5

Density of var1

N = 4000 Bandwidth = 0.05541

91

Page 94: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Boxplots dos valores simulados das variaveis latentes.

1 3 5 7 9 11 13 15 17 19

01

23

45

92

Page 95: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Autocorrelacoes

Os valores simulados segundo um cadeia de Markov sao pordefinicao correlacionados.

Autocorrelacoes amostrais

Considere um sequencia de valores simulados {θ1, . . . , θn}. Asautocorrelacoes de ordem k, ρ(θi , θi+k) podem ser estimadas como,

ρk =

n−k∑

t=1

(θt − θ)(θt−k − θ)

n−k∑

t=1

(θt − θ)2

sendo,

θ =1

n

n∑

t=1

θi .

93

Page 96: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Variancias dos estimadores

Deseja-se estimar uma funcao h(θ) associada a distribuicao alvocom base em uma sequencia de valores simulados {θ1, . . . , θn}.Para uma aproximacao de Monte Carlo,

h =1

n

n∑

i=1

h(θi)

a variancia amostral deve diminuir aumentando-se n.

As autocovariancias amostrais podem ser estimadas como,

γk =1

n

n−k∑

i=1

(h(θi )− h)(θi+k − h).

94

Page 97: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

A variancia de Monte Carlo e dada por,

Var(h) =1

n

(

γ0 + 22δ+1∑

i=1

γi

)

sendo δ o menor inteiro tal que γ2δ + γ2δ+1 > 0.

O tamanho amostral efetivo mede o efeito das autocorrelacoes e edado por,

n =γ0

Var(h).

95

Page 98: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Um metodo simples para estimar o erro padrao

Dados os valores simulados θ1, . . . , θn da distribuicao de equilibrio,

1. Divida a amostra m lotes de tamanho b.

2. Calcule a media hj , j = 1, . . . ,m dentro de cada lote.

3. Estime o erro padrao como,

1

m

1

m − 1

m∑

j=1

(hj − h)2,

sendo h a media global.

96

Page 99: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

• A ideia e considerar as medias dos lotes como sendoindependentes e identidamente distribuidas cujo valoresperado e E (h(θ)).

• O tamanho de cada lote e funcao do numero de valoressimulados. Uma recomendacao pratica e b = ⌊√n⌋.

97

Page 100: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Exemplo. Modelo com variaveis latentes. Lembrando que afuncao de verossimilhanca e dada por,

f (x|θ) =∫

f (x, z|θ)dz

x um vetor aleatorio observavel e z o vetor de variaveis latentes.Alem disso,

f (x, z|θ) = f (x|z, θ)f (z)e

f (z|x, θ) = f (x|z, θ)f (z)f (x|θ) ∝ f (x|z, θ)f (z)

As varıaveis x e z podem ser discretas ou continuas.

98

Page 101: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Se θ tem densidade f (θ) entao pelo teorema de Bayes,

f (θ|x, z) ∝ f (x|z, θ)f (z)f (θ).

Se for possivel simular valores de θ e z destas distribuicoescondicionais completas temos um amostrador de Gibbs.

99

Page 102: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Algoritmos Metropolis-Hastings

Page 103: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Os algoritmos de Metropolis-Hastings usam a mesma ideia dosmetodos de rejeicao. Um valor e gerado de uma distribuicaoauxiliar e aceito com uma dada probabilidade.

Este mecanismo de correcao garante a convergencia da cadeia paraa distribuicao de equilibrio.

100

Page 104: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Definicao

Suponha que a cadeia esteja no estado x e um valor x ′ e gerado deuma distribuicao proposta q(·|x). Note que a distribuicao propostapode depender do estado atual da cadeia, por exemplo q(·|x)poderia ser uma distribuicao normal centrada em x . O novo valorx ′ e aceito com probabilidade

α(x , x ′) = min

{

1,π(x ′) q(x |x ′)π(x) q(x ′|x)

}

. (1)

onde π e a distribuicao de interesse.

101

Page 105: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

• Uma caracterıstica importante e que so precisamos conhecer πparcialmente, i.e. a menos de uma constante ja que nestecaso a probabilidade nao se altera.

• Fundamental em aplicacoes Bayesianas aonde naoconhecemos completamente densidade a posteriori.

• A cadeia pode permanecer no mesmo estado por muitasiteracoes e na pratica costuma-se monitorar isto calculando aproporcao media de iteracoes para as quais novos valores saoaceitos.

102

Page 106: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Em termos praticos, o algoritmo de Metropolis-Hastings pode serespecificado pelos seguintes passos,

1. Inicialize o contador de iteracoes t = 0 e especifique um valorinicial x (0).

2. Gere um novo valor x ′ da distribuicao q(·|x).3. Calcule a probabilidade de aceitacao α(x , x ′) e gere

u ∼ U(0, 1).

4. Se u ≤ α entao aceite o novo valor e faca x (t+1) = x ′, casocontrario rejeite e faca x (t+1) = x (t).

5. Incremente o contador de t para t + 1 e volte ao passo 2.

103

Page 107: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

• A distribuicao proposta possa ser escolhida arbitrariamentemas na pratica deve-se tomar alguns cuidados para garantir aeficiencia do algoritmo.

• Em aplicacoes Bayesianas a distribuicao de interesse e adensidade a posteriori, π = p(θ|x) e a probabilidade deaceitacao assume uma forma particular,

α(θ, θ′) = min

{

1,p(x |θ′)p(x |θ)

p(θ′)

p(θ)

q(θ|θ′)q(θ′|θ)

}

. (2)

O algoritmo sera ilustrado nos exemplos a seguir.

104

Page 108: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Exemplo. Simulando valores de uma distribuicao Beta(0.7,0.7)usando a distribuicao U(0,1) como proposta. Neste caso,

α(x , x ′) = min

{

1,π(x ′) q(x |x ′)π(x) q(x ′|x)

}

= min

{

1,π(x ′)

π(x)

}

sendo π a densidade Beta(0.7,0.7).

105

Page 109: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

> metr0 <- function(niter,a,b,start){

+ x = matrix(NA, nrow=niter)

+ x[1] = start

+ taxa = 0

+ for (i in 2:niter){

+ y = runif(1)

+ A = dbeta(y,a,b)/dbeta(x[i-1],a,b)

+ prob = min(1,A)

+ u = runif(1)

+ if (u < prob) {

+ x[i] = y

+ taxa = taxa + 1

+ }

+ else x[i] = x[i-1]

+ }

+ taxa = taxa/niter

+ return(list(x=x,taxa=round(taxa,2)))

+ }

106

Page 110: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

0.00

0.25

0.50

0.75

1.00

0 1000 2000 3000 4000 5000

va

r1

Traços da cadeia simulada

107

Page 111: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

0.00 0.25 0.50 0.75 1.00

var1

Histograma da cadeia simulada

108

Page 112: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Exemplo. Em uma certa populacao de animais sabe-se que cadaanimal pode pertencer a uma dentre 4 linhagens geneticas comprobabilidades,

p1 =1

2+

θ

4, p2 =

1− θ

4, p3 =

1− θ

4, p4 =

θ

4.

sendo 0 < θ < 1 desconhecido.

Para qualquer θ ∈ (0, 1) e facil verificar que pi > 0, i = 1, 2, 3, 4 e∑4

i=1 pi = 1.

109

Page 113: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Observando-se n animais dentre os quais yi pertencem a linhagemi entao o vetor aleatorio Y = (y1, y2, y3, y4) tem distribuicaomultinomial com parametros n, p1, p2, p3, p4 e portanto,

p(y|θ) =n!

y1!y2!y3!y4!py11 p

y22 p

y33 p

y44

∝ (2 + θ)y1(1− θ)y2+y3θy4 .

110

Page 114: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Atribuindo uma distribuicao θ ∼ U(0, 1) segue que a posteriori eproporcional a expressao acima,

p(θ|y) ∝ p(y|θ) ∝ (2 + θ)y1(1− θ)y2+y3θy4 .

Tomando a distribuicao U(0, 1) como proposta entao q(θ) = 1,∀ θ e a probabilidade de aceitacao se simplifica para,

α(θ, θ′) = min

{

1,p(x |θ′)p(x |θ)

}

= min

{

1,

(

2 + θ′

2 + θ

)y1(

1− θ′

1− θ

)y2+y3(

θ′

θ

)y4}

.

111

Page 115: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

> p <- function(x,y) (2+x)^y[1] * (1-x)^(y[2]+y[3]) * x^y[4]

> metr <- function(n,y,p,start){

+ theta = matrix(NA, nrow=n)

+ theta[1] = start

+ taxa = 0

+ for (i in 2:n){

+ x = runif(1)

+ A = p(x,y)/p(theta[i-1],y)

+ prob = min(1,A)

+ u = runif(1)

+ if (u < prob) {

+ theta[i] = x

+ taxa = taxa + 1

+ }

+ else theta[i] = theta[i-1]

+ }

+ taxa = taxa/n

+ return(list(theta=theta,taxa=round(taxa,2)))

+ }

112

Page 116: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

> m = metr(n=1000,y=c(125,18,20,34),p,start=0.05)

> theta = as.mcmc(m$theta)

> colnames(theta)="theta"

> summary(theta)

Iterations = 1:1000

Thinning interval = 1

Number of chains = 1

Sample size per chain = 1000

1. Empirical mean and standard deviation for each variable,

plus standard error of the mean:

Mean SD Naive SE

0.625057 0.059390 0.001878

Time-series SE

0.004933

2. Quantiles for each variable:

2.5% 25% 50% 75% 97.5%

0.5085 0.5939 0.6261 0.6643 0.7178113

Page 117: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

0.2

0.4

0.6

0 200 400 600 800 1000

the

ta

Traços das cadeias simuladas

114

Page 118: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

0

2

4

6

8

0.2 0.4 0.6

theta

theta

0 5 10 15 20

0.0

0.5

1.0

Lag

Au

toco

rre

latio

n

115

Page 119: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Repetindo com 10000 simulacoes descartando as 5000 primeiras.

0.5

0.6

0.7

0 1000 2000 3000 4000 5000

the

ta

Traços das cadeias simuladas

116

Page 120: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

0

2

4

6

8

0.5 0.6 0.7

theta

theta

0 5 10 15 20

0.0

0.5

1.0

Lag

Au

toco

rre

latio

n

117

Page 121: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Obtendo uma amostra de probabilidades (p1, p2, p3, p4).

> p1 = 1/2+theta/4

> p2 = (1-theta)/4

> p3 = p2

> p4 = theta/4

> prob = as.mcmc(cbind(p1,p2,p3,p4))

> colnames(prob)=c("p1","p2","p3","p4")

118

Page 122: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

> round(summary(prob)$statistics,6)

Mean SD Naive SE Time-series SE

p1 0.655299 0.013002 0.000184 0.000611

p2 0.094701 0.013002 0.000184 0.000611

p3 0.094701 0.013002 0.000184 0.000611

p4 0.155299 0.013002 0.000184 0.000611

> round(summary(prob)$quantile,6)

2.5% 25% 50% 75% 97.5%

p1 0.627293 0.647811 0.65554 0.664693 0.678274

p2 0.071726 0.085307 0.09446 0.102189 0.122707

p3 0.071726 0.085307 0.09446 0.102189 0.122707

p4 0.127293 0.147811 0.15554 0.164693 0.178274

119

Page 123: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Exemplo. No exemplo anterior o algoritmo pode ficar maiseficiente fazendo uma reparametrizacao para a reta. Usando atransformacao logito,

φ = log

(

θ

1− θ

)

∈ R

com transformacao inversa,

θ =exp(φ)

1 + exp(φ).

A distribuicao a priori tambem precisa estar na escalatransformada. Se θ ∼ U(0, 1) a funcao de densidade de φ e,

p(φ) =

=exp(φ)

(1 + exp(φ))2.

120

Page 124: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Valores de φ serao propostos na reta como φ′|φ ∼ N(φ, 1) entao,

q(φ|φ′)

q(φ′|φ) = 1

A probabilidade de aceitacao fica,

α(θ, θ′) = min

{

1,p(x |θ′)p(x |θ)

p(φ′)

p(φ)

}

sendo,

p(x |θ′)p(x |θ) =

[

2 + θ′

2 + θ

]y1[

1− θ′

1− θ

]y2+y3[

θ′

θ

]y4

p(φ′)

p(φ)=

exp(φ′)/(1 + exp(φ′))2

exp(φ)/(1 + exp(φ))2

121

Page 125: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

> prior <- function(phi) exp(phi)/(1+exp(phi))^2

> metr1 <- function(niter,y,p,theta0,tau=1) {

+ phi = matrix(NA, nrow=niter)

+ phi[1] = log(theta0/(1-theta0)); taxa=0

+ for (i in 2:niter) {

+ z = exp(phi[i-1])/(1+exp(phi[i-1]))

+ old = p(z,y)*prior(phi[i-1])

+ x = rnorm(1,mean = phi[i-1], sd=tau)

+ z = exp(x)/(1+exp(x))

+ prob = min(1, p(z,y)*prior(x)/old)

+ u = runif(1)

+ if (u < prob) {

+ phi[i] = x

+ taxa = taxa + 1

+ } else {

+ phi[i] = phi[i-1]

+ }

+ }

+ theta = exp(phi)/(1+exp(phi))

+ taxa = taxa/niter

+ return(list(theta=theta,taxa=taxa))

+ }122

Page 126: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Iterations = 5001:10000

Thinning interval = 1

Number of chains = 1

Sample size per chain = 5000

1. Empirical mean and standard deviation for each variable,

plus standard error of the mean:

Mean SD Naive SE

0.6211405 0.0504233 0.0007131

Time-series SE

0.0016331

2. Quantiles for each variable:

2.5% 25% 50% 75% 97.5%

0.5227 0.5876 0.6220 0.6553 0.7157

123

Page 127: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

0.5

0.6

0.7

0.8

0 1000 2000 3000 4000 5000

the

ta

Traços das cadeias simuladas

124

Page 128: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

0

2

4

6

0.5 0.6 0.7

theta

theta

0 5 10 15 20

0.0

0.5

1.0

Lag

Au

toco

rre

latio

n

125

Page 129: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Iterations = 1:5000

Thinning interval = 1

Number of chains = 1

Sample size per chain = 5000

1. Empirical mean and standard deviation for each variable,

plus standard error of the mean:

Mean SD Naive SE Time-series SE

p1 0.65529 0.01261 0.0001783 0.0004083

p2 0.09471 0.01261 0.0001783 0.0004083

p3 0.09471 0.01261 0.0001783 0.0004083

p4 0.15529 0.01261 0.0001783 0.0004083

2. Quantiles for each variable:

126

Page 130: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

2.5% 25% 50% 75% 97.5%

p1 0.63068 0.64689 0.6555 0.6638 0.6789

p2 0.07108 0.08616 0.0945 0.1031 0.1193

p3 0.07108 0.08616 0.0945 0.1031 0.1193

p4 0.13068 0.14689 0.1555 0.1638 0.1789

127

Page 131: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Exemplo. No exemplo anterior se a distribuicao proposta forφ′|φ ∼ N(φ, τ2) podemos tentar ajustar τ . Para τ = 0.48,

Iterations = 5001:10000

Thinning interval = 1

Number of chains = 1

Sample size per chain = 5000

1. Empirical mean and standard deviation for each variable,

plus standard error of the mean:

Mean SD Naive SE

0.622837 0.051121 0.000723

Time-series SE

0.001686

2. Quantiles for each variable:

2.5% 25% 50% 75% 97.5%

0.5186 0.5898 0.6256 0.6571 0.7186

128

Page 132: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

0.4

0.5

0.6

0.7

0 1000 2000 3000 4000 5000

the

ta

Traços das cadeias simuladas

129

Page 133: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

0

2

4

6

8

0.5 0.6 0.7

theta

theta

0 5 10 15 20

0.0

0.5

1.0

Lag

Au

toco

rre

latio

n

130

Page 134: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Exemplo. Suponha que queremos simular valores X ∼ N(0, 1)propondo valores Y ∼ N(x , σ2). A probabilidade de aceitacao fica,

α(x , y) = min

{

1,fN(y |0, 1)fN(x |0, 1)

fN(x |y , σ2)

fN(y |x , σ2)

}

= min

{

1,fN(y |0, 1)fN(x |0, 1)

}

= min

{

1, exp

(

−1

2(y2 − x2)

)}

.

131

Page 135: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

> metrop <- function(n,sigma){

+ x = matrix(NA,nrow=n)

+ x[1] = taxa = 0

+ for (i in 2:n){

+ y = rnorm(1,x[i-1],sigma)

+ prob = min(1,exp(-0.5*(y^2-x[i-1]^2)))

+ u = runif(1)

+ if (u < prob) {x[i]=y; taxa=taxa+1} else x[i]=x[i-1]

+ }

+ return(list(x=x,taxa=taxa/n))

+ }

132

Page 136: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

xy

0 200 400 600 800 1000

0 200 400 600 800 1000

−2

0

2

−2

−1

0

1

2

sigma = 0.5 e sigma = 10

Traços das cadeias simuladas

133

Page 137: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

• Note que o valor de σ teve um grande impacto na taxa deaceitacao do algoritmo.

• Isto ocorre porque com σ = 0.5 a distribuicao proposta estamuito mais proxima da distribuicao alvo do que com σ = 10.

• Nos dois exemplos anteriores foram ilustrados casos especiaisdo algoritmo nos quais a distribuicao proposta nao dependedo estado atual ou a dependencia e na forma de um passeioaleatorio.

134

Page 138: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

x y

0 5 10 15 20 0 5 10 15 20

0.0

0.5

1.0

Lag

Au

toco

rre

latio

n

Autocorrelações das cadeias simuladas

135

Page 139: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Casos Especiais

Page 140: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Amostrador independente

Um caso particular e quando a distribuicao proposta nao dependedo estado atual da cadeia, i.e. q(x′|x) = q(x′). A probabilidade deaceitacao agora fica,

α(x, x′) = min

{

1,π(x′) q(x)

π(x) q(x′)

}

. (3)

• Em geral, q(·) deve ser uma boa aproximacao de π(·), mas emais seguro se q(·) tiver caudas mais pesadas do que π(·).

• A cadeia nao sera i.i.d. pois a probabilidade de aceitacaoainda depende de x.

136

Page 141: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Teorema

O amostrador independente produz uma cadeia ergodica se existeuma constante c tal que,

f (x)

q(x)≤ c , ∀x .

Lema

Neste caso, a probabilidade de aceitacao esperada e no mınimo1/c se a cadeia for estacionaria.

137

Page 142: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

• O amostrador independente e mais eficiente do que o metodode aceitacao pois em media aceita mais valores propostos.

• No amostrador independente deseja-se uma taxa de aceitacaoalta.

138

Page 143: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Se x ∈ Rd uma distribuicao proposta muita usada e x′ ∼ N(x∗,Σ)

sendo x∗ a moda de π(·) e

Σ = τ

[

−∂2 log π(x)

∂x∂xT

]−1

avaliada em x∗.

Note que a distribuicao proposta e fixa.

139

Page 144: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Algoritmo de Metropolis

Outro caso particular e chamado algoritmo de Metropolis econsidera apenas propostas simetricas, i.e., q(x ′|x) = q(x |x ′) paratodos os valores de x e x ′.

Neste caso a probabilidade de aceitacao se reduz para

α(x , x ′) = min

{

1,π(x ′)

π(x)

}

.

Um algoritmo de Metropolis muito utilizado e baseado em umpasseio aleatorio de modo que a probabilidade da cadeia mover-sede x para x ′ depende apenas da distancia entre eles, i.e.q(x ′|x) = q(|x − x ′|).

140

Page 145: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Neste caso, se usarmos uma distribuicao proposta com variancia σ2

duas situacoes extremas podem ocorrer,

1. se σ2 for muito pequena os valores gerados estarao proximosdo valor atual e quase sempre serao aceitos. Mas levaramuitas iteracoes ate o algoritmo cobrir todo o espaco de x ;

2. valores grandes de σ2 levam a uma taxa de rejeicaoexcessivamente alta e a cadeia se movimenta muito pouco.

Nas duas situacoes o algoritmo fica ineficiente e na pratica temosque tentar varios valores de σ2, monitorando a taxa de aceitacao.

141

Page 146: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Se x ∈ Rd , uma possivel proposta e x′ ∼ N(x,Σ) sendo

Σ = τ

[

−∂2 log π(x)

∂x∂xT

]−1

,

avaliada em x∗.

Note que a distribuicao proposta pode mudar a cada iteracao.

142

Page 147: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Amostragem em Blocos

De um modo geral x = (x1, . . . , xd)′ sera um vetor de parametros

de dimensao d .

Neste caso, pode ser computacionalmente mais eficiente dividir xem k blocos {x1, . . . , xk} e dentro de cada iteracao teremos oalgoritmo aplicado k vezes.

Definimos o vetor x−i = (x1, . . . , xi−1, xi+1, . . . , xk) que contemtodos os elementos de x exceto xi .

143

Page 148: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Suponha que na iteracao t + 1 os blocos 1, 2, . . . , i − 1 ja foramatualizados, i.e.

x−i = (x(t+1)1 , . . . , x

(t+1)i−1 , x

(t)i+1, . . . , x

(t)k ).

Para atualizar a i-esima componente, um valor de xi e gerado dadistribuicao proposta q(·|xi , x−i ) e este valor candidato e aceitocom probabilidade,

α(xi , x′

i ) = min

{

1,π(x′i |x−i ) q(xi |x′i , x−i ))

π(xi |x−i ) q(x′

i |xi , x−i )

}

. (4)

π(xi |x−i ) e a distribuicao condicional completa de xi .

144

Page 149: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Portanto o amostrador de Gibbs e um caso especial do algoritmode Metropolis-Hastings, no qual os elementos de x sao atualizadosum de cada vez (ou em blocos), tomando a distribuicao condicionalcompleta como proposta e probabilidade de aceitacao igual a 1.

145

Page 150: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Exemplo. Seja X1, . . . ,Xn ∼ N(µ, 1/τ) sendo µ ∼ N(0, 1) eτ ∼ Exp(1).

p(x|µ, τ) ∝ τn/2 exp

{

−τ

2

n∑

i=1

(xi − µ)2

}

p(µ) ∝ exp(−µ2/2)

p(τ) = exp(−τ)

p(µ, τ |x) ∝ τn/2 exp

{

−τ

2

n∑

i=1

(xi − µ)2

}

exp(−µ2/2) exp(−τ).

146

Page 151: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Um possivel algoritmo Metropolis-Hastings,

1. gerar (µ′, τ ′),

2. aceitar com probabilidade,

min

{

1,p(x|µ′, τ ′)

p(x|µ, τ)p(µ′) p(τ ′)

p(µ) p(τ)

q(µ, τ |µ′, τ ′)

q(µ′, τ ′|µ, τ)

}

147

Page 152: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Fazendo propostas independentes temos que q(µ, τ) = q(µ)q(τ).

Por exemplo, usando as distribuicoes a priori como propostas,

1. gerar µ′ ∼ N(0, 1) e τ ′ ∼ Exp(1),

2. aceitar com probabilidades,

min

{

1,p(x|µ′, τ ′)

p(x|µ, τ)

}

.

148

Page 153: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

> ll <- function(mu,tau,x) sum (dnorm(x,mean=mu,sd=1/sqrt(tau),lo

> metrop <- function(x,mu0,tau0,s.mu,s.tau,niter) {

+ mu.vec = tau.vec = numeric(niter)

+ mu = mu0; tau = tau0; taxa= 0

+ for (i in 1:niter) {

+ prop.mu = rnorm(1, mean=0, sd=1)

+ prop.tau= rexp (1, rate=1)

+ ratio = ll(prop.mu,prop.tau,x) - ll(mu,tau,x)

+ if (runif(1) < exp(ratio)) {

+ mu = prop.mu

+ tau = prop.tau

+ taxa= taxa + 1

+ }

+ mu.vec [i] = mu

+ tau.vec[i] = tau

+ }

+ cat("taxa de aceitacao:", taxa/niter,"\n")

+ return(cbind(mu=mu.vec,tau=1/tau.vec))}

149

Page 154: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Exemplo. Testando com dados simulados X1, . . . ,Xn ∼ N(2, 4),com n = 50, 10000 iteracoes e 5000 de aquecimento.

> x = rnorm(50, mean=2, sd=2)

> m = metrop(x,mu0=0,tau0=1,s.mu=0.1,s.tau=0.1,niter=10000)

taxa de aceitacao: 0.007

Esta taxa e muito baixa, uma possivel solucao vira no proximo exemplo.

150

Page 155: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Exemplo. No exemplo anterior, reparametrizar θ = (µ, log(τ))sendo

p(log(τ)) = exp(−τ)τ

e a densidade a posteriori e,

p(θ|x) ∝ p(x|µ, τ) p(µ) p(log(τ)).

Propor novos valores como θ′ ∼ N(θ∗,Σ) sendo θ∗ a moda dep(θ|x) e

Σ = τ

[

−∂2 log π(θ)

∂θ∂θT

]−1

avaliada em θ∗.

151

Page 156: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Obter moda e matriz Hessiana avaliada na moda usando a funcaooptim() do R.

> log.posterior <- function(theta,x) {

+ mu = theta[1]

+ tau = exp(theta[2])

+ log.prior = -tau + log(tau) - mu^2/2

+ log.likelihood=(length(x)/2)*log(tau)+sum(-tau*(x-mu)^2/2)

+ log.prior + log.likelihood

+ }

> out=optim(par=c(0,0),fn=log.posterior,control=list(fnscale=-1),

+ hessian=T,x=x)

152

Page 157: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Moda,

[1] 1.975109 -1.039458

matriz Hessiana,

[,1] [,2]

[1,] -18.682323 1.979328

[2,] 1.979328 -26.007433

matriz de variancias,

[,1] [,2]

[1,] 0.053961638 0.004106817

[2,] 0.004106817 0.038763100

153

Page 158: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Distribuicao Normal Multivariada

Seja X um vetor aleatorio d-dimensional com distribuicao normalmultivariada com vetor de medias µ e matriz de variancias ecovariancias Σ. Sua funcao de densidade e dada por,

f (x) =1

(2π)d |Σ|exp

(

−1

2(x− µ)′ Σ−1 (x− µ)

)

, x ∈ Rd

sendo µ ∈ Rd e Σ simetrica e positiva definida, i.e. x′Σx > 0 para

qualquer vetor x 6= 0.

154

Page 159: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Raız Quadrada de uma Matriz

Seja A uma matriz d × d positiva definida. A decomposicao deCholesky encontra uma matriz triangular superior U tal queU ′U = A.

O algoritmo abaixo pode ser utilizado para gerar os valores deX ∼ N(µ,Σ),

1. Calcule U tal que Σ = U ′U.

2. Simule d valores de Zi ∼ N(0, 1), i.e.Z = (Z1, . . . ,Zd)

′ ∼ N(0, Id)

3. Faca Y = U ′Z + µ.

155

Page 160: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

O vetor aleatorio Y tem distribuicao normal multivariada com,

E (Y ) = U ′E (Z) + µ = µ

Var(Y ) = U ′ Var(Z) U = U ′U = Σ.

156

Page 161: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Funcoes do R para simular valores das distribuicoes normal et-Student multivariadas.

> library(mnormt)

> args(rmnorm)

function (n = 1, mean = rep(0, d), varcov, sqrt = NULL)

NULL

> library(MASS)

> args(mvrnorm)

function (n = 1, mu, Sigma, tol = 1e-06, empirical = FALSE, EISPACK = FALSE)

NULL

157

Page 162: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

> metrop1 <- function(x,mu0,tau0,theta,Sigma,niter) {

+ mu.vec = tau.vec = numeric(niter)

+ mu = mu0; tau = log(tau0)

+ for (i in 1:niter) {

+ prop = rmnorm(n = 1, mean=theta,varcov=Sigma)

+ prop.mu = prop[1]

+ prop.tau= prop[2]

+ ratio = log.posterior(c(prop.mu,prop.tau),x)-

+ log.posterior(c(mu,tau),x)

+ if (runif(1) < exp(ratio)) {

+ mu = prop.mu

+ tau= prop.tau

+ }

+ mu.vec[i] = mu

+ tau.vec[i] = exp(tau)

+ }

+ return(cbind(mu=mu.vec,sigma2=1/tau.vec))

+ }

158

Page 163: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

> m=metrop1(x,mu0=0,tau0=1,theta=theta,Sigma=Sigma,niter=10000)

> m=as.mcmc(m[5001:10000,])

> summary(m)

Iterations = 1:5000

Thinning interval = 1

Number of chains = 1

Sample size per chain = 5000

1. Empirical mean and standard deviation for each variable,

plus standard error of the mean:

Mean SD Naive SE Time-series SE

mu 1.976 0.1629 0.002303 0.002944

sigma2 2.880 0.4052 0.005730 0.007093

2. Quantiles for each variable:

2.5% 25% 50% 75% 97.5%

mu 1.651 1.867 1.974 2.087 2.294

sigma2 2.193 2.593 2.837 3.125 3.769159

Page 164: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

mu sigma2

0 1000 2000 3000 4000 5000 0 1000 2000 3000 4000 5000

2

3

4

5

1.6

2.0

2.4

Traços das cadeias simuladas

160

Page 165: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

sigma2

mu

1 2 3 4 5

Distribuições a posteriori com medianas e intervalos de 95%

161

Page 166: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

mu sigma2

0 5 10 15 20 0 5 10 15 20

0.0

0.5

1.0

Lag

Au

toco

rre

latio

n

162

Page 167: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

2

3

4

5

1.6 2.0 2.4

mu

sig

ma

2

163

Page 168: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Exemplo. Seja uma variavel aleatoria X cuja funcao de densidadee,

π(x) =1

2√2π

[

exp

(

−x2

2

)

+ exp

(

−(x − 10)2

2

)]

ou equivalentemente,

X ∼ N(0, 1), com probabilidade 0.5 ou,

X ∼ N(10, 1), com probabilidade 0.5,

164

Page 169: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

> metrop4 <- function(sigma,niter) {

+ u= runif(1,0,1)

+ if (u < 0.5) x=rnorm(1,0,1) else x=rnorm(1,10,1)

+ for(i in 1:(niter-1)) {

+ y=rnorm(1,x[i],sd=sigma)

+ acc=(exp(-y^2/2)+exp(-(y-10)^2/2))/(exp(-x[i]^2/2)+ex

+ u=runif(1,0,1)

+ if (u<acc) x[i+1]=y else x[i+1]=x[i]

+ }

+ return(x)

+ }

165

Page 170: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

> niter=10000

> nburn=5000

> m1= metrop4(sigma=1,niter)

> m1= as.mcmc(m1[(nburn+1):niter])

> m2= metrop4(sigma=10,niter)

> m2= as.mcmc(m2[(nburn+1):niter])

> m = cbind(m1,m2)

166

Page 171: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

m1 m2

0 1000 2000 3000 4000 5000 0 1000 2000 3000 4000 5000

0

5

10

8

10

12

14

Traços das cadeias simuladas

167

Page 172: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

• A primeira cadeia parece estacionaria mas explora somente aregiao em torno da moda 10.

• A segunda cadeia tem taxa de aceitacao menor mas exploraregioes em torno das 2 modas.

168

Page 173: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

m2

m1

−5 0 5 10 15

169

Page 174: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Reparametrizacao

Dados Y1, . . . ,Yn e um parametro θ suponha que serao simuladosvalores de φ = g(θ) com inversa θ = g−1(φ) = h(φ).

A distribuicao alvo e p(φ|y) sendo que,

p(φ|y) = p(θ|y)∣

∂θ

∂φ

= p(θ|y)|J(φ)|

170

Page 175: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Portanto,

α(φ, φ′) = min

{

1,p(φ′|y) q(φ|φ′)

p(φ|y) q(φ′|φ)

}

= min

{

1,p(θ′|y)|J(φ′)| q(φ|φ′)

p(θ|y)|J(φ)| q(φ′|φ)

}

= min

{

1,p(y|θ′)p(θ′) q(φ|φ′)

p(y|θ)p(θ) q(φ′|φ)|J(φ′)||J(φ)|

}

171

Page 176: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Algoritmos de Langevin

Alternativa ao Metropolis-Hastings com passeio aleatorio baseadaem teoria da difusao.

Difusao de Langevin

dLt =1

2∇f (Lt)dt + dBt

sendo Bt um movimento Browniano e f (·) a distribuicaoestacionaria.

A ideia consiste em discretizar este processo e escrever a transicaocomo um tipo de passeio aleatorio.

172

Page 177: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Para um vetor aleatorio x ∈ Rd com densidade f (x),

x(t+1) = x(t) +σ2

2∇x log f (x

(t)) + σǫt (5)

sendo ǫt ∼ N(0, Id) e σ2 o parametro de discretizacao.

173

Page 178: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

E preciso definir uma probabilidade de aceitacao para garantirconvergencia para a distribuicao de equilibrio.

1. Um valor x′ e gerado da distribuicao normal multivariada commedia,

µ(x(t), σ2) = x(t) +σ2

2∇x log f (x

(t))

e matrix de variancia-covariancia σ2Id .

2. Este valor e aceito com probabilidade,

min

{

1,f (x′)

f (x(t))

exp{−||x′ − µ(x(t), σ2)||2/2σ2}exp{−||x(t) − µ(x′, σ2)||2/2σ2}

}

pois a distribuicao proposta e N(µ(x(t), σ2), σ2Id).

174

Page 179: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

• Este algoritmo e conhecido como MALA (Metropolis Adjusted

Langevin Algorithm).

• O uso do gradiente avaliado no valor atual da cadeia deve emprincipio tornar o algoritmo mais eficiente, mas isto nao egarantido na pratica.

• Em cenarios mais complicados pode-se usar derivadasnumericas como aproximacao.

Em aplicacoes Bayesianas,

∇x log f (x) =∂ log f (θ|x)

∂θ

=∂ log f (x|θ)

∂θ+

∂ log p(θ)

∂θ.

175

Page 180: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Exemplo. Suponha que a distribuicao de Y |θ e tal que,

p(y |θ) = f (y , θ)

Z (θ),

sendo Z (θ) uma constante normalizadora,

Z (θ) =

f (y , θ)dy , ou

Z (θ) =∑

y

f (y , θ).

Assume-se que Z (θ) nao tem forma analitica fechada ou nao podeser obtida com recursos computacionais finitos.

A funcao de verossimilhanca e dada por,

p(x|θ) =n∏

i=1

f (xi , θ)

Z (θ),

176

Page 181: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Dada uma distribuicao p(θ) e propondo um novo valor θ′ comdensidade q(·|θ), a probabilidade de aceitacao fica,

α(θ, θ′) = min

{

1,

∏ni=1 f (xi , θ

′)∏n

i=1 f (xi , θ)

p(θ′)

p(θ)

q(θ|θ′)q(θ′|θ)

Zn(θ)

Zn(θ′)

}

.

Note que esta probabilidade e impossivel de ser calculada.

177

Page 182: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Problemas de Dimensao Variavel

Page 183: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

“One of the things we do not know is the number of things we donot know”Peter Green

Em muitas aplicacoes praticas e razoavel assumir que existeincerteza tambem em relacao ao modelo que melhor se ajusta aum conjunto de dados.

• Crie uma variavel aleatoria discreta k (o indicador de modelo)e atribua probabilidades p(k).

• Para cada k existe um vetor de parametros θ(k) ∈ Rnk com

• uma funcao de verossimilhanca: p(y|θ(k), k),

• uma distribuicao a priori: p(θ(k)|k).

178

Page 184: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

A distribuicao de interesse agora e dada por,

π(θ, k|y) ∝ p(y|θ, k) p(θ|k) p(k)

e temos que simular valores desta distribuicao.

• A dimensao de θ pode variar ao longo dos modelos.Precisamos construir uma cadeia com espaco de estados quemuda de dimensao ao longo das iteracoes.

• Os algoritmos de Metropolis-Hastings e o amostrador deGibbs nao podem ser utilizados ja que sao definidos apenaspara distribuicoes com dimensao fixa.

• Embora existam outras possibilidades iremos estudar osalgoritmos MCMC com saltos reversıveis.

179

Page 185: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Exemplo. Sejam Y1, . . . ,Yn os tempos de vida de componenteseletronicos sorteados ao acaso e existe incerteza em relacao adistribuicao dos dados. Sabe-se que

Yi ∼ Exp(λ) (Modelo 1), ou

Yi ∼ Gama(α, β) (Modelo 2),

i = 1, . . . , n.

• k = 1 (Modelo 1), θ(1) = λ,

p(y|λ, k = 1) = λne−λ∑

yi

• k = 2 (Modelo 2), θ(2) = (α, β),

p(y|α, β, k = 2) =βnα

Γn(α)

yα−1i e−β

∑yi .

180

Page 186: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Seja M conjunto de todos os possıveis modelos.

• As probabilidades a posteriori de cada possıvel modelo saodadas por,

π(k|y) = p(k) p(y|k)∑

k∈M

p(k) p(y|k), k ∈ M

• p(y|k) e a verossimilhanca marginal obtida como,

p(y|k) =∫

p(y|θ, k)p(θ|k)dθ.

• Esta ultima integral so e analiticamente tratavel em algunscasos restritos.

• Se o numero de modelos candidatos for muito grande calcular(ou aproximar) p(y|k) pode ser inviavel na pratica.

181

Page 187: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

• Se for especificada a distribuicao de interesse como,

π(θ, k|y) ∝ p(y|θ, k) p(θ|k) p(k)

e conseguirmos simular valores desta distribuicao entaoautomaticamente teremos uma amostra aproximada de π(k|y)e π(θ|k, y).

182

Page 188: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

MCMC com Saltos Reversiveis (RJMCMC)

• Proponha um novo valor para a cadeia e defina umaprobabilidade de aceitacao.

• Os movimentos podem ser entre espacos de dimensoesdiferentes.

• Em cada iteracao atualize os parametros, dado o modelo,usando os metodos MCMC usuais.

• Atualize a dimensao.

183

Page 189: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Suponha que o estado atual da cadeia e (k,θ), i.e. estamos nomodelo k com parametros θ

Um novo modelo k ′ com parametros θ′ e proposto comprobabilidade rk,k ′ . Em geral isto significa incluir ou retirarparametros do modelo atual.

Vamos assumir inicialmente que o modelo proposto tem dimensaomaior, i.e. nk ′ > nk e que θ′ = g(θ,u) para uma funcaodeterministica g e um vetor aleatorio u ∼ q(u) com dimensaonk ′ − nk .

184

Page 190: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Entao o seguinte algoritmo e utilizado,

1. proponha (k,θ) → (k ′,θ′) com probabilidade rk,k ′

2. gere u ∼ q(u) com dimensao nk ′ − nk

3. faca θ′ = g(θ,u),

4. aceite (k ′,θ′) com probabilidade min(1,A) sendo

A =π(k ′,θ′)

π(k,θ)× rk ′,k

rk,k ′ q(u)

∂g(θ,u)

∂(θ,u)

.

185

Page 191: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Exemplo. Sejam Y1, . . . ,Yn os tempos de vida de componenteseletronicos sorteados ao acaso e existe incerteza em relacao adistribuicao dos dados. Sabe-se que,

Yi ∼ Exp(λ) (Modelo 1) ou Yi ∼ Gama(α, β) (Modelo 2),

i = 1, . . . , n. O objetivo e estimar qual modelo explica melhor osdados.

186

Page 192: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Distribuicoes a priori,

p(k) = 1/2

λ|k = 1 ∼ Gama(2, 1)

α|k = 2 ∼ Gama(4, 2)

β|k = 2 ∼ Gama(4, 2)

Funcoes de verossimilhanca,

p(y|λ, k = 1) = λne−λ∑

yi

p(y|α, β, k = 2) =βnα

Γn(α)

yα−1i e−β

∑yi

187

Page 193: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Distribuicoes condicionais completas,

p(λ|y, α, β, k = 1) ∝ p(y|λ, k = 1)p(λ)

∝ λne−λ∑

yiλe−λ

∝ λn+1e−λ(1+∑

yi )

Portanto,

λ|y, α, β, k = 1 ∼ Gama(n + 2, 1 +∑

yi )

p(β|y, α, λ, k = 2) ∝ p(y|α, β, k = 2)p(β)

∝ βnαe−β∑

yiβ3e−2β

∝ βnα+3e−β(2+∑

yi )

Portanto,

β|y, α, λ, k = 2 ∼ Gama(nα+ 4, 2 +∑

yi )

188

Page 194: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

p(α|y, β, λ, k = 2) ∝ βnα

Γn(α)

yα−1i α3e−2α

A distribuicao condicional completa de α nao e conhecida entaovamos usar o algoritmo de Metropolis-Hastings propondo valoresα′ ∼ U[α− ǫ, α+ ǫ].

A probabilidade de aceitacao e,

min

{

1,p(y|α′, β, k = 2)

p(y|α, β, k = 2)

p(α′|k = 2)

p(α|k = 2)

}

ja que q(α′|α) = q(α|α′) = 1/2ǫ.

189

Page 195: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

> mh.alpha <- function(y,n,alpha,beta,eps) {

+ z = runif(1, alpha - eps, alpha + eps)

+ if (z <= 0){

+ acc=0

+ } else {

+ t1=prod(y)

+ num = beta^(n*z) * t1^(z-1)/(gamma(z)^n)

+ den = beta^(n*alpha) * t1^(alpha-1)/(gamma(alpha)^n)

+ num = num * exp(-2*z)*z^3

+ den = den * exp(-2*alpha)*alpha^3

+ }

+ aceita = min(1,num/den)

+ u = runif(1)

+ newalpha = ifelse(u < aceita, z, alpha)

+ return(newalpha)

+ }

190

Page 196: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Suponha que o modelo atual e Exp(λ) e queremos propor omodelo Gama(α, β). Um possivel esquema de atualizacao e oseguine,

1. gere u ∼ Gama(a, b)

2. defina (α, β) = g(λ, u) = (u, λu)

3. calcule o Jacobiano,∣

0 1u λ

= u

4. aceite o novo modelo com probabilidade min(1,A) sendo

A =p(y | α, β, k = 2)

p(y | λ, k = 1)

p(α)p(β)

p(λ)

u

q(u)

191

Page 197: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

• A transformacao no item (2) preserva a media, ou sejaE (Y ) = 1/λ sob o modelo exponencial eE (Y ) = u/λu = 1/λ sob o modelo gama.

• Se o modelo atual for Gama(α, β) e propomos o modeloExp(λ) o esquema reverso consiste em fazer

(λ, u) = g−1(α, β) = (β/α, α).

• A probabilidade de aceitacao e simplesmente min(1, 1/A)substituindo u = α,

A =p(y | λ, k = 1)

p(y | α, β, k = 2)

p(λ)

p(α)p(β)

q(α)

α.

192

Page 198: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

> rj.modelo <- function(y,n,lambda,alpha,beta,model,a,b) {

+ if (model == 1) {

+ u = rgamma(1,a,b)

+ alpha1 = u

+ beta1 = lambda*u

+ lambda1 = lambda

+ } else {

+ lambda1 = beta/alpha

+ alpha1 = alpha

+ beta1 = beta

+ u = alpha

+ }

+ t1 = prod(y); t2 = sum(y)

+ num=beta1^(n*alpha1)*t1^(alpha1-1)*exp(-beta1*t2)/(gamma(alpha1)^n)

+ num=num * 2^4 * alpha1^3 * exp(-2*alpha1)/gamma(4)

+ num=num * 2^4 * beta1^3 * exp(-2* beta1)/gamma(4) * alpha1

+ den=(lambda1^n) * exp(-lambda1*t2)

+ den=den * lambda1 * exp(-lambda1)/gamma(2)

+ den=den * b^a * u^(a-1) * exp(-b*u)/gamma(a)

+ u = runif(1,0,1)

193

Page 199: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

+ if (model == 1) {

+ aceita = min(1,num/den)

+ if (u < aceita) {

+ model = 2

+ alpha = alpha1

+ beta = beta1

+ }

+ } else {

+ aceita = min(1,den/num)

+ if (u < aceita) {

+ model = 1

+ lambda = lambda1

+ }

+ }

+ if (model == 1) return(list(model=model, lambda=lambda))

+ else return(list(model=model, alpha=alpha, beta=beta))

+ }

194

Page 200: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

> rjmcmc <- function(niter,nburn,y,n,a,b,eps=0.25){

+ x = matrix(0, nrow=niter+1, ncol=3)

+ x1 = matrix(0, nrow=niter-nburn, ncol=3)

+ nv = nv1= array(0,2)

+ mod= array(0,niter-nburn)

+ x[1,] = c(1,1,1)

+ model = 1

+ t1 = prod(y)

+ t2 = sum(y)

+ for (i in 1:niter){

+ if (model==1){

+ x[nv[1]+1,1] = rgamma(1, n + 2, t2 + 1)

+ } else {

+ x[nv[2]+1,3] = rgamma(1, 4 + n*x[nv[2],2], t2 + 2)

+ x[nv[2]+1,2] = mh.alpha(y,n,x[nv[2],2],x[nv[2]+1,3],eps)

+ }

+ new = rj.modelo(y,n,x[nv[1]+1,1],x[nv[2]+1,2],x[nv[2]+1,3],model,a,b

+ model = new$model

+ if (i>nburn) mod[i-nburn]= model

+ if (model == 1) {

195

Page 201: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

+ x[nv[1]+1,1] = new$lambda

+ nv[1] = nv[1] + 1

+ if (i > nburn) {

+ x1[nv1[1]+1,1] = new$lambda

+ nv1[1] = nv1[1] + 1

+ }

+ } else {

+ x[nv[2]+1,2] = new$alpha

+ x[nv[2]+1,3] = new$beta

+ nv[2] = nv[2] + 1

+ if (i > nburn) {

+ x1[nv1[2]+1,2] = new$alpha

+ x1[nv1[2]+1,3] = new$beta

+ nv1[2] = nv1[2] + 1

+ }

+ }

+ }

+ cat("Probabilidades a posteriori dos modelos","\n")

+ print(nv1/(niter-nburn))

+ cat("Medias a posteriori dos parametros","\n")

+ somas = apply(x1,2,sum)

196

Page 202: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

+ print(somas/c(nv1[1],nv1[2],nv1[2]))

+ return(list(x=x,nv=nv, x1=x1, nv1=nv1, model=mod))

+ }

197

Page 203: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Exemplo. Testando o algoritmo com saltos reversiveis para oexemplo anterior. Os dados foram simulados como

Y1, . . . ,Yn ∼ Exp(3), com n = 10.

Total de 10000 iteracoes com 5000 de aquecimento.

Distribuicao proposta: u ∼ Gamma(1, 1).

Probabilidades de propor os saltos: rk,k ′ = rk ′,k = 1/2.

198

Page 204: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Probabilidades a posteriori dos modelos

[1] 0.7924 0.2076

Medias a posteriori dos parametros

[1] 3.514171 1.072862 3.306942

O modelo exponencial tem probabilidade a posteriori bem maiorque o modelo Gama.

199

Page 205: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Analise sob o modelo Exponencial.

Iterations = 1:3962

Thinning interval = 1

Number of chains = 1

Sample size per chain = 3962

1. Empirical mean and standard deviation for each variable,

plus standard error of the mean:

Mean SD Naive SE

3.51417 1.00962 0.01604

Time-series SE

0.01604

2. Quantiles for each variable:

2.5% 25% 50% 75% 97.5%

1.797 2.790 3.407 4.136 5.704

200

Page 206: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

2.5

5.0

7.5

0 1000 2000 3000 4000

va

r1

Traço da cadeia simulada

201

Page 207: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

0.0

0.1

0.2

0.3

0.4

2 4 6 8

var1

var1

0 5 10 15 20

0.0

0.5

1.0

Lag

Au

toco

rre

latio

n

202

Page 208: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Analise sob o modelo Gama.

Iterations = 1:1038

Thinning interval = 1

Number of chains = 1

Sample size per chain = 1038

1. Empirical mean and standard deviation for each variable,

plus standard error of the mean:

Mean SD Naive SE Time-series SE

alpha 1.073 0.2977 0.009239 0.009662

beta 3.307 1.0730 0.033304 0.031201

2. Quantiles for each variable:

2.5% 25% 50% 75% 97.5%

alpha 0.5659 0.8551 1.053 1.252 1.740

beta 1.6413 2.5053 3.175 3.955 5.717

203

Page 209: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

alpha beta

0 200 400 600 800 1000 0 200 400 600 800 1000

2

4

6

0.5

1.0

1.5

2.0

Traço da cadeia simulada

204

Page 210: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

beta

alpha

2 4 6

0.5 1.0 1.5 2.00.00

0.25

0.50

0.75

1.00

1.25

0.0

0.1

0.2

0.3

beta

alpha

0 5 10 15 20

0.0

0.5

1.0

0.0

0.5

1.0

Lag

Au

toco

rre

latio

n

205

Page 211: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Estimando os tempos medios de vida (E (Y )) sob o modelo 1(1/λ) e sob o modelo 2 (α/β).

Iterations = 1:3962

Thinning interval = 1

Number of chains = 1

Sample size per chain = 3962

1. Empirical mean and standard deviation for each variable,

plus standard error of the mean:

Mean SD Naive SE

0.310509 0.099184 0.001576

Time-series SE

0.001576

2. Quantiles for each variable:

2.5% 25% 50% 75% 97.5%

0.1753 0.2418 0.2935 0.3584 0.5564

206

Page 212: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Iterations = 1:1038

Thinning interval = 1

Number of chains = 1

Sample size per chain = 1038

1. Empirical mean and standard deviation for each variable,

plus standard error of the mean:

Mean SD Naive SE

0.341420 0.096746 0.003003

Time-series SE

0.002866

2. Quantiles for each variable:

2.5% 25% 50% 75% 97.5%

0.1903 0.2749 0.3265 0.3961 0.5626

207

Page 213: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Outros exemplos

Exemplo. Sejam Y1, . . . ,Yn independentes tais que,

Yi ∼ Bernoulli(pi )

pi = F (α+ βxi), i = 1, . . . , n.

Poderiamos considerar diferentes funcoes de ligacao F (·): logito,probito, Gumbel, etc.

208

Page 214: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Exemplo. Sejam Y1, . . . ,Yn independentes tais que,

Yi ∼ N(µi , σ2), i = 1, . . . , n

sendo,

µi =

β0, (Modelo 0) ou,β0 + β1xi , (Modelo 1) ou,β0 + β1xi + β2x

2i , (Modelo 2)

para uma covariavel x .

209

Page 215: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Exemplo. Sejam Y1, . . . ,Yn independentes tais que,

Yi ∼ N(µi , σ2), i = 1, . . . , n

sendo,

µi = β0 +k∑

j=1

βjxij

Quais covariaveis devem entrar no modelo?

210

Page 216: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Para k = 3 e assumindo que β0 6= 0 temos 23 = 8 possiveismodelos: M0, M1, M2, M3, M12, M13, M23 e M123, sendo M0 omodelo sem convariaveis e M123 o modelo completo.

• Remover uma covariavel consiste em fazer seu coeficienteigual a zero.

• incluir uma covariavel consiste em gerar um novo coeficiente.Por exemplor, podemos propor um salto de M1 para M12

gerando u ∼ N(0, γ2) e fazendo a transformacao,

(β1, β2) = (β1, u)

• Neste caso a funcao g(·) e a identidade e o Jacobiano e iguala 1.

211

Page 217: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Exemplo. Velocidades de 82 galaxias em Km/seg na constelacaode Coroa Boreal.

0 10 20 30 40

0.0

00

.05

0.1

00

.15

0.2

0

velocity of galaxy (1000km/s)

de

nsity

212

Page 218: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

• Podemos considerar o modelo k como uma mistura de k

distribuicoes normais.

• Se X representa a velocidade de uma galaxia em Km/segundo,

f (x |Θ) =k∑

j=1

pj fN(x |µj , σ2j ),

sendo

pj > 0 e

k∑

j=1

pj = 1.

e fN(x |µ, σ2) denota uma densidade normal com media µ evariancia σ2.

• Em geral nao e razoavel assumir um valor fixo para k.

213

Page 219: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Exemplo. Modelo autoregressivo de ordem p.

Xt = θ1pXt−1 + · · ·+ θppXt−p + ǫt

=

p∑

j=1

θjpXt−j + ǫt , ǫt ∼ N(0, σ2p).

O modelo pode ser reescrito como,

p∏

j=1

(1− λjB)Xt = ǫt , ǫt ∼ N(0, σ2p)

sendo λj as raizes inversas (reais ou complexas) e B o operador deretardo.

A condicao de estacionariedade e obtida se |λj | < 1, j = 1, . . . , p.

214

Page 220: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Distribuicao das raizes inversas.

• Se λj ∈ R, λj ∼ U(−1, 1).

• Se λj /∈ R,

λj = a+ bi = r cos θ + (r sin θ)i

sendo,

r ∼ U(0, 1)

θ ∼ U(0, π)

215

Page 221: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

A funcao de densidade conjunta fica,

fp(λ1, . . . , λp) =1

Np

λj∈R

1

2I (|λj | < 1)

λj /∈R

1

πI (|λj | < 1)

sendo Np o numero de diferentes configuracoes das raizes,

Np =⌊p

2

+ 1.

Por exemplo, se p = 2 ha 2 raizes reais ou 2 complexas (Np = 2).Se p = 4 ha 4 raizes reais, ou 2 reais e 2 complexas, ou 4complexas (Np = 3).

Se a cadeia esta no modelo AR(p) e propomos um salto paraAR(p′) os termos Np e Np′ nao se cancelam na probabilidade deaceitacao.

216

Page 222: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Usando todos os modelos

• Assume-se que ∆ e uma quantidade de interesse bem definidaem todos os modelos.

• A distribuicao a posteriori de ∆ e,

p(∆|x) =k∑

i=1

p(∆|Mi , x)p(Mi |x)

com probabilidades a posteriori,

p(Mi |x) =p(x|Mi ) p(Mi)

p(x)

sendo

p(x|Mi) =

p(x|θi ,Mi)p(θi |Mi )dθi

217

Page 223: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

Criterios de Informacao

Criterio de Informacao Deviance (DIC)

DIC (θi ,Mi ) = 2E [D(θi )|x]− D(E [θi |x])= −2 log p(x|θi ,Mi ) + 2pD

= D + pD

sendo

D(θi ) = −2 log p(x|θi ,Mi)

θi = E (θi |x)D = E (D(θi )|x)

O numero efetivo de parametros no modelo e,

pD = D − D(θi ).218

Page 224: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

• Na pratica θi e D sao desconhecidos.

• Se θ1i , . . . ,θ

mi e uma amostra de p(θi |x) entao

D ≈ 1

m

m∑

k=1

D(θki ) e D(θi ) ≈ D

(

1

m

m∑

k=1

θki

)

• O DIC pode ser facilmente calculado durante as simulacoes dacadeia.

219

Page 225: Monte Carlo via Cadeias de Markov (MCMC) · Monte Carlo via Cadeias de Markov (MCMC) Ricardo Ehlers ehlers@icmc.usp.br Departamento de Matem´atica Aplicada e Estat´ıstica Universidade

• O valor individual do DIC nao e relevante, somente diferencassao importantes.

• O DIC nao pode ser calculado se a verossimilhanca dependede algum parametro discreto (e.g. modelos de mistura).

220