graphics.off()
windows(record=TRUE)

n = 12  # mulh 20-24 em Bodó
Y = 3   # num. de filhos nasc. no ultimo ano

# distrib. a priori (vem de uma análise da distrib. das taxas mais
#  confiáveis dos munic. com >200 mulheres na faixa etária 20-24
mu.priori    = -2.4
sigma.priori = 0.32
tau.priori   = 1/(sigma.priori^2)

# a priori normal + veross. Poisson
logP = function(log.theta, n, Y) {
  theta = exp(log.theta)
  -tau.priori/2 * (log.theta - mu.priori)^2 - n * theta + Y*log(theta)
}

log.theta = rep(NA, 10000)

log.theta[1]   = -1
logP.atual = logP(log.theta[1], n, Y)

for (i in 2:10000) {
  log.theta.propo = log.theta[i-1] + rnorm(1,0,.25) 
  logP.propo  = logP(log.theta.propo, n, Y)
  logu        = log(runif(1,0,1))
  if (logu < logP.propo - logP.atual) {
    log.theta[i]   = log.theta.propo
    logP.atual = logP.propo 
  }
  else {log.theta[i] = log.theta[i-1]}
}

# taxa de aceitação 
mean(diff(log.theta) != 0)


# plot(  sort(theta), seq(theta)/length(theta), type='l', lwd=3, col='blue')
# lines( sort(theta), pnorm( log(sort(theta)), -2.3, 0.4), col='grey', lwd=3)

plot( log.theta, type='l', main='Realizações de log.theta')

theta = exp(log.theta)
plot(theta, type='l', main='Realizações de theta')

val.theta = seq(0,.40,.001)

# density of log theta
plot(density(log.theta),main='Dist. a priori e a posteriori',
         xlab='log.theta', ylab='dens', col='blue',lwd=4)
lines(log(val.theta), dnorm(log(val.theta) ,mu.priori, sigma.priori) ,lwd=2, col='grey')

# density of theta
de = density(log.theta)
plot(exp(de$x), de$y, type='l',xlim=c(0,.40), col='blue',lwd=4,
     main='Dist. a priori e a posteriori', xlab='theta',ylab='dens',
     bty='l', ylim=c(0,1.5))
abline(v=seq(0,.40,.05),lty=2,col='grey')

lines(val.theta, dnorm(log(val.theta) ,mu.priori, sigma.priori) ,lwd=2, col='grey')