O algoritimo de Metropolis–Hastings aplicado a distribuição normal

Bem a gente ja viu aqui como é um algorítimo de Metropolis–Hastings usado para a distribuição binomial. Mas para muitas coisas a gente usa a distribuição normal como ele fica para a distribuição normal?

Eu peguei aqui um exemplo da página do Brian Neelon aqui. Tem bastante coisa legal para olhar, mas vamos la.

A primeira parte interessante, é que como já vimos nos modelos usando openbugs que fazemos, é que o desvio-padrão, a gente estima indiretamente, usando o tau que é 1/variância. Aqui tem uma explicação melhor e mais detalhada, mas é simples pensar que a média a gente pode usar uma distribuição normal para representar, agora o desvio-padrão não, já que ele não pode assumir valores negativos, assim como uma distribuição de poisson, mas pode assumir qualquer valor entre 0 e + infinito, dai precisamos de algo como a distribuição gamma que é o conjugate prior do desvio padrão da distribuição normal, nesse post aqui tem algo sobre conjugate prior se alguém quiser olhar.

Mas então a gente ta fazendo isso aqui.

normal

Beleza. A idéia é que o que falei ali em cima é o que esta nessa figura, eu espero, se não estiver entendendo nada errado.

Então precisamos do nosso Prior

###########
# Priors  #
###########
mu.prior<-0
sigma2.prior<-10
tau.prior<-1/sigma2.prior
a<-0.001
mtau<-a+n/2

Além disso começamos a cadeia em qualquer lugar, eu comecei aqui em média igual a 0.1 e tau igual a 0.1 também, só para ver aquela caminhadinha inicial, mas normalmente queremos começar com um valor aleatório aqui. Ja que se de qualquer lugar que começamos, terminamos no mesmo lugar, é muito mais convincente que se delimitamos sempre os lugares onde começamos.

####################
# Inicio da cadeia #
####################
 
mu<-0.1
tau<-0.1
nsim<-5000
 
##########
# Cadeia #
##########
cadeia.mu<-c(mu,rep(0,nsim-1))
cadeia.tau<-c(tau,rep(0,nsim-1))
A<-0

Agora é so rodar o MCMC. Mas antes um pequeno comentário. Os valores de máxima verosimilhança ficam bem pequenos, apesar de que para a distribuição binomial ficou algo com 4 casas depois do zero, aqui o likelihood vai ser um produtório de vários números pequenininhos, que vai resultar num número incrivelmente pequeno.

Então gerando vários valores possíveis para média e calculando o likelihood vemos o seguinte.

01

Agora olhe os valores no eixo y, 3.502837e-75 no topo la, é um número muito muito pequeno. Números grandões gastam assim (grande no sentido de número de casas) gastam muito bytes de memória para serem guardados, inclusive temos um limite que podemos usar, quando tentamos usar mais que o limite temos um estouro da memória ou overflow.

Agora o que a gente vê normalmente é as pessoas trabalhando com loglikelihood, que é meramente calcular o log natural do valor usado ali em cima. E quando a gente faz isso, veja o que acontece com o mesmo gráfico.

02

Os valores são bem menores, já que a gente muda do produtório para o um somatório. Agora séra que isso da problema na estimativa de maximum likelihood? Que nada mais é que o pico desses gráficos ae!

Veja só, a gente fez um vetor chamado valores, e cada elemento desse vetor teve o likelihood calculado, o maior valor que encontramos foi:

> max(likelihood) [1] 3.502837e-75

Agora a gente pode olhar qual elemento do vetor é isso.

> which(likelihood==max(likelihood)) [1] 101

Agora se a gente olhar qual valor representa esse elemento do vetor de valores vemos que:

> valores[101] [1] 5

E olha a média dos dados são, a média é exatamente a estimativa de maxima verosimilhança;

> mean(dados) [1] 5.038861

E se olhamos o proximo, o elemento 102, o elemento 101 é muito mais próximo da média.

> valores[102] [1] 5.1

Agora o valor em log é um valor diferente, ja que fizemos o log do valor acima:

> max(loglikelihood) [1] -171.4403

Mas o pico continua sendo no mesmo lugar:

> which(loglikelihood==max(loglikelihood)) [1] 101

Agora algo legal é como, devido aos números serem pequenos, não da para calcular exatamente, ou seja, o contrario de log, que é o exponencial do loglikelihood não é igual o likelihood.

> exp(loglikelihood)[100]==likelihood[100] [1] FALSE

Mas isso é por causa de aproximação, não tem como representar no computador um número como pi por exemplo, ja que não da para ter infinitas casas, mas para fazer comparações assim, sem sofrer com erros de aproximação, comum quando entramos nesse tipo de conta, da pra usar a função all.equal do R.

> all.equal(exp(loglikelihood)[100],likelihood[100]) [1] TRUE

Mais comentários sobre isso estão no faq do R.

Mas vamos a como fazer o MCMC usando Metropolis-Hastings para uma distribuição normal.

############
#   MCMC   #
############
for (t in 2:nsim){
    mus<-mu+0.3*rt(1,3)
 
    r<- (sum(dnorm(dados,mus,1/sqrt(tau),log=T) - dnorm(dados,mu,1/sqrt(tau),log=T)) +
         dnorm(mus,mu.prior,sqrt(sigma2.prior),log=T) -  dnorm(mu,mu.prior,sqrt(sigma2.prior),log=T))
 
    if (log(runif(1))<r) {
        mu<-mus
        A<-A+1
        }
 
    cadeia.mu[t]<-mu
 
    tau<-rgamma(1,mtau,b+sum((dados-mu)^2)/2)
    cadeia.tau[t]<-tau
}

Bem curtinho não. Mas seguimos o mesmo esquema, a cada iteração, sugerimos uma nova média para a cadeia, aqui usando uma distribuição t, depois disso calculamos uma chance de aceitar, que aqui parece um pouco mais complicado, mas note que estamos usando o log em tudo, log=T, então os produtos viram soma, divisão subtração, depois vemos isso com calma, mas é a mesma coisa feita para a distribuição binomial, o likelihood * prior dividido pelo do valor anterior da cadeia, mas a gente ve a chance de aceitar o novo valor ou não, dai sorteamos se vamos aceitar ele ou não, se sim adicionamos ele a cadeia, se não repetimos o ultimo valor, depois disso usamos a distribuição gamma para ajustar a precisão tau com a nova ou antiga média. Essa parte precisa de um pouco mais de estudo mas vamos indo.

Na verdade ja acabou, temos nossa cadeia final, ou cadeia, uma para média e uma para o desvio padrão. Podemos olhar os nossos resultados.

03

E as cadeia convergem bonitinho, em baixo nos histogramas eu tentei olhar varios pedacinhos da cadeia para a média e o desvio padrão, e vemos como ambas vão para intervalos que cobrem os valores reais usados para gerar os dados.

Podemos também olhar como a cadeia caminha em relação aos dois parametros, média e desvio padrão.

04

Novamente vemos que existem muitas nuances, que espero fazer um post experimentando o que acontece de acordo com a estrategia usada para gerar novos valores de parâmetros propostos, mas ja da para começar a ter uma ideia de como é um MCMC para uma distribuição normal que tem dois parâmetros, diferente do que fizemos para a distruição binomial que so tinha um parâmetro.

Bem é isso, o repositório recologia no github ta atualizado se alguém quiser alguma coisa, e ficamos por aqui. Quem quiser pode olhar o script original desse MCMC na página do Brian Neelon.

#Original
#http://people.duke.edu/~neelo003/r/
set.seed(250)
 
#######
#Dados#
#######
 
n<-100
mu.real<-5
sigma2.real<-2
tau.real<-1/sigma2.real
dados<-rnorm(n,mu.real,sqrt(sigma2.real))
 
a<-0.001
b<-0.001
 
###########
# Priors  #
###########
mu.prior<-0
sigma2.prior<-10
tau.prior<-1/sigma2.prior
mtau<-a+n/2
 
####################
# Inicio da cadeia #
####################
 
mu<-0.1
tau<-0.1
nsim<-5000
 
##########
# Cadeia #
##########
cadeia.mu<-c(mu,rep(0,nsim-1))
cadeia.tau<-c(tau,rep(0,nsim-1))
A<-0
 
#Valores para avaliar a maxima verossimilhança
valores<-seq(-5,15,by=0.1)
likelihood<-vector()
 
for(i in 1:length(valores)) {
    likelihood[i]<-prod(dnorm(dados,valores[i],sqrt(2),log=F))
    }
 
#likelihood
plot(likelihood~valores,type="b",frame=F,ylab="likelihood")
 
loglikelihood<-vector()
 
for(i in 1:length(valores)) {
    loglikelihood[i]<-sum(dnorm(dados,valores[i],sqrt(2),log=T))
    }
 
#Loglikelihood
plot(loglikelihood~valores,type="b",frame=F,ylab="loglikelihood")
 
exp(loglikelihood)[100]==likelihood[100]
all.equal(exp(loglikelihood)[100],likelihood[100])
 
max(likelihood)
 
which(likelihood==max(likelihood))
valores[101]
mean(dados)
valores[102]
 
max(loglikelihood)
which(loglikelihood==max(loglikelihood))
 
 
############
#   MCMC   #
############
for (t in 2:nsim){
    mus<-mu+0.3*rt(1,3)
 
    r<- (sum(dnorm(dados,mus,1/sqrt(tau),log=T) - dnorm(dados,mu,1/sqrt(tau),log=T)) +
         dnorm(mus,mu.prior,sqrt(sigma2.prior),log=T) -  dnorm(mu,mu.prior,sqrt(sigma2.prior),log=T))
 
    if (log(runif(1))<r) {
        mu<-mus
        A<-A+1
        }
 
    cadeia.mu[t]<-mu
 
    tau<-rgamma(1,mtau,b+sum((dados-mu)^2)/2)
    cadeia.tau[t]<-tau
}
 
###########
# Gráfico #
###########
 
# Cadeias
layout(matrix(c(1,1,1,2,2,2,1,1,1,2,2,2,3,4,5,6,7,8),nrow=3,byrow=T))
 
plot(seq(1,nsim) ,cadeia.mu[1:nsim],type="l",xlab="Iteration", ylab="mu",frame=F,
main="Média",col=2,lty=1,ylim=c(0,10))
abline(h=mu.real)
plot(seq(1,nsim) ,cadeia.tau[1:nsim],type="l",xlab="Iteration", ylab="mu",frame=F,
main="Tau",col=2,lty=1,ylim=c(0,1))
abline(h=1/sigma2.real)
 
par(mar=c(2,2,2,2),font=3,cex=0.7)
 
hist(cadeia.mu[1:100],xlab="",ylab="Frequência",main="",cex=0.7,xlim=c(3,6),breaks=seq(0,6,by=0.02))
hist(cadeia.mu[2000:2500],xlab="",ylab="Frequência",main="",cex=0.7,xlim=c(3,6),breaks=seq(3,6,by=0.02))
hist(cadeia.mu[4500:5000],xlab="",ylab="Frequência",main="",cex=0.7,xlim=c(3,6),breaks=seq(3,6,by=0.02))
hist(cadeia.tau[1:500],xlab="",ylab="Frequência",main="",cex=0.7,xlim=c(0,1),breaks=seq(0,1,by=0.01))
hist(cadeia.tau[2000:2500],xlab="",ylab="Frequência",main="",cex=0.7,xlim=c(0,1),breaks=seq(0,1,by=0.01))
hist(cadeia.tau[4500:5000],xlab="",ylab="Frequência",main="",cex=0.7,xlim=c(0,1),breaks=seq(0,1,by=0.01))
 
#Evolução da cadeia
plot(cadeia.mu[1:500],sqrt(1/cadeia.tau[1:500]),type="b",frame=F,xlab="Média",ylab="Desvio Padrão",cex=0.5,pch=19,
     xlim=c(0,6),ylim=c(0,6))
library(MASS)
contornos<-kde2d(cadeia.mu[1:500],sqrt(1/cadeia.tau[1:500]))
contour(contornos,add=T,lwd=2,nlevels = 6,col=rev(heat.colors(6)))

Leave a Reply

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