Como diabos funciona um Gibbs Sampler?

Vamos dar uma olhada aqui num algorítimo de Gibbs Sampler. Bem o Openbugs, seu irmão winbugs e o Jags todos usam Gibbs Sampler para ajustar nossos modelos. O que esses programas fazem é fazer um sampler parecido com o que a gente vai ver aqui, so que ele monta todo o código sozinho, somente com o modelo que a gente manda pra ele e dados.
Definitivamente um grande corte de caminho, já que nos permite um bocado de abstração, e se focar no modelo que nos interessa, e em questões numéricas.

Mas nem por isso não vamos deixar de pelo menos entender um pouco o que está acontecendo, ter um pouco de intuição em como funciona o Gibbs Sampler. Apesar de ser um caso especial do Metropolis-Hasting, que a gente viu aqui, aqui e aqui, ele é um pouco mais complicado de entender. Enquanto no Metropolis-Hasting a gente tinha uma chance de aceitar cada novo valor para a distribuição posterior, ou repetir o anterior, aqui a gente sempre aceita o novo valor para a distribuição posterior.

Mas vamos olhar um exemplo funcionando. Primeiro vamos criar uns dados de exemplo. Vamos simular 1000 valores com média -2 e desvio padrão de 2.

Então ficamos com os seguintes dados:

01

Então vamos lembrar o que queremos fazer, queremos inferir a distribuição posterior, saber quais parâmetros realmente geraram os dados a partir dos dados e de um pressuposto que temos, a nossa distribui a priori, nosso prior.

Então queremos:

p(\mu,\sigma \mid Dados) \propto p(Dados \mid \mu,\sigma) \cdot p(\mu,\sigma)

O que são essas coisas?

Queremos saber distribuição posterior p(\mu,\sigma \mid dados), ou seja, quais são os valores de média e desvio padrão que geraram os nossos dados.

Sendo que nos temos os valores de média e desvio padrão que mais provavelmente gerarão os dados que estamos olhando p(Dados \mid \mu,\sigma) , o likelihood, máxima verossimilhança, que vimos até como conseguir ela aqui.

E temos também um presuposto, algo que achamos que deve ser, dado a literatura, os nossa intuição, o que temos p(\mu,\sigma).

Então é isso. Como aqui é um exemplo vamos usar um Prior bem ruim, chutando que a média deveria ser 5 e o desvio 3.162278. Veja que o desvio ta no prior como o tau.prior, e como a variância é 10, então colocar o tau.prior como 0.1 é a mesma coisa que 1/10, ou seja estamos colocando a variância como 10, e a raiz de 10 é 3 e alguma coisa.

# Valores do Prior  #
media.prior<-5
tau.prior<-0.1

Além disso setamos valores iniciais para começar o Gibbs Sampler.

# Valores iniciais #
tau<-1
mtau<-n/2

Criamos também vetores para receber os resultados

nsim<-5000
cadeia.media<-rep(0,nsim)
cadeia.variancia<-rep(0,nsim)

E temos todos os ingredientes, temos os dados, o nosso prior, valores iniciais para iniciar, e agora vamos ao Gibbs Sampler propriamente dito.

for (i in 1:nsim) {
 
  v<-1/(tau*n + tau.prior)
  m<-v*(tau*sum(dados)+tau.prior*media.prior)
  cadeia.media[i]<-rnorm(1,m,sqrt(v))
 
  tau<-rgamma(1,mtau,sum((dados-cadeia.media[i])^2)/2)
  cadeia.variancia[i]<-1/tau
 
}

E acabou, 5 linha de código, e se olharmos o resultado:

O valor estimado para média é de:

> mean(cadeia.media[1501:nsim]) [1] -1.999777

Com uma precisão:

> sd(cadeia.media[1501:nsim]) [1] 0.06566074

E no vaso de desvio padrão:

> sqrt(mean(cadeia.variancia[1501:nsim])) [1] 2.039191

Com uma precisão:

> sd(cadeia.variancia[1501:nsim]) [1] 0.1877103

Realmente funciona bem, já que os dados foram gerados com média -2 e desvio padrão 2, os distribuição posterior está bem próximo, englobando os valores reais.

Mas a pergunta é, o que aconteceu durante aquele loop que gerou esses resultados, o que o Gibbs Sampler fez?

Realmente é difícil entender aquele monte de multiplicação e somas, pelo menos pra mim é complicado. Mas vamos olhar um pouco melhor as variáveis, primeiro vamos ver como é a relação entre elas no grafo.

Vamos olhar ele e o script ao mesmo tempo.

01

v<-1/(tau*n + tau.prior)
m<-v*(tau*sum(dados)+tau.prior*media.prior)

No grafo, as variáveis são vértices, e as arestas aponta a direção em que uma variável influencia outra, se o vértice não recebe nenhuma aresta, ele não é influenciado por ninguém.

Veja que começamos gerando valores para v e m, para gerar esses valores, usamos os dados, os priores e o valor de tau. Com exceção de tau, todos são valores “estáticos”. Eles não mudam, veja que no grafo, só tem setinhas saindo deles a indo para os valores que eles influenciam, mas eles não recebem nenhuma setinha.

cadeia.media[i]<-rnorm(1,m,sqrt(v))

Agora esses dois valores são usados para gerar uma amostra da distribuição posterior, que vai ser guardado nesse vetor chamado cadeia.media.

tau<-rgamma(1,mtau,sum((dados-cadeia.media[i])^2)/2)

Depois usamos ele, mais um parâmetro externo chamado mtau, que depois vamos comentar, mas que serve para regularmos os tamanhos dos passos que vamos dar na caminhada aleatória do Gibbs Sampler, mas foco na hora que usamos a distribuição gamma para gerar um novo valor de tau, usamos o valor recém adicionado a distribuição posterior, usamos ele para calcular os mínimos quadrados, que é o que? Uma função de custo, se o nosso valor novo é “bom”, perto do real, a soma aqui sera algo mais próximo de zero, quanto pior o novo valor adicionado, mais sera o valor aqui, mas por enquanto vamos continuar olhando o grafo e se ater ao fato que ele é usado para estabelecer o novo valor de tau.

cadeia.variancia[i]<-1/tau

Usamos tau também para estabelecer o novo valor de variância, nosso outro parâmetro da distribuição normal. E assim acabou uma iteração, e agora começamos de novo.

Mas o detalhe é olhar esse ciclo que é formado aqui.

01

Um valor da distribuição posterior é usado para gerar o próximo valor de tau, que é usado para gerar os próximos valores de v e m que são usados para gerar o próximo valor da distribuição posterior, que é usado para gerar o próximo valor de tau, e assim por quantas iterações quisermos, para o tamanho da cadeia de mcmc que definimos.

Agora pensando um pouco mais nos valores, vamos olhar os primeiros 50 valores de tau, v e m.

01

Veja que quando o valor da média para o próximo valor que vamos sortear sobe, a valor de tau desce, quando tau desce o valor de m sobe, e assim vai.
Se a gente olhar com mais carinho essa linha

tau<-rgamma(1,mtau,sum((dados-cadeia.media[i])^2)/2)

Pense que se estamos com um valor muito bom para o parâmetro média, essa soma dará um valor muito próximo de zero, e a implicação disso ao se gerar um número aleatório para tau é o seguinte:

01

Ou seja, próximo de zero, deixamos tau ter a possibilidade de ser grande, já se o valor é grande, pegamos um valor muito ruim para a média, deixamos o valor sempre pequeno, próximo de zero.

Agora ao olhar a linha que gera o valor que será usado para propor a próxima amostra para a distribuição posterior da média a gente entende isso.

m<-v*(tau*sum(dados)+tau.prior*media.prior)

Veja que se, ao extremo, tau for zero, ele tira a influencia dos dados, do likelihood e deixa os valores de prior comandando o que deve ser a próxima amostra da média. O contrario, se tau é um valor grande, ele aumenta a influencia do likelihood no que será o a próxima amostra da média. E lembre-se que o valor do posterior é algo entre a máxima verossimilhança e o nosso prior. E assim estamos caminhando, entre esses valores.
Outra coisa legal a se notar é que temos uma multiplicação pelo n, e na linha anterior

v<-1/(tau*n + tau.prior)

Veja que o tau é multiplicado pelo n, e n é o tamanho da amostra, e daqui sai a nossa influencia da quantidade de dados, no que será a distribuição posterior, lembra quando comentamos que quanto mais dados, menor a importância relativa do prior, em algum post que não lembro mais hehe?

Mas é isso ai, assim funciona o Gibbs sampler, sempre caminhando sem parar, diferente do Metropolis-Hastings.

Podemos olhar que o resultado e ver como a convergência foi rápida

01

Não é a toa, usamos um valor gigantesco de amostra, 1000 amostras.
E veja como os valores de média e desvio que conseguimos como a distribuição posterior ficam próximos ao valor real, mais precisamente, o intervalo de confiança de 95% engloba o valor real de média e desvio padrão usado para produzir os dados.

01

E esse é o Gibbs Sampler, que tanto falamos. E o que o OpenBugs faz quando a gente chama ele? Ele constrói um algorítimo desse tipo, mas com o modelo que definimos no model. Mas a idéia geral é por aqui.

É um bocado mais complicado que o Metropolis-Hasting, principalmente até captar o que ele está fazendo com os dados. Eu precisei ficar rodando ele com menos dados, e mudando vários dos parâmetros para conseguir entender mais ou menos o que está acontecendo e espero não ter entendido nada errado, por favor me corrijam se eu estou falando alguma abobrinha.

Outra coisa, se olharem no script original, do Brian Neelon, verão que ele tem 2 parâmetros, a e b, que deixa em zero. Eu suprimi esses parâmetros aqui, mas veja que eles influenciam nos dois números que servem de parâmetro para a função rgamma que faz o novo valor de tau a cada iteração. Isso possibilita o controle do tamanhos dos passos que damos, e vai influenciar na distribuição posterior. Esse tipo de coisa pode ser otimizada para valores que dão soluções melhores, mas por enquanto, como mal entendemos o Gibbs Sampler, não vamos mexer neles. Mas vendo o algorítimo eu acho que ajuda a desmistificar um pouco e nos deixar mais confiantes com a estatística bayesiana, ja que não tem nada muito milagroso aqui.

Bem vamos ficar por aqui, que esse post ja está meio grande, mas o script está no repositório recologia e no topo dele existe um link para o script original da onde eu copie o Gibbs Sampler mostrado aqui.

# Original:  http://people.duke.edu/~neelo003/r/Normal%20Mean%20Gibbs.r
# Gibbs sampler para uma distribuição normal
###################################
set.seed(250)
 
########
# Data #
########
n<-1000
media.real<--2
desvio.padrao.real<-2
dados<-rnorm(n,media.real,desvio.padrao.real)
 
hist(dados,main="Histograma",ylab="Frequência",xlab="Valores")
abline(v=mean(dados),lty=2,col="red",lwd=2)
 
#####################
# Valores do Prior  #
#####################
media.prior<-5
tau.prior<-0.1
 
####################
# Valores iniciais #
####################
tau<-1
mtau<-n/2
 
####################################
# Vetor para guardar os resultados #
####################################
 
nsim<-5000
cadeia.media<-rep(0,nsim)
cadeia.variancia<-rep(0,nsim)
 
#################
# Gibbs Sampler #
#################
vc<-rep(0,nsim)
mc<-rep(0,nsim)
tauc<-rep(0,nsim)
 
 
for (i in 1:nsim) {
 
  v<-1/(tau*n + tau.prior)
  m<-v*(tau*sum(dados)+tau.prior*media.prior)
  cadeia.media[i]<-rnorm(1,m,sqrt(v))
 
  tau<-rgamma(1,mtau,sum((dados-cadeia.media[i])^2)/2)
  cadeia.variancia[i]<-1/tau
 
}
 
library(igraph)
grafo<-graph.formula(n-+v,tau-+v,tau.prior-+v,v-+m,dados-+m,tau.prior-+m,media.prior-+m,m-+cadeia.media,v-+cadeia.media,
                   cadeia.media-+tau,dados-+tau,tau-+cadeia.variancia,tau-+m)
tkplot(grafo,vertex.size=60,vertex.color="gray",vertex.label.color="black",edge.arrow.size=2,edge.color="black",
       vertex.label.cex=1.2)
 
grafo<-graph.formula(m-+cadeia.media,v-+cadeia.media,cadeia.media-+tau,tau-+cadeia.variancia,tau-+m,tau-+v)
tkplot(grafo,vertex.size=60,vertex.color="gray",vertex.label.color="black",edge.arrow.size=2,edge.color="black",
       vertex.label.cex=1.2)
 
layout(matrix(c(1,2,3,4),nrow=2,ncol=2,byrow=T))
curve(dgamma(x,shape=1,rate=0.01),0,10,main="Rate = 0.01",frame=F)
curve(dgamma(x,shape=1,rate=0.1),0,10,main="Rate = 0.1",frame=F)
curve(dgamma(x,shape=1,rate=1),0,10,main="Rate = 1",frame=F)
curve(dgamma(x,shape=1,rate=10),0,10,main="Rate = 10",frame=F)
 
##########################
# Distribuição Posterior #
##########################
mean(cadeia.media[1501:nsim])
sd(cadeia.media[1501:nsim])
 
sqrt(mean(cadeia.variancia[1501:nsim]))
sd(cadeia.variancia[1501:nsim])
 
 
#Evolução da cadeia
plot(1:nsim,cadeia.media,type="l",col="red",frame=F)
abline(h=mean(cadeia.media),lwd=2,lty=2)
abline(v=1500,lwd=4,lty=3,col="gray")
text(750,-1.76,"Burn-in")
 
#Distribuição posterior
plot(cadeia.media[1:500],sqrt(cadeia.variancia[1:500]),type="b",frame=F,xlab="Média",ylab="Desvio Padrão",cex=0.5,pch=19)
library(MASS)
contornos<-kde2d(cadeia.media[1:500],sqrt(cadeia.variancia[1:500]))
contour(contornos,add=T,lwd=4,nlevels = 6,col=rev(heat.colors(6)))

Leave a Reply

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