Universidade de São Paulo
Escola Superior de Agricultura “Luiz de Queiroz”
Departamento de Ciências Exatas
Ítalo Marcus da Mota Frazão
Luzia Pedroso de Oliveira
Seminário da disciplina de Inferência Bayesiana
Professora: Dra. Roseli Aparecida Leandro
SUMÁRIO
 Introdução
 Proposta do Artigo
 Metropolis Adjusted Langevin Algorithm – MALA
 Hamiltonian Monte Carlo – HMC
 Conceitos Geométricos em MCMC
 Riemann Manifold Metropolis Adjusted Langevin Algorithm – mMALA
 Riemann Manifold Hamiltoniano Monte Carlo, RM – HMC
 Exemplos
INTRODUÇÃO
 Seja
~
p ( )
uma função densidade de probabilidade não-normalizada, com
  D , então a densidade normalizada é dada por:
~
p( ) 
p(  )
~
 p (  ) d
 Acontece que, em muitos problemas estatísticos a integral em não pode ser
resolvida analiticamente.
 Métodos de Monte Carlo para estimativas de integrais são, portanto,
necessários.
 Metodologia predominante para amostragem de uma densidade de
probabilidade é o método de Monte Carlo via Cadeia de Markov (MCMC).
 Metropolis-Hastings: Indiscutivelmente um dos mais influentes e bem
sucedidos algoritmos de Monte Carlo.
Este algoritmo propõe transições    *
então aceitas com probabilidade:
com densidade q( * |  )
que são
 p( * ) q( |  * ) 
 ( , )  min1,

*
p
(

)
q
(

|

)


*
que garante que a cadeia de Markov é reversível com relação à densidade
alvo p ( )
 Em dimensões altas, quando D é grande, o passeio aleatório se torna
ineficiente, resultando em uma baixa taxa de aceitação e amostras
altamente correlacionadas.
 Conseqüência: pequeno tamanho da amostra efetiva da cadeia.
 Uma série de sugestões foram propostas para superar tal ineficiência,
garantindo o equilíbrio e a ergodicidade da cadeia, porém colocam
restrições sobre o que pode ser alcançado em aliviar este problema.
 Grande avanços foram feitos neste sentido com as propostas:
 Metropolis Adjusted Langevian Algorithm (MALA)
 Hamiltonian Monte Carlo method (HMC)
 Apesar dos potenciais ganhos de eficiência obtidos com os mecanismos
MALA e HMC, a afinação desses métodos MCMC continuam a ser uma
grande questão, especialmente desafiadora, para os problemas de
inferência.
 O artigo em questão pretende abordar estas questões de uma forma
sistemática por meio quadro geométrico global para o desenvolvimento de
métodos gerais MCMC.
PROPOSTA DO ARTIGO
 Métodos de amostragem Metropolis Langevin Ajustadas e Monte Carlo
Hamiltoniano definidos no espaço de Riemann para resolver as
deficiências existentes em algoritmos Monte Carlo quando amostrado a
partir de uma densidade alvo que pode ser de alta dimensão e apresentar
fortes correlações
 A metodologia proposta explora a geometria Riemanniana do espaço
paramétrico de modelos estatísticos e assim automaticamente se adapta à
estrutura local quando simulando caminho através deste espaço,
proporcionando a convergência altamente eficiente e exploração da
densidade alvo.
METROPOLIS ADJUSTED LANGEVIN ALGORITHM - MALA
 Processo derivado de uma difusão Langevin discretizada, com um termo
final baseado no gradiente de informação da densidade alvo.
 Considere um vetor aleatório   D com densidade p( ) . Seja então,
L( )  log[ p( )] o logaritmo da densidade.
 O algoritmo Metropolis Langevin Ajustada é baseado em uma difusão
Langevin, com distribuição estacionária p( ) definido por uma equação
diferencial estocástica (EDE), em que a discretização de primeira ordem de
Euler desta EDE dá o seguinte mecanismo proposto:
 n 1   n 
2
2
 L( n )   z n
Em que:
z ~ N (z | 0, I )
e

é a dimensão do passo de integração.
 Problema: A convergência para a distribuição invariante, p ( ) , não é
garantida para
finito, devido ao erro de integração de primeira ordem
introduzido.

 Esta discrepância pode ser corrigida empregando uma probabilidade de
aceitação Metropolis após cada passo de integração, garantindo assim a
convergência.
HAMILTONIAN MONTE CARLO - HMC
 Foi proposto, da mesma forma que no caso MALA, na literatura física
estatística como um meio eficiente de simular estados de um sistema físico
que foi então aplicado a problemas de inferência estatística.
 Consideremos novamente a variável aleatória   D com densidade p( )
D
 Neste método, uma variável auxiliar independente p   é introduzida,
com densidade p( p)  N ( p | 0, M )
 A densidade conjunta segue uma forma fatorada como:
p( , p)  p( ) p( p)  p( ) N ( p | 0, M )
 Denotando o logaritmo da densidade desejada como
L( )  log[ p( )]
 O negativo do logaritmo da verossimilhança conjunta é :
1
1 T 1
D
H ( , p )   L( )  log[( 2 ) M ]  p M p
2
2
H ( , p)
 LHamiltoniano,
( )
 A analogia física de
é um
que descreve a soma de
T
M 1 p
uma função energia ppotencial
definido na posição
e um termo
2
de energia cinética
em que a variável auxiliar p é interpretada
como uma variável cinética e a matriz de covariância M denota uma matriz
de massa.
 O mecanismo de atualização é baseado no integrador Stormer-Verlet ou
Integrador Leapfrog (Duane et al. 1987) e é dado por:
que nada mais é que uma difusão Langevin discreta pré-condicionada, como
empregado em MALA.
 A escolha da matriz de massa M, vai ser crítico para o desempenho do
HMC, porém não existe um princípio orientador de como ela deve ser
escolhida.
CONCEITOS GEOMÉTRICOS EM MCMC
 A relação entre geometria diferencial e estatística tem sido recentemente
empregada no desenvolvimento, principalmente assintótico, da teoria
estatística.
 Conceitos geométricos, como por exemplo, distância e curvatura, são de
interesse natural na metodologia estatística.
 Rao (1945) define, formalmente, a distância entre duas funções densidade.
 Mostrou-se ainda que a matriz de informação de Fisher, G ( ) , é positiva
definida e uma métrica do espaço Riemanniano.
 Portanto, o espaço das funções densidade de probabilidade é dotado de
uma geometria natural de Riemann.
 Ainda em 1945 Rao mostrou que expressões para a curvatura do espaço
entre duas densidades poderia, em principio, ser encontradas.
 Estas idéias foram formalizadas no estudo Information Geometry (Amari
and Nagaoka, 2000)
 Numa perspectiva Bayesiana, a verossimilhança conjunta dos dados e
parâmetros define um tensor métrico, que é a informação de Fischer mais o
Hessiano negativo do logaritmo da distribuição a priori.
 Em suma, o espaço paramétrico de um modelo estatístico é um espaço
Riemanniano. Portanto, uma estrutura geométrica natural da densidade p( )
é definida pelo espaço de Riemann e um tensor métrico associado.
 Então, dada a estrutura geométrica do espaço paramétrico de modelos
estatísticos, a adoção adequada de uma métrica de posição especifica
dentro de um esquema MCMC deve render transições mais eficazes no
algoritmo geral.
RIEMANN MANIFOLD METROPOLIS ADJUSTED LANGEVIN ALGORITHM
 Dada a estrutura geométrica para modelos de probabilidade, uma difusão
Langevin com medida invariante p ( ) ,   D , pode ser definida
diretamente sobre o espaço Riemanniano com tensor métrico G ( )
 A equação diferencial estocástica definindo a difusão Langevin é dada por:
~
1~
d (t )   L[ (t )]dt  d b(t )
2
~
em que, o gradiente natural é  L[ (t )]  G 1 ( (t )) L( (t )) e o
movimento Browniano no espaço Riemanniano é dado por:
 O primeiro termo da equação acima refere-se à curvatura local do espaço e
se reduz a zero se esta curvatura for constante em todos os lugares.
 Algoritmo de atualização:
RM HAMILTONIAN MONTE CARLO
 Aqui o Hamiltoniano que forma a base do HMC será definido de uma
forma geral sobre o espaço de Riemann.
 A definição de um Hamiltoniano em um espaço de Riemann é simples e é
uma técnica empregada em mecanismos geométricos para resolver
equações diferenciais parciais.
 Um Hamiltoniano definido no espaço de Riemann segue como:
e é a base para o método RM HAMILTONIAN MONTE CARLO.
EXEMPLO
 Para N observações retiradas de uma distribuição
métrico baseado na informação de Fisher é :
N 2
G (  ,  )  
 0
N ( x | ,  )
o tensor


2
2N  
0
 Isto define um espaço de Riemann com curvatura constante que é um
espaço hiperbólico na metade superior direita do plano definido por
coordenadas horizontal e vertical  e 
Figura 1: Os contornos acima representam a amostra estimada de
p( ,  | X ) em que a amostra de tamanho N=30 foi
retirada de uma N ( X |   0,  10) .
 Ambas difusões MALA e mMALA foram simuladas a partir de pontos
iniciais 0  0 e  0  40 com tamanho do passo   0.75 para 200
passos .
 A Figura da esquerda, mostra o “caminho” da amostra de um processo
MALA. Como o espaço é hiperbólico e uma métrica Euclidiana é
empregada, a proposta leva à passos ineficientes, de comprimentos quase
que iguais.
 A Figura da direita mostra a proposta mMALA, que é baseada na métrica
do espaço hiperbólico com constante de curvatura negativa e, como tal, as
distâncias percorridas por cada etapa refletem as distâncias naturais neste
espaço, resultando em muito mais eficiência na “travessia” do espaço.
Figura 2: Os mesmos dados amostrais são usados e os valores iniciais
são 0  15 e  0  2 . O tamanho do passo é reduzido para
  0 .2
 A Figuras 3 e 4 fornece uma demonstração intuitiva visual das diferenças
em HMC e RM-HMC quando convergindo para uma amostra de uma
densidade alvo.
Figura 3: Contornos plotados de um modelo de volatilidade estocástica.
Figura 4: Aqui temos uma vista dos “caminhos” da cadeia de Markov. É claro que
RM-HMC efetivamente normaliza os gradientes em cada direção. Entretanto, HMC
apresenta fortes gradientes na direção horizontal comparado com a direção vertical,
logo leva mais tempo para explorar o espaço inteiro.
Exemplo: Modelo Bayesiano de Regressão Logística
Yi ~ Ber(pi)
0 fracasso
yi  
 1 sucesso
P(Yi  y i )  pi i (1 pi )1yi
y
i=1,2,.., n
Função de ligação logística
 p 
i  log it (p i )  log i    0  1x i1  ...   p x ip  Χ iβ
 1 pi 
 pi 
e Χ iβ
  Χ iβ  p i 
 log
1  e Χ iβ
 1 pi 
 x11
x
21

 

 x1n
x 21  x p1 
x 22  x p 2 
   

x 2n  x pn 
Χ i  i - ésima linha da matriz X
Exemplo: Modelo Bayesiano de Regressão Logística
Função de Verossimilhança
n
L(p )
1 y i
i
p
(
1

p
)
 i
i
y
i1
Χ iβ
 e
L(β )   
Χ iβ
i1  1  e
n
Χ iβ
 e
  
Χ iβ
1

e
i1 
n



yi
Χ iβ

e
1 
Χ iβ
 1 e
yi
  1 
 

Χ iβ
  1 e 
n
  e Χ iβyi (1  e Χ iβ ) 1
i1
n
e
 Χ iβyi
i 1
n
i
 (1  e
i1
Χ iβ 1
)



1 y i
1 y i
Exemplo: Modelo Bayesiano de Regressão Logística
Distribuição Conjunta a Priori
β ~ N( 0, Ι)
1
(β )

( 2 )
 e
p 1
2
Ι
1
2
e
1
-1
 β T ( Ι )
β
2
1
-1
 β T ( Ι )
β
2
e

1 T
β β
2
Distribuição Conjunta a Posteriori
n
 Χ iβyi
(β | dados )  e i 1
n
i
 (1  e
Χ iβ 1
) e

1 T
β β
2
i1
log (β | dados ) 
expressão do artigo
n
n
Χβ
 Χ iβy i   log(1  e i ) i 
i 1
i1
=
1 T
β β
2
Que corresponde a
expressão do artigo
Exemplo: Modelo Bayesiano de Regressão Logística
(log joint likelihood?)
X matriz do delineamento de dimensão (
)?
linha da matriz X
N tamanho da amostra
D no de covariáveis
t variável resposta binária com distribuição Bernoulli
função de ligação: logística
com priori apropriada (no ex N(0,) ) com  conhecido
Termos necessários para a aplicação dos métodos
RM-HMC e mMala
• derivada de
:
• 2a. derivada de
:
• tensor métrico: matriz de informação de Fisher + negativo da
Hessiana do logaritmo da priori
em que a matriz
tem elementos
• derivadas do tensor métrico
em que a matriz diagonal
tem elementos
Exemplo: Modelo Bayesiano de Regressão Logística
• São analisados 6 conjuntos de dados (Michie et al, 1994;
Ripley, 1996), com características bem abrangentes
Desvios padrões das distribuições marginais a posteriori variando
de 0,0004 a 9,9
Exemplo: Modelo Bayesiano de Regressão Logística
• Objetivo: obter para cada conjunto de dados uma amostra da
distribuição conjunta a posteriori (o, 1, 2,… k| dados), em
que o é o intercepto e i , i=1,…k são os coeficientes das variáveis
regressoras no modelo logístico com variável resposta Bernoulli.
• Para os conjuntos de dados, exceto Ripley, um modelo de regressão
linear logístico com intercepto foi ajustado. No caso do Ripley
foram considerados também os termos quadrático e cúbico.
• Utilizou-se as prioris i ~ N(0,100)
Exemplo: Modelo Bayesiano de Regressão Logística
• Para os 6 conjuntos de dados são comparados os
métodos RM-HMC e mMALA propostos no artigo,
juntamente com alguns já existentes.
Exemplo: Modelo Bayesiano de Regressão Logística
• Em cada conjunto de dados os métodos foram rodados 10
vezes e os resultados médios foram gravados.
• Os autores reproduziram os resultados de Holmes and
Held(2005), considerando um burning de 5000. As próximas
500o iterações foram utilizadas.
• O tempo de CPU necessário para obter as amostras foi
também registrado para cada caso.
• Os métodos foram implementados no MATLAB.
• Os autores comparam a eficiência desses métodos através dos
tamanhos efetivos das amostras (ESS) das posterioris marginais.
soma das autocorrelações amostrais.
Ex: cadeia gerada pelo método M_H
Dados Heart, “Trace plot” n=1000 e gráficos autocorrelacões
para um dos parâmetros
com média ao redor de -7
Métodos
Metropolis
IWLS
Aux. Var.
HMC
mMALA
Simplified mMala
RM-HMC
REFERÊNCIAS
 Duane. S. Kennedy. A. D. Pendleton. B. J. and Roweth. D. (1987) Hybrid
Monte Carlo, Physics Letters. B., 55, pp. 2774–2777.
 Rao. C. R. Information and Accuracy Attainable in the Estimation of
Statistical Parameters. Bulletin of the Calcutta Mathematical Society. 37.
pp 81 – 91.
 Amari. S. and Nagaoka. H. (2000) Methods of Information Geometry,
Oxford University Press
• Michie, D., Spiegelhalter, D. J., and Taylor, C. C. (1994). Machine
Learning, Neural and Statistical Classification. Prentice Hall,
Englewood Cliffs, N.J.
• Ripley. B. D. (1996) Pattern Recognition and Neural Networks.
Cambridge University press.
• Seefeld, K. Linder, E. Statistics Using R with Biological Examples.
2007