########################################################
# Criar distribuições a priori para a taxa de fecund.
# em todas as 7 faixas etárias 15-19...45-49
# a partir de dados que vêm dos municípios com 
# pelo menos 500 mulheres na faixa etária
#
# Usa o Teorema de Bayes para estimar o perfil da fecund. 
#  para um município pequeno (Bodó RN) e para um
#  município grande (Natal RN)
########################################################

graphics.off()   # desliga o painel gráfico
windows(record=TRUE)  # abre uma nova janela, fora do R studio, para os gráficos

# local dos dados na internet
url = 'http://ufrn.schmert.net/dados/fecund.2010.csv'

# baixando os dados para um data.frame
D = read.csv( url)

#################################################################

mu.priori    = rep(0,7)
sigma.priori = rep(0,7)

names(mu.priori)    = seq(15,45,5)
names(sigma.priori) = seq(15,45,5)

for (x in seq(15,45,5)) {

  # quais são os nomes dos variáveis que contêm nascimentos e mulheres para esta faixa etária? 
  nome.var.mulh = paste0("Mu",x)
  nome.var.nasc = paste0("Cr",x)
  
  # cria vetores mulh e nasc com as informações para esta faixa etária
  #  usa [[  ]] para indexar elementos de D
  mulh = D[[nome.var.mulh]]
  nasc = D[[nome.var.nasc]]
  
  # vetor de variáveis lógicas para (n > 500)
  grande = (mulh > 500)
  
  print( table(D$Região, grande))
  
  # calcular a estimativa nascimentos no último ano / mulher
  #   para todos os munic. do Brasil
  tx = nasc / mulh
  
  # distribuicao de estimativas para os *grandes* municipios 
  plot( density( tx[grande]), lwd=3, col='red', main=paste('Faixa Etária',x))
  
  # acrescenta ao plot a distibuicão para todos os municípios 
  #  nota: para acrescentar, usa-se o comando 'lines'
  lines( density( tx), lwd=1, col='black')
  
  ##
  plot( density( log(tx[grande]), na.rm=TRUE), lwd=6, col='red',main=paste('Faixa Etária',x))
  
  zero = (tx ==0) # quais estimativas =0?
  
  m = mean(log(tx[grande & !zero]))
  s = sd(log(tx[grande & !zero]))

  mu.priori[paste(x)]    = m  
  sigma.priori[paste(x)] = s  
  
  val = seq(-9,-1, 0.1)
  lines( val, dnorm(val, m, s), col='blue', lwd=3)
} # for x

####################################################################
# estimativas de max. veross. (MV) e max. a posteriori (MAP)
# para um municip. pequeno (Bodó RN) e um munic. grande (Natal RN)
####################################################################

## log. da distrib. a posteriori, se a distrib. a priori for normal na
##  escala logorítmica e a veross. for Poisson

logP = function(theta,n,Y,mu.pri, sigma.pri) {
  tau.pri = (1/sigma.pri)^2
  resultado = -n*theta + Y*log(theta) - 1/2 * tau.pri * (log(theta) - mu.pri)^2
  resultado[is.nan(resultado)] = -Inf
  return(resultado)
}


#---------------------------------------
# dados para Bodó
#---------------------------------------
n = c(12,12,10,9,10,7,6)
Y = c(1,3,0,1,0,0,0)

# dá nomes aos índices
names(n) = seq(15,45,5)
names(Y) = seq(15,45,5)

estimativas.MV = Y/n

theta.apost = rep(0,7)
names(theta.apost) = seq(15,45,5)

for (x in seq(15,45,5)) {
  xchar     = as.character(x)
  val.theta = seq(0,.40,.001)
  log.post  = logP(val.theta, n[xchar], Y[xchar], mu.priori[xchar], sigma.priori[xchar])
  
  # muda escala (+c)
  log.post = log.post - max(log.post)
  
  Pmesmo = exp(log.post) / sum(exp(log.post))   # soma a 1
    
  plot( val.theta, Pmesmo, main=paste('A posteriori Bodó, Faixa Etária',x),
         type='l', col='blue', lwd=3)  
  
  abline(v=estimativas.MV[xchar], lty=2, lwd=2)
  
  # guardar o theta que tem a prob. max. a posteriori
  ix = which.max(Pmesmo)   # qual elemento de log.post é o maior
  
  theta.apost[xchar] = val.theta[ix]
  
} # for x

plot( seq(15,45,5), estimativas.MV, type='o', pch="V", main='Bodó RN \n Estimativas MV e MAP')
lines(seq(15,45,5), theta.apost, type='o', col='blue', lwd=4, pch='P')

#---------------------------------------
# dados para Natal
#---------------------------------------
n = c(3920,4400,4181,3653,3256,3279,3062)
Y = c(156,341,281,184,97,46,6)

# dá nomes aos índices
names(n) = seq(15,45,5)
names(Y) = seq(15,45,5)

estimativas.MV = Y/n

theta.apost = rep(0,7)
names(theta.apost) = seq(15,45,5)

for (x in seq(15,45,5)) {
  xchar     = as.character(x)
  val.theta = seq(0,.40,.001)
  log.post  = logP(val.theta, n[xchar], Y[xchar], mu.priori[xchar], sigma.priori[xchar])
  
  # muda escala (+c)
  log.post = log.post - max(log.post)
  
  Pmesmo = exp(log.post) / sum(exp(log.post))   # soma a 1
  
  plot( val.theta, Pmesmo, main=paste('A posteriori Natal, Faixa Etária',x),
        type='l', col='blue', lwd=3)  
  
  abline(v=estimativas.MV[xchar], lty=2, lwd=2)
  
  # guardar o theta que tem a prob. max. a posteriori
  ix = which.max(Pmesmo)   # qual elemento de log.post é o maior
  
  theta.apost[xchar] = val.theta[ix]
  
} # for x

plot( seq(15,45,5), estimativas.MV, type='o', pch="V", main='Natal RN \n Estimativas MV e MAP')
lines(seq(15,45,5), theta.apost, type='o', col='blue', lwd=4, pch='P')