oeste_potiguar.R

graphics.off()
windows(record=TRUE)

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

## distrib de f20 nos munic bras com > 500 mulheres na faixa etária 20-24
Brasil = read.csv( file='http://ufrn.schmert.net/dados/fecund.2010.csv',
            stringsAsFactors = FALSE) %>%
            filter(Mu20 > 500) %>%
            mutate(f20 = Cr20/Mu20)

plot(density(Brasil$f20), 
     main='5f20 nos municípios com 500+ mulheres 20-24',
     xlab='5f20', ylim=c(0,25), lwd=2)

x = seq(from=0, to=.30, by=.001 )
lines(x, dbeta(x, 17, 180), col='red', lwd=2)
text(.15,  5, 'distrib empírica', col='black')
text(.15, 15, 'approx beta(17,180)', col='red')

######## estimativas para 3 municípios no RN

munic.sel = c('Augusto Severo', 'Pilões', 'Natal')

RN = read.csv( file='http://ufrn.schmert.net/dados/fecund.2010.csv',
               stringsAsFactors = FALSE) %>%
        filter(Sigla=='RN', Nome_Munic %in% munic.sel) %>%
        mutate( f20=Cr20/Mu20) %>%
        select(Nome_Munic, Mu20, Cr20, f20)

RN

##       Nome_Munic Mu20 Cr20    f20
## 1 Augusto Severo   40    1 0.0250
## 2          Natal 4400  341 0.0775
## 3         Pilões   16    0 0.0000

theta = seq(from=0, to=.30, by=.001)

densidade.a.priori = dbeta(theta, 17, 180)

for (i in 1:nrow(RN)) {
  
  este.nome   = RN$Nome_Munic[i]
  este.Cr     = RN$Cr20[i] 
  este.Mu     = RN$Mu20[i]
  
  Vrel        = dbinom(x=este.Cr, size=este.Mu, prob=theta )
  este.veross = 1000*Vrel/sum(Vrel)
  
  densidade.a.posteriori = dbeta(theta, 17+este.Cr, 180+este.Mu-este.Cr)
  
  
  plot( theta, este.veross, type='l', lty=2, lwd=2,main=este.nome,
          ylim=range(0,densidade.a.posteriori, densidade.a.priori,
                     este.veross))
  lines( theta, densidade.a.priori, type='l', col='red', lwd=2)
  lines( theta, densidade.a.posteriori, type='l', col='seagreen', lwd=2)
  
  abline(v= este.Cr/este.Mu)
  }