Técnicas computacionais em probabilidade e estatística IImbranco/AulaMCMC1_TecnicasII2015.pdf · Uma boa referˆencia introdutoria ´e Chib and Greenberg ... de probabilidade π

Tecnicas computacionais em probabilidade e

estatıstica II

Marcia D’Elia Branco

Universidade de Sao PauloInstituto de Matematica e Estatıstica

http:www.ime.usp.br/ mbranco

Metodos de Monte Carlo baseados em Cadeias deMarkov: Algoritmos de Gibbs e Metropolis-Hastings.

Marcia D’Elia Branco Tecnicas computacionais em probabilidade e estatıstica II

Metodos de Monte Carlo baseados em Cadeias de Markov

O objetivo e construir um algoritmo de simulacao de umaCadeia de Markov (CM) cuja distribuicao limite (estacionaria)e alguma distribuicao de interesse π.

π e multidimensional e usualmente conhecemos somente onucleo desta distribuicao.

O interesse principal e obter amostras das distribuicoesmarginais de π.

Se deixarmos a CM evoluir, a partir de um certo tempo T0

obtemos uma amostra aproximada da distribuicao conjunta πe tambem das distribuicoes marginais de interesse.

Principais algoritmos: Amostrador de Gibbs eMetropolis-Hastings.


O algoritmo de Gibbs

O nome provem do artigo de Geman e Geman (1984)publicado na IEEE, numa aplicacao de restauracao deimagem, onde a distribuicao a posteriori de interesse e adistribuicao de Gibbs.

O algoritmo foi difundido na area de estatıstica a partir doartigo de Gelfand e Smith (1990), os autores percebem ageneralidade do algoritmo GS e comparam com os algoritmosde dados aumentados e SIR.

A construcao do algoritmo depende do conhecimento dasdistribuicoes condicionais completas. A partir da simulacaodestas distribuicoes de forma iterativa, podemos obteramostras aproximadas das distribuicoes marginais.

Uma boa referencia introdutoria e Casella e George (1992) emThe American Statistician.


O algoritmo de Gibbs

Considere o vetor aleatorio θ = (θ1, . . . , θd) com funcao densidade(distribuicao no caso discreto) de probabilidade π e

π(θj | θ1, . . . , θj−1, θj+1, . . . , θd)

a densidade da j-esima componente condicional as demaiscomponentes do vetor, j = 1, 2, . . . , d.

Essas distribuicoes sao denominadas condicionais completas.

Suponha que e possıvel obter um algoritmo de simulacao eficientepara simular amostras dessas distribuicoes.


O algoritmo GS consiste em:

(a) Fixar um conjunto de valores iniciais θ(0)1 , . . . , θ

(0)d .

(b) Obter um novo valor θ(j) = (θ(j)1 , . . . , θ

(j)d ) a partir de um valor

θ(j−1) usando as etapas

θ(j)1 ∼ π(θ1 | θ

(j−1)2 , . . . , θ

(j−1)d )

θ(j)2 ∼ π(θ2 | θ

(j)1 , θ

(j−1)3 , . . . , θ

(j−1)d )

...................................................................

...................................................................

θ(j)d ∼ π(θ1 | θ

(j)1 , θ

(j)2 , . . . , θ

(j)d−1)

(c) Mudar o contador de j para j + 1 e voltar a etapa (b) ateobter a convergencia.


Exemplo basico

Exemplo 1: Casella e George (92)Considere a seguinte distribuicao conjunta

π(x, y) ∝ Cn,xyx+a−1(1− y)n−x+b−1

x = 0, 1, . . . , n e 0 ≤ y ≤ 1.Interesse: obter uma amostra da densidade marginal de x.Considera: (θ1, θ2) = (X,Y ).Note que

x | y ∼ Binomial(n, y)

y | x ∼ Beta(x+ a, n− x+ b)


Exemplo basico

Algortimo GS:

(i) Fixar y(0);(ii) simular x(j) da Binomial(n, y(j−1));(iii) simular y(j) da Beta(x(j) + a, n− x(j) + b);(iv) Fazer j = j + 1 e voltar a (ii).


Ilustracao das trajetorias da CM para o caso d = 2


Aplicacoes na Inferencia Bayesiana

Usualmente em IB e possıvel obter o nucleo da distribuicaoconjunta a posteriori, mas e dificel obter as distribuicoes aposteriori marginais.

No caso das distribuicoes a posteriori condicionais teremformas conhecidas (de facil simulacao), o algoritmo GS podeser utilizado.

Exemplo 2: Carlin, Gelfand e Smith (92)Modelo Poisson com ponto de mudanca.

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

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

As quantidades desconhecidas de interesse sao: θ = (λ, φ,m).



Vamos considerar as seguintes distribuicoes a priori:λ ∼ Gama(a, b) , φ ∼ Gama(c, d) em e uma uniforme em {1, 2, . . . , n}.A distribuicao a posteriori conjunta e proporcional a

λa+Sm−1e−(b+m)λφc+Sn−Sm−1e−(d+n−m)φ

com Sk =k∑

i=1yi .

Condicional a m, λ e φ sao independentes com as respectivasdistribuicoes condicionais

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

φ | m, y ∼ Gama(c+ Sn − Sm, d+ n−m)



Alem disso, condicional a λ e φ a distribuicao a posteriori de m euma discreta em {1, 2, . . . , n} com probabilidades

pm =λa+Sm−1e−(b+m)λφc+Sn−Sm−1e−(d+n−m)φ

n∑

k=1

λa+Sk−1e−(b+k)λφc+Sn−Sk−1e−(d+n−k)φ

Amostras destas distribuicoes sao relativamente simples de seremsimuladas.



O modelo de Poisson com ponto de mudanca foi aplicadopara modelar o numero de desastres em minas de carvao naGra-Bretanha entre os anos de 1851 e 1962.

As n=112 observacoes estao dadas na Tabela 5.1 (retirada deGamerman e Lopes, 2006)

As estatısticas a posteriori obtidas para os parametros domodelo usando o algoritmo GS sao apresentadas na Tabela5.2 (retirada de Gamerman e Lopes, 2006).

A Figura 5.2 apresenta os graficos das distribuicoes aposteriori para os parametros λ, φ e m .

Os resultados indicam um ponto de mudanca no ano de 1890.




O algoritmo de Metropolis-Hastings

Metropolis et.al.(1953), Journal of Chemical Physics.

Hastings (1970), Biometrika, apresenta uma versao mais geraldo algoritmo e introduz o seu uso na area de estatıstica.

O algoritmo depende de um nucleo de transicao proposto Q etambem de uma probabilidade de aceitacao α(θ, φ) para ovalor simulado φ, dado que a cadeia esta em θ.

Uma boa referencia introdutoria e Chib and Greenberg (1995)da revista The American Statistician.


O algoritmo de M-H

Considere o vetor aleatorio θ = (θ1, . . . , θd) com funcao densidade(distribuicao no caso discreto) de probabilidade π e um nucleo detransicao Q associada a uma CM da qual sabemos simular.

(a) Fixar um conjunto de valores iniciais θ(0) = (θ(0)1 , . . . , θ

(0)d ).

(b) Obter um novo valor φ simulado de Q(θ(j−1), .).(c) Calcular a probabilidade de aceitacao

α(θ, φ) = min

{

1,π(φ)Q(φ, θ)

π(θ)Q(θ, φ)

}

.(d) Mover a cadeia para θ(j) = φ se u < α, em que u ∼ U(0,1).

Caso contrario, fazer θ(j) = θ(j−1).(e) Mudar o contador de j para j + 1 e voltar a etapa (b) ateobter a convergencia.


O algoritmo de M-H

O nucleo de transicao da CM gerada pelo algoritmo M-H edado por

P (θ,A) =

∫

A

Q(θ, φ)α(θ, φ)dφ+I(θ∈A)

[

1−

∫

Q(θ, φ)α(θ, φ)dφ

]

Para todo subconjunto A do espaco de estados.

Se Q e irredutıvel e aperiodica e α(θ, φ) > 0, ∀(θ, φ), entao oalgoritmo de M-H define uma CM irredutıvel e aperiodica comtransicao dada por P e distribuicao limite π (Roberts andSmith, 1994).

Usualmente π e Q sao multidimensionais o que pode dificultara simulacao.


Escolha da distribuicao proposta

1. Proposta simetrica: Metropolis.

Q(x, y) = Q(y, x) , ∀x, y ∈ S.

Neste caso

α(θ, φ) = min

{

1,π(φ)

π(θ)

}

.Um exemplo e o uso de um Passeio aleatorio.

θ(j) = θ(j−1) + ε

Se θ esta em Rd e usual considerar ε ∼ Np(0, cV ). Para definir amatriz V podemos utilizar alguma aproximacao da matriz decovariancias a posteriori. O valor c e denominado tuning constant

(constante de afinacao) e pode ser monitorada.


Exemplo 3: Suponha que a distribuicao de interesse seja umamistura de duas normais bivariadas, cuja densidade e dada por

φ(θ) = 0.7fN (θ;µ1,Σ1) + 0.3fN (θ;µ2,Σ2)

com fN a f.d.p. da Normal bivariada,

µ1 =

(

45

)

, µ2 =

(

0.73.5

)

,

Σ1 =

(

1.0 0.70.7 1.0

)

e Σ2 =

(

1.0 −0.7−0.7 1.0

)


Proposta P.A.: q(θ, φ) = fN (φ; θ, cI2) .

Considerou-se dois valores iniciais: (4, 5) e (0, 7). Ainda,c = 0.01 (pequeno), c = 1 (medio) e c = 100 (grande).

Observou-se

O comportamento da CM nao e bom para valores pequeno egrande da constante c.

Valores pequenos de c indicam alta taxa de aceitacao, mas acadeia tem dificuldade e se mover em todo espacoparametrico.

Valores grandes de c indicam baixa taxa de aceitacao,resultando num algoritmo pouco eficiente.

O valor inicial teve pouca influencia no comportamento daCM.



2. Proposta independente.

A distribuicao proposta independe da posicao atual da CM, isto e,Q(θ, φ) = f(φ). Resultando na seguinte probabilidade de aceitacao

α(θ, φ) = min

{

1,w(φ)

w(θ)

}

.w(x) = π(x)

f(x) representa o peso associado ao valor x.

Obs: Se usarmos com proposta a distribuicao a priori, o pesos wserao dados pela razao de verossimilhanca.


Proposta Independente:q(θ, φ) = fN (φ;µ3, cI2) com µ3 = (3.01, 4.55).

Considerou-se dois valores iniciais: (4, 5) e (0, 7). Aindac = 0.5 (pequeno), c = 5 (medio) e c = 50 (grande).

Observou-se

O comportamento da CM nao e bom para valores pequeno egrande da constante c.

Valores pequenos de c indicam alta taxa de aceitacao, mas acadeia tem dificuldade e se mover em todo espacoparametrico.

Valores grandes de c indicam baixa taxa de aceitacao,resultando num algoritmo pouco eficiente.

O valor inicial teve pouca influencia no comportamento daCM.



Tamanho Efetivo da Amostra e Autocorrelacoes

A autocorrelacao de comprimento k (lag k) da CM

ρk =Cov(t(n), t(n+k))

σ2

σ2 e a variancia de t(n) e t e alguma funcao de interesse de θ.Vamos denotar V ar(tn) = τ2n/n, em que tn e a media amostralate instante n. E possıvel mostrar que

τ2n = σ2

(

1 + 2

n−1∑

k=1

n− k

nρk

)

e quando n → ∞ entao τ2n → τ2 com

τ2 = σ2

(

1 + 2∞∑

k=1

ρk

)

se existe a soma da serie.


Tamanho Efetivo da Amostra e Autocorrelacoes

O Tamanho Efetivo da Amostra e dado por

nEF =n

1 + 2∞∑

k=1

ρk

.

Interpretacao: Numa amostra independente de tamanho nEF

obtemos o mesmo valor de variancia para o estimador tn do queaquela obtida pela amostra correlacionada de tamanho n obtida.

Esta quantidade e muito util para comparar algoritmos.

No nosso exemplo n = 5000. Considerando t(θ) = θ2 na melhorsituacao da proposta com Passeio Aleatorio nEF ficou por volta de200. Por outro lado, na melhor situacao da proposta IndependentenEF variou entre 600 e 1200. Evidenciando que neste caso amelhor proposta foi a Independente.




Documents

Técnicas computacionais em probabilidade e estatística IImbranco/AulaMCMC1_TecnicasII2015.pdf · Uma boa referˆencia introdutoria ´e Chib and Greenberg ... de probabilidade π