Olhando o MCMC mais de perto…

Algumas vezes pode ser frustrante tentar usar estatística bayesiana.
Mas vamos continuar usando o Metropolis-Hasting para um modelo de regressão linear. E para tanto simulamos mais alguns dados

01

Primeiro vamos ver o seguinte, eu alterei o prior para sempre usar uma distribuição normal agora

prior <- function(param){
alfa = param[1]
beta = param[2]
sd = param[3]
alfaprior = dnorm(alfa, sd=2, log = T)
betaprior = dnorm(beta, sd=1, log = T)
sdprior = dnorm(sd, sd=5, log = T)
return(alfaprior+betaprior+sdprior)
}

E alterei um pouquinho a função que propõe um novo passo.

proposalfunction <- function(param){
return(rnorm(3,mean = param, sd= c(1,0.1,1)))
}

Aqui reside a seguinte decisão, se você quer passos grandes ou passos pequenos. Como assim? A cada novo passo, uma condição é testada, se vamos ficar no mesmo lugar ou vamos para esse novo lugar, se esse novo lugar tem um likelehood ruim, ficamos onde estamos, lembre-se que tem um if que faz essa restrição, sempre é sorteado uma chance de aceitar ou não o novo passo, mas quanto melhor posicionado estamos perto daquilo que deve ser a verdade, menor a chance de aceitar um passo para muito longo.
O ponto também são que usamos alfa e beta, os parâmetros do modelo, o anterior da cadeia, então ali nos desvios, quando estamos falando sd=número, é o quão longe podemos sair da onde estamos.

metropolis-hastings

Mas pagaremos um preço se dermos passos curtos, que é de ter que caminhar mais, meio óbvio não.

Mas vamos olhar num plano, como mudam os valores de intercepto e inclinação ao longo da cadeia. Então o modelo sempre tem um ponto que ele esta , e vamos unindo os pontos com linhas, então a linha vermelha é a união de como esta indo a cadeia. e vamos apagar os pontos mais antigos para facilitar a visualização, clareando eles com o tempo. Como é uma simulação, vamos deixar um pontinho azul que representa a “verdade”, que infelizmente nunca saberemos numa análise com dados reais, mas como aqui o intuito é entender melhor como a análise ta se comportando, temos esse privilegio de saber a verdade.

E vemos o seguinte

rw

Se alguém lembrou da figurinha do random walk não é por acaso, só que aqui a probabilidade da uma direção, e o random walk é em torno dos parâmetros reais.

Olhe como começamos de um ponto bem longe e a cadeia caminha para perto do ponto azul, que são os parâmetros reais usados na simulação, e la começa a ficar rodando. Usamos uma cadeia de 50 mil passos aqui, e fui colocando exponencialmente os passos no gráfico, cada vez mais e mais passos de uma vez.
Outra coisa é que não importa da onde começamos, a cadeia sempre vai convergir para a “verdade”. Por isso vai ser comum ao invés de usar um ponto inicial fico, começar de um ponto aleatório, já que se o ponto inicial influenciar no resultado final, algo não deve estar certo.

Mas vamos olhar como esta a distribuição desses 50 mil pontos ai nesse plano dos dois parâmetros (lembrando que a cadeia tem 3 partes, existe o desvio que não esta ai representado).

03

Fazendo curvas de nível para representar as densidades, igual curvas de altitude em mapinhas, vemos que o pico esta bem próximo da verdade.

Outra coisa que a gente pode ver é como parece que um parâmetro compensa o outro.
Veja a cadeias em relação aos parâmetros reais.

02

Veja como os valores de inclinação ficam um pouco abaixo da linha vermelha (o parâmetro real) e o intercepto faz o contrario, fica acima, mas quando um muda outro muda junto, vemos isso na figura anterior também, já que a figura parece uma bolinha achatada, mas que tem uma inclinação.

Por isso é importante olhar como ficam as cadeias no final de uma análise e não somente computar alguma estatística do tipo aceitação (acceptance), que apesar de ficar alta (alto ou baixo sempre são conceitos abstratos sem algo para comparar), quase 55%, ainda podemos ver esse pequeno desviozinho.

Realmente todo esse processo é bem difícil de fazer, e ficar mudando os algoritmos para qualquer alteraçãozinha, mas agora nos sabemos (apesar que bugs usa um esquema ligeiramente diferente, o gibbs sampler) o que o BUGS vai fazer com os dados, ele pula toda essa parte sofrida de fazer o algoritmos, e só vamos falar o likelehood e o prior e dar como queremos as cadeias (aqui só estávamos fazendo uma, mas o ideal é fazer mais de uma) e ele faz essa parte chata, e nos preocupamos em interpretar os resultados.

Bem por último, para aqueles que tem dificuldade de entender as curvas de nível, podemos montar aquele gráfico em perspectiva, onde a montanha é o empilhamento dos pontos e quanto mais alto a montanha, maior era a densidade de pontos ali.
Gráficos assim são ruins de visualizar, principalmente representar numeração em 3d é, mas fica tão legal ele girando que não quis jogar ele fora.

per

Agora podemos voltar a usar o Openbugs ou Winbugs com um pouco mais de conforto quanto ao que diabos esta acontecendo. Ou mesmo o stan 🙂

Aqui também estão o blog original onde vi a ideia de como fazer no R o Metropolis-hasting, o Theoretical Ecology, e se alguém quer ver um gibbs sampler, olhem o Darren Wilkinson’s research blog. É isso 🙂

###################################
# Metropolis-Hastings MCMC Denovo #
###################################
#dados de exemplo
set.seed(5)
 
alfa.verdadeiro <- 4
beta.verdadeiro <- 8
desv.verdadeiro <- 6
tamanho.amostra <- 30
 
# Criando valores do preditor
x <-runif(30,0,10)
# Criando valores dependentes a partir de "beta*x + alfa + N(0,sd)"
y <-alfa.verdadeiro + beta.verdadeiro*x + rnorm(
n=tamanho.amostra,mean=0,sd=desv.verdadeiro)
 
#Grafico
plot(x,y, main="Dados Simulados",xlab="Preditor",ylab="Resposta",
frame.plot=F,pch=19)
abline(lm(y~x))
 
#Likelihood
likelihood <- function(param){
alfa = param[1]
beta = param[2]
sd = param[3]
pred = alfa + beta*x
singlelikelihoods = dnorm(y, mean = pred, sd = sd, log = T)
sumll = sum(singlelikelihoods)
return(sumll)
}
 
#Prior distribution (Agora com distribuição normal para todos os priors)
prior <- function(param){
alfa = param[1]
beta = param[2]
sd = param[3]
alfaprior = dnorm(alfa, sd=2, log = T)
betaprior = dnorm(beta, sd=1, log = T)
sdprior = dnorm(sd, sd=5, log = T)
return(alfaprior+betaprior+sdprior)
}
 
posterior <- function(param){
return (likelihood(param) + prior(param))
}
 
######## Metropolis algorithm ################
#Note que aqui voce vai configurar o tamanho dos passos.
proposalfunction <- function(param){
return(rnorm(3,mean = param, sd= c(1,0.1,1)))
}
 
run_metropolis_MCMC <- function(startvalue, iterations){
chain = array(dim = c(iterations+1,3))
chain[1,] = startvalue
for (i in 1:iterations){
proposal = proposalfunction(chain[i,])
probab = exp(posterior(proposal) – posterior(chain[i,]))
if (runif(1) < probab){
chain[i+1,] = proposal
}else{
chain[i+1,] = chain[i,]
}
}
return(chain)
}
 
#vamos começar de algum lugar aleatorio agora
#mas cuidado, um start value ruim, pode não ser bom
startvalue <- runif(3,0,10)
chain <- run_metropolis_MCMC(startvalue, 50000)
burnIn <- 25000
 
#Agora a cadeia é maior, mas temos um acceptance menor
#Mas lembre-se que o acceptance vai depender do prior
#e do proposal
acceptance <- 1-mean(duplicated(chain[-(1:burnIn),]))
acceptance
 
#olhando estabilidade das cadeias
par(mfrow = c(3,1))
plot(chain[-(1:burnIn),1], type = "l",
xlab="" , main = "Cadeia dos valores de intercepto", )
abline(h=alfa.verdadeiro,lwd=2,col="red")
plot(chain[-(1:burnIn),2], type = "l",
xlab="" , main = "Cadeia dos valores da inclinação", )
abline(h=beta.verdadeiro,lwd=2,col="red")
plot(chain[-(1:burnIn),3], type = "l",
xlab="" , main = "Cadeia dos valores do desvio", )
abline(h=desv.verdadeiro,lwd=1,col="red")
 
#Grafico animado do MCMC para os parametros
#intercepto e inclinação
passos<-c(round(seq(1,223,length=200)^2),50000)
 
for(i in 1:200) {
jpeg(sprintf("RW%03d.jpg",i), width = 300, height = 300)
plot(0,0,type="n",xlab="Intercepto",ylab="Inclinação",
xlim=c(-2,12),ylim=c(0,10),main=paste("Passo",passos[i+1]))
cores<-rev(heat.colors(i, alpha =0.5))
for(j in 1:i) {
lines(chain[passos[j]:passos[j+1],1:2],type="l",col=cores[j])
}
points(alfa.verdadeiro,beta.verdadeiro,pch=19,col="blue")
dev.off()
}
 
system("convert RW*.jpg -delay 5 rw.gif")
system("rm RW*.jpg")
 
#preparando a densidade dos pontos
library(MASS)
z<-kde2d(chain[,1],chain[,2])
 
#contorno da densidade de pontos (passos)
plot(chain[,1],chain[,2],col=heat.colors(nrow(chain),alpha=0.4))
contour(z,lwd=3,add=T)
points(alfa.verdadeiro,beta.verdadeiro,pch=19,col="blue")
 
#perspectiva animada da densidade de pontos
theta<-seq(0,360,2)
for(i in 1:length(theta)) {
jpeg(sprintf("per%03d.jpg",i), width = 350, height = 350)
persp(z,phi = 15,theta=theta[i],xlab="Intercepto",
ylab="Inclinação",zlab="Densidade")
dev.off()
}
system("convert per*.jpg -delay 30 per.gif")
system("rm per*.jpg")
 
#lembre-se que trabalhar em outro diretorio apenas para os graficos
#pode evitar transtornos.

Leave a Reply

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