Um exemplo do algoritimo de Metropolis–Hastings aplicado a distribuição binomial

E la vamos nos denovo olhar o algoritmo Metropolis–Hastings

Bem depois de olhar o exemplo do como estimar o valor de pi usando o método de monte carlo, devemos progredir um pouco no entendimento desse negócio de mcmc.

Ja tivemos posts aqui e aqui, onde olhamos como o algoritmo funcionava para fazer uma regressão linear.

Mas vamos tentar um exemplo mais simples para estimar a probabilidade encontrar ou não aranhas num arbusto.

Certo então existe uma probabilidade \theta (letra grega theta) na qual eu encontro ou não as aranhas nos arbustos.

Então eu fui la e olhei 14 arbustinhos e vi o seguinte:

> meusdados [1] 0 1 1 0 1 0 0 1 0 1 1 1 1 1

Ok 1 significam que nesse arbusto eu encontrei uma aranha, e 0 significa que não tinha aranha.

Eu posso estimar a verosimilhança para várias probabilidades, para vários valores de \theta, ou mais precisamente, o p(D|\theta) onde D são meus dados

likelihood <- function( theta , data ) {
  z <- sum( data == 1 )
  N <- length( data )
  pDataGivenTheta <- theta^z * (1-theta)^(N-z)
  pDataGivenTheta[ theta > 1 | theta < 0 ] <- 0
  return( pDataGivenTheta )
}

Ou seja qual a chance de determinado \theta dado meus dados (p(D|\theta)). Aqui tem uma revisão rápida de probabilidade condicional.

Bom mais fácil que ficar nisso é fazer um gráfico para entender o que essa função ta fazendo.

01

Então temos a como são as chances dado um valor de \theta, ou seja olhando os dados la tivemos 9 arbustos com aranhas, dos 14 que olhamos, sabemos que o maximum likelihood nesse caso é 9/14, ou seja o pico do likelihood é quanto \theta esta num valor igual a 0.64, ou seja um \theta de 0.64 é o mais provável de ter produzido esses resultado, dado os dados que temos, veja que se você pegar um valor de vamos supor 0.2, quase não tem como achar 9 arbustos com aranhas de 14 arbustos verificados, teria que ser uma cagada muito grande.

Mas lembramos que em inferência bayesiana, ponderamos nossos dados, ou seja multiplicamos esse likelihood por um prior.

Então fazemos uma função para gerar esse prior.

prior <- function( theta ) {
  prior <- dbeta(theta,1,1)
  #exemplo de um prior bimodal
  #prior < dbeta( pmin(2*theta,2*(1-theta)),2,2 )
  prior[ theta > 1 | theta < 0 ] <- 0
  return( prior )
}

Aqui, o prior ta definido ali como uma distribuição beta com parâmetros a=1 e b=1. Isso vai dar uma gráfico assim:

02

Veja que todas as possibilidades estão em 1, ou seja eu estou atribuindo uma chance igual para todos os valores possíveis de \theta.

Veja que podemos mudar nossa crença inicial mudando os parâmetros a e b para outros valores, se você não sabe o que acontece com diferentes valores, é mais ou menos isso, veja alguns exemplos:

01

Certo, mas para facilidade do exemplo vamos assumir aquele prior e não ficar discutindo muito isso agora.

Agora se lembrarmos o teorema de bayes usado na inferência bayesiana:

p(\theta|D) = {p(D|\theta) * p(\theta) \over {\int p(D|\theta)p(\theta)}}

Então o que falta é sabe qual o valor desse divisor. Na verdade para esse exemplo da pra resolver analiticamente como feito aqui ou aproximar como aqui, mas vamos usar o algorítimo de Metropolis–Hastings para construir a probabilidade posterior.

A ideia é a mesma de quando estimamos o valor de pi. Relembrando os passos:

1. Defina um domínio de todas as possibilidades 2. Gere aleatoriamente entradas de uma distribuição de probabilidade que cobre o domínio. 3. Faça uma computação determinística das entradas 4. Reuna os resultados

Então o domínio são as probabilidades de achar aranha em um arbusto, que são valores entre 0 e 1.

Eu tenho que gerar valores ao acaso de probabilidade, que cubram todas as possibilidades e fazer uma computação determinística para levar ele em conta ou não, no caso do pi o que a gente gerava um ponto ao acaso dentro de um quadradinho e tínhamos um círculo dentro desse quadrado, ai olhávamos se ele estava ou não dentro do círculo, para no final fazer uma razão.

Aqui o que temos são valores possíveis de \theta e eu quero saber qual gerou os dados, bom isso eu nunca vou saber, mas eu quero ter uma estimativa, uma distribuição que se o valor de \theta seja bem razoável, ele tem que ser alto, senão ele é baixo.
Talvez eu esteja falando bobeira aqui, mas o que acontece é o seguinte, no caso do Metropolis–Hastings, ficamos escolhendo um valor de \theta ao acaso, e colocamos ele ou não na nossa cadeia numa probabilidade que depende de o quão bom o novo \theta é dividido pelo \theta anterior.

É difícil explicar isso, acho que entender também, mas fazemos uma função assim:

targetRelProb <- function( theta , data ) {
  targetRelProb <-  likelihood( theta , data ) * prior( theta )
  return( targetRelProb )
}

Essa função multiplica o likelihood pelo prior, dai geramos uma valor de theta vai ter essa multiplicação do valor novo de theta pela multiplicação com o antigo.

Vamos olhar alguns exemplos, bem eu usei 0.7 como chance para gerar os dados. Então vamos supor que o valor anterior na cadeia é 0.7 e fosse sugerido fosse de 0.4, qual a chance de aceitar esse valor e colocar na cadeia?

> targetRelProb(0.4,meusdados)/targetRelProb(0.7,meusdados) [1] 0.2078775

20%, baixinho certo, ou seja eu quero “visitar” esse valor pouco, quero poucas ocorrências dele na cadeia, assim no histograma da distribuição posterior ele vai ter uma barrinha baixinha, mas impossível ver um valor desse? Não uai, mas é pouco provável.

Agora vamos mais pra perto, se propusermos um valor de 0.5.

> targetRelProb(0.5,meusdados)/targetRelProb(0.7,meusdados) [1] 0.6224313

Olha ai, temos 60% de chance de enfiar esse valor na cadeia, porque ta mais pertinho do valor real. E quanto mais próximo vamos do valor real, agora com 0.55.

> targetRelProb(0.55,meusdados)/targetRelProb(0.7,meusdados) [1] 0.8666388

Temos mais e mais chance de enfiar esses valores pra dentro da cadeia.
Agora até aqui eu sempre tava dividindo com um theta de 0.7, mas e se sei la estamos no começo da cadeia, e estamos em valores bem distantes desses. Por exemplo estamos em 0.2 e temos uma proposta de ir para 0.1

> targetRelProb(0.1,meusdados)/targetRelProb(0.2,meusdados) [1] 0.003519595

Meu, 0.2 ta longe pra caramba das probabilidade real, mas 0.1 é muito mais longe, então a chance de aceitar esse valor é muito muito baixa. Agora num caso ao contrário, estamos em 0.9, que é longe e recebemos uma proposta de 0.8, que é bem mais próximo, qual a chance?

> targetRelProb(0.8,meusdados)/targetRelProb(0.9,meusdados) [1] 11.08606

Puxa, um valor muito maior que 1, veja que isso pode acontecer, na verdade aqui o que fazemos e com certeza aceitar o valor proposto, como se estamos em 0.5 e recebemos a proposta de 0.6

> targetRelProb(0.6,meusdados)/targetRelProb(0.5,meusdados) [1] 1.690757

Mas para fugir desses valores e ficar com valores sempre entre 0 e 1 na cadeia temos que fazer umas correções nos valores, mas basicamente é isso o algoritimo, gera um theta novo e aceita ele ou não, um bocado de vezes, mas um bocado são muitas e muitas vezes mesmo.

Então o algoritimo vai ser isso aqui:

for ( t in 1:(Ncadeia-1) ) {
    valoratual = cadeia[t]
    valorproposto = rnorm( 1 , mean = 0 , sd = 0.1 )
    chanceaceitar = min( 1,
        targetRelProb( valoratual + valorproposto, meusdados ) /
        targetRelProb( valoratual , meusdados )
        )
    if ( runif(1) < chanceaceitar ) {
        cadeia[ t+1 ] <- valoratual + valorproposto
        } else {
            cadeia[ t+1 ] <- valoratual
	}
}

A gente salva numa variável o valor atual, na prática a gente gera uma quantidade, que é o valor proposto, que é uma quantidade para mudar o valor de theta atual, dai calculamos como comentamos acima a chance de aceitar esse novo valor, e se ele for aceito colocamos ele na cadeia, senão repetimos o último valor. E milagrosamente vamos visitar cada valor possível de theta um número de vezes proporcional ao valor de theta que gerou os dados. Na verdade não é milagre, essa estabilidade é prevista para cadeias de markov dado a chance de transição entre estados, mas pra mim é mágico como funciona 🙂

Então ao usarmos nosso algoritimozinho, o que vemos como resultado?

03

Uma boa estimativa de valor real de theta, o intervalo de 95% da distribuição cobre o valor real de theta.

Mas temos uma distribuição continua de thetas, não quero ficar olhando as coisas na forma de histograma, podemos olhar a mesma distribuição na forma de densidade contínua.

04

Lembrando que aqui a gente ja arrancou o início da cadeia, que como começamos num lugar qualquer, traz estimativas ruim. Vamos dar uma olhadinha na cadeia, para o início e um pedaço la na frente.

05

Como começamos num valor arbitrário qualquer, começamos num valor de theta que não tem nada a ver com o valor real e vamos caminhando ao valor próximo do real, ai chegando la ficamos rodeando em volta dele. Veja que o tamanho dos passos aqui depende do valor proposto, se propomos valores pouco diferentes do anterior, vamos andando devagarinho na cadeia, agora poderíamos propor passos mais largos, mais ai também, assim como poderíamos chegar num valor próximo do real mais rápido, poderíamos derrepente sair dele pra muito longe. Esse tipo de ajustamento, como o tamanho do passo, como propor uma esse passo que é a parte difícil, aqui tudo funciona mas em algum momento alguém deve ter ficado testando possibilidades e ficou frustado em como as coisas não funcionavam. Outra coisa é que pode ser proposto valores maiores que 1 ou negativos, menor que zero, veja que ficamos mudando esses valores para os extremos, 0 e 1, mas eu não sei qual o impacto de fazer isso se o valor de theta fosse 1 por exemplo, mas todas essas coisas influenciam no funcionamento e eficácia do uso de MCMC. Nisso pacotes como OpenBugs ou Stan salvam a vida, ao possibilitar a gente abstrair dessa coisas e focar em entender melhor nossas observações da natureza. Mas ainda sim é importante entender pelo menos um pouco do de tudo que estamos fazendo.

Bem acho que ja esta bom por aqui, depois olhamos mais sobre MCMC. Esse script como outros do blog estão sempre disponíveis no repositório do github aqui.

Referência:
John K. Kruschke 2011 Doing Bayesian Data Analysis: A Tutorial with R and BUGS. Academic Press / Elsevier

#Adaptado de uma função do livro Bayesian Data Analysis: A Tutorial with R and BUGS
#Semente para replicar os resultado do post.
set.seed(51)
 
#Gerando uma amostra de dados
probabilidade<-0.70
meusdados<-rbinom(n=14, size=1, prob=probabilidade)
 
#Definindo a função de likelihood binomial.
likelihood <- function( theta , data ) {
  z <- sum( data == 1 )
  N <- length( data )
  pDataGivenTheta <- theta^z * (1-theta)^(N-z)
  pDataGivenTheta[ theta > 1 | theta < 0 ] <- 0
  return( pDataGivenTheta )
}
 
plot(seq(0.01,0.99,0.01),likelihood(seq(0.01,0.99,0.01),meusdados),frame=F,
     xlab="Probabilidade",ylab="Likelihood",type="l")
 
#Função para definir o prior p(D),
prior <- function( theta ) {
  prior <- dbeta(theta,1,1)
  #exemplo de um prior bimodal
  #prior < dbeta( pmin(2*theta,2*(1-theta)),2,2 )
  prior[ theta > 1 | theta < 0 ] <- 0
  return( prior )
}
 
plot(seq(0.01,0.99,0.01),prior(seq(0.01,0.99,0.01)),frame=F,
     xlab="Probabilidade",ylab="Likelihood",type="l")
 
#Função auxiliar para ajudar a calcular a probabilidade relativa
targetRelProb <- function( theta , data ) {
  targetRelProb <-  likelihood( theta , data ) * prior( theta )
  return( targetRelProb )
}
 
# Tamanho da cadeia
Ncadeia = 100000
 
#Inicializando um vetor para salvar os resultados do MCMC
cadeia = rep( 0 , Ncadeia )
 
# Iniciando o MCMC com um valor arbritario
cadeia[1] = runif(1,min=0,max=1)
 
# Especificando um valor de burnning, ou seja o quanto do inicio da cadeia vamos jogar fora
# Aqui vamo jogar fora 10%
burnIn = ceiling( .1 * Ncadeia )
 
# Vamos contar quantas vezes aceitamos ou não os valores propostos
nAceito = 0
nRejeitado = 0
 
# Agora iniciando o "Random Walk"
for ( t in 1:(Ncadeia-1) ) {
    #Valor Atual
    valoratual = cadeia[t]
    #Propomos uma quantidade para alterar esse valor
    valorproposto = rnorm( 1 , mean = 0 , sd = 0.1 )
    #Computamos a chance de aceitar a mudança
    #Essa chance é igual a Likelihood*prior com o novo valor de theta
    #dividido pele Likeliood*prior antigo
    chanceaceitar = min( 1,
        targetRelProb( valoratual + valorproposto, meusdados ) /
        targetRelProb( valoratual , meusdados )
        )
    #Agora aqui a gente ve na sorte se aceitamos ou não o novo valor.
    #Veja que temos uma chance de fazer a mudança ou manter o ultimo valor
    #Se a gente sortei um número qualquer de 0 a 1 e compara se ele é
    #menor que a chance de aceitar, dessa forma mudamos numa probabilidade
    #da chance de mudar
    if ( runif(1) < chanceaceitar ) {
        #Se aceitarmos mudamos ele na cadeia
        cadeia[ t+1 ] <- valoratual + valorproposto
        #So vamos começar a contar depois do burnIn que escolhemos
        #Afinal, com o começo arbitrario, esperamos aceitar bastante
        #a não ser que escolhemos um valor proximo do final no começo
        if ( t > burnIn ) { nAceito <- nAceito + 1 }
 
        } else {
            # Quando não aceitamos o valor proposto mantemos o anterior
            cadeia[ t+1 ] <- valoratual
            # E incrementamos aqui apos o BurnIn, como acima
            if ( t > burnIn ) { nRejeitado <- nRejeitado + 1 }
	}
}
 
#Agora vamos olhar como ficou a cadeia final, ou seja, tirando os valores iniciais
cadeiafinal <- cadeia[ (burnIn+1) : length(cadeia) ]
 
#Alguns Gráficos
 
#-----------------------------------------------------------------------
HDIofMCMC = function( sampleVec , credMass=0.95 ) {
    sortedPts = sort( sampleVec )
    ciIdxInc = floor( credMass * length( sortedPts ) )
    nCIs = length( sortedPts ) - ciIdxInc
    ciWidth = rep( 0 , nCIs )
    for ( i in 1:nCIs ) {
        ciWidth[ i ] = sortedPts[ i + ciIdxInc ] - sortedPts[ i ]
    }
    HDImin = sortedPts[ which.min( ciWidth ) ]
    HDImax = sortedPts[ which.min( ciWidth ) + ciIdxInc ]
    HDIlim = c( HDImin , HDImax )
    return( HDIlim )
}
 
 
#------------------------------------------------------------------------
intconf<-HDIofMCMC(cadeiafinal)
 
hist(cadeiafinal,ylab="Frequência na cadeia final",xlab="Probabilidades",xlim=c(0,1))
lines(c(intconf[1],intconf[2]),c(500,500),lwd=4,lty=3,col="red")
abline(v=mean(cadeiafinal),lwd=4,col="red")
abline(v=probabilidade,lwd=2,col="blue")
 
legend("topleft",lty=c(3,1,1),lwd=c(4,4,2),bty="n",cex=0.8,col=c("red","red","blue"),
       legend=c("Intervalo de confiança","Média das estimativas","Probabilidade usada nas estimativas"))
 
plot(density(cadeiafinal),ylab="Frequência",xlab="Probabilidades",xlim=c(0,1),frame=F)
lines(c(intconf[1],intconf[2]),c(0.25,0.25),lwd=3,lty=2,col="red")
abline(v=mean(cadeiafinal),lwd=1,col="red")
abline(v=probabilidade,lwd=1,col="blue")
 
legend("topleft",lty=c(3,1,1),lwd=c(2,1,1),bty="n",cex=0.8,col=c("red","red","blue"),
       legend=c("Intervalo de confiança","Média das estimativas","Probabilidade usada nas estimativas"))
 
layout(matrix(c(1,2,3,4),ncol=2))
hist(cadeia[1:100],ylab="Frequência na cadeia",xlab="Probabilidades",main="100 primeiros valores",xlim=c(0,1))
plot(0,0,xlim=c(0,100),ylim=c(0,1),type="n",frame=F,main="Cadeia",
     xlab="Posição na cadeia",ylab="Probabilidade")
points(cadeia[1:100],type="l")
 
hist(cadeia[10000:10500],ylab="Frequência na cadeia",xlab="Probabilidades",main="Do 10000 ao 10500",xlim=c(0,1))
plot(0,0,xlim=c(0,500),ylim=c(0,1),type="n",frame=F,main="Cadeia",
     xlab="Posição na cadeia",ylab="Probabilidade",xaxt="n")
axis(1,at=seq(0,500,by=50),labels=seq(10000,10500,by=50))
points(cadeia[10000:10500],type="l")
 
# Evidencia para o modelo, p(D).
 
media<-mean(cadeiafinal)
desvio<-sd(cadeiafinal)
 
a =   media   * ( (media*(1-media)/desvio^2) - 1 )
b = (1-media) * ( (media*(1-media)/desvio^2) - 1 )
 
wtdEvid = dbeta( cadeiafinal , a , b ) /
    ( likelihood( cadeiafinal , meusdados ) * prior( cadeiafinal ) )
pData = 1 / mean( wtdEvid )
 
format(pData,scientific = F)

One thought on “Um exemplo do algoritimo de Metropolis–Hastings aplicado a distribuição binomial

Leave a Reply

Your email address will not be published. Required fields are marked *