Métodos estatísticos II
Almir R. Pepato
(Aula preparada com a ajuda daquelas
disponibilizadas por Fred(rik) Ronquist)
Resolução do exemplo numérico
0
1
0
0
0
0
1
0
0
1
0
0
0
1
0
0
Resolução do exemplo numérico
0
1
0
0
0
0
1
0
0
1
0
0
0.00097
0.02828
0.02828
0.00097
0
1
0
0
Resolução do exemplo numérico
0
1
0
0
0
0
1
0
0
1
0
0
0.00097
0.02828
0.02828
0.00097
0.0000026
0
1
0
0
Resolução do exemplo numérico
0
1
0
0
0
0
1
0
0
1
0
0
0.00097
0.02828
0.02828
0.00097
0.0000026
0.0218338
0
1
0
0
Resolução do exemplo numérico
0
1
0
0
0
0
1
0
0
1
0
0
0.00097
0.02828
0.02828
0.00097
0.0000026
0.0218338
0.0000259
0
1
0
0
Resolução do exemplo numérico
0
1
0
0
0
0
1
0
0
1
0
0
0.00097
0.02828
0.02828
0.00097
0.0000026
0.0218338
0.0000259
0.0000026
0
1
0
0
Resolução do exemplo numérico
0
1
0
0
0
0
1
0
0
1
0
0
0.00097
0.02828
0.02828
0.00097
0.0000026
0.0218338
0.0000259
0.0000026
0
1
0
0
Inferência Bayesiana
Exemplo Simples, comparando dois modelos.
Há dois sapos de origami, Joe e Herman. Por
experiências anteriores sabe-se que Joe cai 60%
das vezes em pé, enquanto Herman cai apenas
20% das vezes. O nome dos sapos foi apagado.
Como podemos inferir qual é Joe apenas fazendoos saltar?
Primeiro lançamento, caiu em pé:
Inferência Bayesiana
Segundo lançamento, caiu em pé:
Terceiro lançamento, caiu de costas:
Inferência Bayesiana aplicada à
filogenias
Grupo externo:
Inferência Bayesiana aplicada à
filogenias
C
A
B
Inferência Bayesiana aplicada à
filogenias
Probabilidade
a priori
Probabilidade
Dados
Probabilidade
Probabilidade
a posteriori
Inferência Bayesiana aplicada à
filogenias
Probabilidade posterior
f ( | X )
tree 1
tree 2
Espaço paramétrico
tree 3

Inferência Bayesiana aplicada à
filogenias
D = Dados
 = Parâmetros do modelo
Probabilidade
Posterior
Prior
f ( | D) 
”Verossimilhança”
f ( ) f ( D |  )
 f ( ) f ( D |  ) d
Constante Normalizadora
Monte Carlo-Cadeia de Markov
1-Inicia-se em um ponto arbitrário (θ)
2-Faz-se uma pequena modificação propondo um novo estado (θ*)
3-Calcula-se a razão r entre novo estado θ*, e θ:
(a) r>1: novo estado é aceito.
(b) R<1: novo estado é aceito com uma probabilidade r.
f ( * | D) f ( |  * ) f ( * ) f ( D |  * ) / f ( D) f ( |  * ) f ( * ) f ( D |  * ) f ( |  * )
r






*
*
f ( | D) f ( |  )
f ( ) f ( D |  ) / f ( D) f ( |  ) f ( ) f ( D |  ) f ( * |  )
Monte Carlo-Cadeia de Markov
1-Inicia-se em um ponto arbitrário (θ)
2-Faz-se uma pequena modificação propondo um novo estado (θ*)
3-Calcula-se a razão r entre novo estado θ*, e θ:
(a) r>1: novo estado é aceito.
(b) R<1: novo estado é aceito com uma probabilidade r.
Sempre aceito
2a
1
2b
20 %
tree 1
48 %
tree 2
Aceito às vezes
32 %
tree 3
O tempo que a MCMC passa
amostrando uma região do
espaço paramétrico é uma
estimativa da densidade da
probabilidade posterior naquela
região.
Regulando a cadeia de Markov
• Tipicamente um ou poucos parâmetros são
modificados por vez.
• Uma geração é um ciclo completo ou uma
nova proposta tomada ao acaso.
Novos valores são retirados uniformemente de uma janela de
tamanho δ e centrada em x.
Para lances mais “ousados”: aumente δ, mas isso também
diminuirá as chances de novos estados serem aceitos...
Regulando a cadeia de Markov
”burn-in”
“Mixing”: capacidade da cadeia de explorar
adequadamente as regiões de maior probabilidade
posterior do espaço paramétrico
Não adianta amostrar todas as gerações. As mais próximas estão muito
correlacionadas.
Regulando a cadeia de Markov
Distribuição
esperada
Valores
amostrados
Lances muito acanhados:
taxa de aceitação dos
novos estados altos.
“Mixing” deficiente.
Lances muito ousados:
taxa de aceitação muito
baixa. “Mixing”
deficiente.
Lances “na medida”
Bom “mixing”
Convergência
Convergência é o grau em que a cadeia convergiu para a distribuição de máxima
probabilidade posterior.
Trocando em miúdos: MCMC é uma técnica heurística, precisamos algo que nos dê
segurança a respeito da busca.
Indicadores de convergência:
1- A cadeia atingiu um platô.
2- O comportamento da busca parece adequado:
Através do ESS (Effective Sample Size ):
O número de amostras realmente independentes da distribuição posterior à que a
cadeia de Markov é equivalente.
Convergência
Telas do programa TRACER
Convergência entre corridas
• Topologias:
– Compara as probabilidades dos clados (”split
frequencies”), a diferença entre o desvio padrão das duas
ou mais corridas deve tender a zero.
• Variáveis contínuas
– ”Potential scale reduction factor” (PSRF). Compara
variância dentro e entre as corridas. Deve tender a zero na
medida em que as corridas convergem.
Convergência
Telas do programa AWTY (Are We There Yet)
Comparação das probabilidades posteriores
dos clados de duas corridas.
Esta análise funciona como que parando a
corrida em pontos a intervalos regulares e
verificando as probabilidades posteriores
até aquele ponto.
MC3: Metropolis Coupling Markov
Chain Monte Carlo
T é a temperatura,  é o coeficiente de aquecimento
i  0,1,...,n  1


T

1
/
1


i
A idéia consiste em introduzir uma série de cadeias rodando
em paralelo e acopladas, ou seja, trocando valores entre si.
= 0.2:
Algumas dessas Exemplo
cadeias’para
sãoaquecidas,
isto é: a sua
probabilidade posterior é elevado a um número menor que
i deTprobabilidades
Distr. aparece como que
1. Assim o espaço
1.00
Cadeia fria


0
1
.
00
f

|
X
aplainado.
0.83
1 0.83 f  | X 
0.71
Determinar a melhor
temperatura
Cadeia aquecida
2 0.71
f  | X  é crucial.
3 0.62
f  | X 
0.62
MC3: Metropolis Coupling Markov
Chain Monte Carlo
Cadeia fria
Cadeia aquecida
Cadeia fria
Cadeia aquecida
Cadeia fria
Cadeia aquecida
Cadeia fria
Cadeia aquecida
Cadeia fria
Troca mal sucedida
Cadeia aquecida
Cadeia fria
Cadeia aquecida
Cadeia fria
Cadeia aquecida
Cadeia fria
Troca bem sucedida
Cadeia aquecida
Sumarizando as árvores
• Árvore de Maior Probabilidade Posterior
– Pode ser difícil de encontrar
– Pode ter baixa probabilidade para alguns clados (não reflete suporte)
• Árvore de consenso de Maioria
– Reflete melhor a probabilidade posterior dos clados
– Distribuição de comprimento de ramos pode ser multimodal
• Intervalo de credibilidade de árvores
– Incluí as árvores em ordem decrescente de probabilidade até obter
um intervalo de credibilidade de, e.g., 95 %
Consenso de maioria
Frequências representam
a probabilidade posterior
dos clados
Sumarizando os parâmetros
• Média, mediana, variância são os mais
comuns
• intervalo de credibilidade de 95 %: descarte os
2.5 % superiores e inferiores
• Intervalo de 95 % de maior densidade
posterior: encontre a menor região contendo
95 % da probabilidade posterior
Média e o intervalo de
credibilidade de 95% para
os parâmetros do modelo.
Priors
Antes de falar dos priors é necessário revisar as principais distribuições contínuas
e discretas.
Distribuições contínuas
Distribuições discretas
•
•
•
•
•
•
•
•
•
•
•
Normal
Beta
Gama
Dirichlet
Exponencial
Uniforme
Lognormal
Uniforme
Binomial
Multinomial
Poisson

Distribuição uniforme discreta
Distribuições uniformes são utilizadas quando quer se expressar ausência completa
de conhecimento a respeito de um parâmetro que tem impacto uniforme sobre a
verossimilhança. A uniforme discreta é utilizada para as topologias, por exemplo.
1
2
3
4
5
6
Espaço amostral
  {1, 2, , k }


m(
)





Função da distribuição
Distribuição contínua
Disco com
circumferência 1
Espaço Amostral
(um intervalo)
  0,1
f (x)
Função da densidade
de probabilidades
(e.g. Uniforme (0,1))

a

b
E  a,b
Pr(E) 

 f (x) dx
x E
Evento (um subespaço do
espaço amostral)
Probabilidade
Distribuição exponencial
X~
Lembram dessas equações?
Exp(  )
Parametros:

= taxa de decaimento
f (x)  e x
Média: 1/ 


Nelas percebemos que a
probabilidade, base do calculo
da verossimilhança é uma
função exponencial negativa do
comprimento do ramo. Nada
mais natural portanto que usar
uma distribuição exponencial
para seu prior.
Distribuição Gama
X~
Gamma( ,  )
Parâmetros:
 = formato  = escalar
f (x)  x  1e x
Média:  / 


Gama escalonado:   

Gama escalonado
Como vimos na aula sobre
modelos, a distribuição gama é
utilizada para descrever a
variação na taxa de evolução
entre sítios.
Na verdade, aqui temos um
Hiperprior , isto é, α dita a
distribuição a priori das taxas de
variação e é retirado de uma
distribuição (uniforme por
exemplo) .
Distribuição Beta
X~
Beta( 1, 2 )
Parâmetros:
1, 2 = formato
f (x)  x 1 1 (1 x) 2 1
1 1
Modo
:

  i 1
i

É utilizada para parâmetros que
descrevem proporções de um todo,
com apenas dois eventos possíveis.
Por exemplo: proporção de
invariáveis e razão de
Transversões/Transições.
Distribuição Dirichet
X~
Dir( ) :   1, 2,..., k 
Parâmetros:
 = vetor de k shapes
f (x)   x i i 1
i

Definida
como k proporções de um todo
Semelhante à Beta, mas
para várias classes de
eventos: descreve a
frequência de nucleotídeos
e as taxas no GTR por
exemplo.
Dir(1,1,1,1)
Dir(300,300,300,300)
Porque usar análises Bayesianas
Nós podemos focar em qualquer parâmetro de
interesse (não existem parâmetros “sem uso”)
marginalizando a probabilidade posterior por
sobre outros parâmetros (integrando a incerteza
dos outros parâmetros)
20%
tree 1
48%
tree 2
32%
tree 3
(Porcentagens mostram a probabilidade marginal das árvores)
Porque usar análises Bayesianas
Comprimentos dos ramos
árvores
1



1
2
3
2
3
Probabilidades conjuntas
0.10 0.07 0.12 0.29
0.05 0.22 0.06 0.33
0.05 0.19 0.14 0.38
0.20 0.48 0.32
Probabilidades marginais
Porque usar análises Bayesianas
•Capaz de implementar modelos altamente parametrizados.
•A estimativa da incerteza da árvore e a hipótese filogenética são obtidas ao
mesmo tempo.
•As probabilidades posteriores são de interpretação intuitiva
•Pode incorporar conhecimento prévio a respeito do problema (através do
Prior)
Possível problema
Os Priors!