CE-227 - Primeiro semestre de 2018

No quadro abaixo será anotado o conteúdo dado em cada aula do curso.
São indicados os Capítulos e Sessões correspondentes nas referências bibliográficas, bem como os exercícios sugeridos.

Veja ainda depois da tabela as Atividades Complementares.

Observação sobre exercícios recomendados os exercícios indicados são compatíveis com o nível e conteúdo do curso.
Se não puder fazer todos, escolha alguns entre os indicados.

1. Problema 1: catapora ou varíola? Formalizar e responder pergunta de interesse
   1. 90% dos que tem varíola apresentam os sintomas reportados
   2. 80% dos que tem catapora apresentam os sintomas reportados
   3. prevalência de varíola na população: 1/1000
   4. prevalência de catapora na população: 1/1000
2. Problema 2:
   1. Assistir a parte inicial do //sketch// four candles/fork handles e notar o problema do reconhecimento de voz
   2. Relacionar elementos com o problema anterior
   3. Formalizar notações, atribuir probabilidades e calcular quantidades de interesse
3. Problema 3: (está na apostila mas tentar resolver de alguma forma antes de ler o material)
   1. Uma caixa possui 6 bolas. Retiram-se 3 que são todas pretas. Qual a probabilidade de não haver mais bolas pretas na caixa?
4. Ler e resolver exercícios do Capítulo 1 da apostila

1. Revisitar Problema 2: acima
2. Considerar o exemplo de aula com dados n=200, y=75 e as prioris escolhidas. Obter gráficos das prioris e posterioris. Verificar o efeito do tamanho de amostra aumentando n e y na mesma proporção e repetindo os gráficos.
3. Escrever algum código para obtenção a priori a partir da opinião sobre a proporção cf discutido em aula. Verificar seu código com a sua opinião e as ilustradas em aula (ver tabela a seguir)

1. Completar problemas propostas nas aulas anteriores após as discussões em aula
2. Escrever um código para o Exemplo da Poisson (2.3 do material), que permita desenhas as funções e avaliar efeitos de prioris e dados
3. Ler e resolver exercícios do Capítulo 2 da apostila

1. Exercícios do Cap 3
2. Escrever um código que receba: modelo, dados, priori conjugada e retorne posteriori (parâmetros e gráficos)

1. Reler o material do Exercício para estudo e discussão
2. Montar um algoritmo para aproximação por discretização do exemplo
3. Tome algum outro modelo de um parâmetro e desenvolva os resultados análogos aos vistos em aula.

1. Visualizar, experimentar e comentar o aplicativo shiny construído por Bruna, Hektor e Alcides

1. Refazer exemplos e fazer Exercício 5.1 do Cap 5
   1. Código para o Exemplo 5.1:

      # Priori [\theta]
      (th <- c(th1=0.6, th2=0.4))
       
      # Verossimilhança  [Y|\theta]
      y.th1 <- c(y1=0.35, y2=0.30, y3=0.21, y4=0.14) 
      y.th2 <- c(y1=0.09, y2=0.17, y3= 0.25, y4=0.49) 
      (y.th <- rbind(y.th1, y.th2))
       
      # Conjunta [Y, \theta] = [\theta] \cdot [Y|\theta]
      (yth <- th * y.th)
      rownames(yth) <- c("yth1", "yth2")
      yth 
       
      # Marginal [Y]
      (y <- drop(crossprod(th, y.th)))
      # ou ...
      colSums(th * y.th)
       
      # Posteriori
      (th.y <- t(t(yth)/drop(y)))
      rownames(th.y) <- c("th1.y", "th2.y")
      th.y
       
      ## Fc Perda
      L <- diag(c(8,20))
      rownames(L) <- paste("a", 1:2, sep="")
      L
      (L.th.y <- L %*% th.y)
       
      ## Decisão
      D.f <- function(x) ifelse(x[1] < x[2], "a1:Vacina", "a2:Não Vacina")
      apply(L.th.y, 2, D.f)
       
      ## Perda baseada na regra de decisão para cada resultado de exame
      apply(L.th.y, 2, min)
       
      ## Risco de Bayes (associado a uma determinada regra adotada aqui)
      ## -- perda média esperada
      sum(apply(L.th.y, 2, min) * y)
       
      ## Outra regra: Vacina todo mundo
      L <- diag(c(8,0))
      (LT.th.y <- LT %*% th.y)
      sum(apply(LT.th.y, 2, sum) * y)
      ## ou simplesmente...
      sum(th * c(8,0))
       
      ## Ainda outra regra: não vacina ninguém
      LN <- diag(c(0,20))
      (LN.th.y <- LN %*% th.y)
      sum(apply(LN.th.y, 2, sum) * y)
      ## ou simplesmente...
      sum(th * c(0,20))

1. Na questão 1 verificar como a mudança na priori (proporção de motoristas acima do limite) afeta os resultados
2. Na questão 2 refazer com a parametrização alternativa da exponencial e gama
3. Ainda na questão 2 fazer utilizando o resultado genérico de prioris conjugadas para família exponencial
4. Na questão 4 supor uma amostra de 5 valores (7, 3, 4, 5, 2), obter a posteriori e fazer os gráficos de priori, verossimilhança e posteriori

1. Refazer exemplos e fazer Exercício 5.2 a 5.4 do Cap 5
2. Refazer exemplos e fazer Exercícios cap 6
3. Escrever funções mostrando como média, mediana e quartis podem ser calculados a partir de minimização de função perda:
   1. para um conjunto de dados
   2. para uma distribuição discreta
   3. para uma distribuição contínua

1. Fazer um código (com operações matriciais) para os cálculos do Exemplo 1. O código deve permitir definir diferentes prioris e verossimilhanças. Experimentar com valores diferentes do exemplo.
2. Especificar valores para os hiperparâmetros p e q no Exemplo 2 e simular um conjunto de dados. Obter a posteriori e maginais. Fazer gráficos conjuntos e marginais da priori e posteriori.
3. No Exemplo 3 obter a marginal e a posteriori condicional
4. Ainda no exemplo 3 definir os hiperparâmetros de obter uma simulação de dados do modelo
5. Com os dados simulados obter as expressões da posteriori conjunta, marginais (do material) e condicional conjunta (item anterior)
6. Obter uma simulação da posteriori. Comparar a conjunta e marginais teórica e simulada.

1. Para os exemplos e exercícios do Cap 6, obter a preditiva pelas 3 formas discutidas em aula. Escrever códigos que mostrem e comparem as preditivas (simlar ao visto em aula)
2. Experimentar diferentes
3. Segue código visto para ex 6.1

   ## Exercício 6.1
   ## Adicional:
   ## Seja uma amostra 7,5,8,9,3
   ## n=5 , soma = 33
   ## Seja a priori G(2, 2)
   ## A posteriori é G(2+33, 2+5) 
   ## A preditiva analítica é BN(2+33, (2+5)/(2+5+1))
    
   ## 1. Obtendo 1 simulação da preditiva
   ## Passo 1: simula valor do parâmetro da posteriori
   th <- rgamma(1, 35, 7)
   ## Passo 2: simula valor predito da verossimilhança 
   yp <- rpois(1, lam=th)
    
   ## 2. Obtendo 1000 simulações da preditiva
   ## Passo 1: simula valores do parâmetro da posteriori
   th <- rgamma(1000, 35, 7)
   ## Passo 2: simula valores predito da verossimilhança 
   yp <- rpois(1000, lam=th)
    
   ## Preditiva estimada por simulação
   table(yp)
   yp.sim <- table(yp)/1000
    
   ## Preditiva exata
   yp.teo <- dnbinom(0:14, size=35, prob=7/8)
    
   ## comparando 
   rbind(yp.sim, yp.teo)
   ## Pode-se aumentar o número de simulações para uma melhor predição
   th <- rgamma(1000, 35, 7)
   th <- rgamma(10000, 35, 7)
   yp <- rpois(10000, lam=th)
   yp.sim <- table(yp)/10000
   yp.teo <- dnbinom(0:max(yp), size=35, prob=7/8)
   rbind(yp.sim, yp.teo)
    
   ## Gráficos
   ## preditiva teórica (analítica)
   plot((0:17)-0.05, yp.teo, type="h")
   ## simulação da preditiva
   lines((0:17)+0.05, yp.sim, type="h", col=2)
   ## preditiva não bayesiana (plug-in)
   yp.nonB <- dpois(0:17, lam=33/5)
   lines((0:17)+0.15, yp.nonB, type="h", col=4)
   ## aproximação normal da preditiva
   curve(dnorm(x, m=5, sd=sqrt(5+35/49)), add=T)

Código visto em aula

##
## Inferência na distribuição normal
##
## Conjunta:
##f(\mu, \sigma^2|y) = (\sigma^2)^{\frac{n}{2}-1} \exp\left{-\frac{1}{2\sigma^2} (S^2 + n(\theta - \overline{y})) \right\}
##
## Condicionais
##    [\mu|\sigma^2, y] \sim {\rm N}(\overline{y}, \sigma^2/n)
##    [\sigma^2|\mu, y] \sim {\rm IG}(\frac{n}{2}, \frac{2}{A})
##
## Marginais
##    [\mu|y] \sim {\rm t}_{n-1}(\overline{y}, S^2/n)
##    \frac{\mu - \overline{y}}{\sqrt{sigma^2/n}} \sim {\rm t}_{n-1}
##    
##    [\sigma^2|y] \sim {\rm IG}(\frac{n-1}{2}, \frac{2}{S^2})
##    \frac{S^2}{\sigma^2} \sim \chi^2_{n-1}
##    
##    S^2 = \sum_{i=1}^{n} (y_i - \overline{y})^2
##    A = S^2 + n(\theta - \overline{y})^2
 
## Nos códigos abaixo S^2 é denotado por SQ
set.seed(20180419)
(y <- rnorm(12, mean=50, sd=8))
dados <- list(n=length(y), m=mean(y), v = var(y), SQ = sum((y-mean(y))^2))
##
## Amostra (exata) da posteriori
##
## para amostrar de pode-se explorar a fatoração:
## [\mu, \sigma^2|y] = [\sigma^2|y] \cdot [\mu|\sigma^2,y] = 
## ou, alternativamente
## [\mu, \sigma^2|y] = [\mu|y] \cdot [\sigma^2|\mu,y] = 
##
## Vamos adotar aqui a primeira fatoração:
## Obtendo uma amostra
##  (i) Amostrar \sigma^2 de [\sigma^2|y]
(sigma2.sim <- with(dados, 1/rgamma(1, shape=(n-1)/2, scale=2/SQ)))
## (ii) Amostrar \mu de [\mu |\sigma^2,y]
(mu.sim <- with(dados, rnorm(1, mean=m, sd=sqrt(sigma2.sim/n))))
## Obtendo 25.000 amostras
N <- 25000
sigma2.sim <- with(dados, 1/rgamma(N, shape=(n-1)/2, scale=2/SQ))
mu.sim <- with(dados, rnorm(N, mean=m, sd=sqrt(sigma2.sim/n)))
 
## Gráficos das amostras (correespondem às marginais)
par(mfrow=c(1,2))
t.sim <- with(dados, (mu.sim - m)/sqrt(v/n))
curve(dt(x, df=dados$n-1), from=-4, to=4)
lines(density(t.sim), col=4)
## note a diferença para uma distribuição normal:
curve(dnorm(x), from=-4, to=4, col=2, lty=3, add=TRUE)
 
chi.sim <- with(dados, SQ/sigma2.sim)
curve(dchisq(x, df=dados$n-1), from=0, to=40)
lines(density(chi.sim), col=4)
 
##
## Amostra (Gibbs) da posteriori
##
## A estratégia de Gibbs é alternar as simulações entre **as distribuições condicionais**
## o que "parece" errado ,as provouse que a cadeia de valores assim simulados **converge** para a distribuição conjunta 
##    [\mu|\sigma^2, y] \sim {\rm N}(\overline{y}, \sigma^2/n)
##    [\sigma^2|\mu, y] \sim {\rm IG}(\frac{n}{2}, \frac{2}{A})
## Obtendo uma amostra
## Como a distribuição de um parâmetro depende da distribuição do outro, 
## é necessário fornecer/arbitrar um valor para inicial o algoritmo
mu0 <- 50
##  (i) Amostrar \sigma^2 de [\sigma^2|\mu, y]
A <- with(dados, SQ + n*(mu0 - m)^2)
(sigma2.simG <- with(dados, 1/rgamma(1, shape=n/2, scale=2/A)))
## (ii) Amostrar \mu de [\mu |\sigma^2,y]
(mu.simG <- with(dados, rnorm(1, mean=m, sd=sqrt(sigma2.sim/n))))
 
## Gerando agora 25.000 amostras
N <- 25000
mu.simG <- sigma2.simG <- numeric(N)
mu.simG[1] <- 30
sigma2.simG[1] <- 100
 
{for(i in 2:N){
    A <- with(dados, SQ + n*(mu.simG[i-1]-m)^2)
    sigma2.simG[i] <- with(dados, 1/rgamma(1, shape=n/2, scale=2/A))
    mu.simG[i] <- with(dados, rnorm(1, mean=m, sd=sqrt(sigma2.simG[i]/n)))
 }
}
 
plot(mu.simG, type="l")
plot(mu.simG[-(1:1000)], type="l")
 
plot(sigma2.simG, type="l")
plot(sigma2.simG[-(1:1000)], type="l")
 
plot(log(sigma2.simG), type="l")
plot(log(sigma2.simG[-(1:1000)]), type="l")
 
par(mfrow=c(1,2))
t.sim <- with(dados, (mu.sim - m)/sqrt(v/n))
curve(dt(x, df=dados$n-1), from=-4, to=4)
lines(density(t.sim), col=4)
##curve(dnorm(x), from=-4, to=4, col=2, add=TRUE)
t.simG <- with(dados, (mu.simG - m)/sqrt(v/n))
lines(density(t.simG), col=3, lwd=2)
 
chi.sim <- with(dados, SQ/sigma2.sim)
curve(dchisq(x, df=dados$n-1), from=0, to=40)
lines(density(chi.sim), col=4)
chi.simG <- with(dados, SQ/sigma2.simG)
lines(density(chi.simG), col=3, lwd=2)

1. Implementar modelo semelhante ao visto em aula porém com <math>log(lambda ~Normal). (ver detalhes na versão revisada do Cap 8 do material do curso.
2. Implementar a regressão linear via algoritmo de Gibbs. Usar dados simulados de uma regressão linear simples. Incluir amostras da preditiva no algoritmo
3. Código para o modelo visto em aula:

   ## Simulando dados do modelo sendo estudado
   set.seed(2018)
   ctes <- list(a=3, c=2.5, d=0.8, n=50)
   with(ctes, EVIG(c, d))
   betas <- with(ctes, 1/rgamma(n, shape=c, scale=d))
   c(mean(betas),var(betas))
   lambdas <- with(ctes, rgamma(n, shape=a, rate=betas))
   (ctes$y <- rpois(ctes$n, lambda=lambdas))
   with(ctes, c(media=mean(y), var=var(y)))
   with(ctes, plot(prop.table(table(y)), type="h", ylim=c(0,0.3)))
   with(ctes,lines((0:max(y))+0.1, dpois(0:max(y), lambda=mean(y)), type="h", col=2))
   ##
   ## Iniciando inferência a ser feita via amostrador de Gibbs 
   ##
   ctes$sumY <- sum(ctes$y)
   ##
   N <- 11000  # número de simulação no algorítmo
   B <- 1000   # bunr-in - amostras s serem descartadas no início da cadeia
   beta.sam <- lambda.sam <- numeric(N) 
   beta.sam[1] <- lambda.sam[1] <- 10
   {
       for(i in 2:N){
           beta.sam[i] <- with(ctes, 1/rgamma(1, shape=a+c, scale=d/(d*lambda.sam[i-1]+1)))
           lambda.sam[i] <- with(ctes, rgamma(1, shape=ctes$a+sumY, scale=beta.sam[i]/(n*beta.sam[i]+1)))
       }
   }
    
   ## Explorando simulações
   par(mfrow=c(2,1))
   plot(beta.sam, type="l")
   plot(lambda.sam, type="l")
   ## retirando amostras consideradas aquecimento
   beta.sam <- beta.sam[-(1:B)]
   lambda.sam <- lambda.sam[-(1:B)]
   plot(beta.sam, type="l")
   plot(lambda.sam, type="l")
   plot(log(beta.sam), type="l")
   plot(lambda.sam, type="l")
    
   par(mfrow=c(1,2))
   plot(density(beta.sam)); abline(v=mean(betas)); rug(betas)
   plot(density(lambda.sam)); abline(v=mean(lambdas)); rug(lambdas)
   summary(ctes$y)
   summary(betas)
   summary(beta.sam)
   summary(lambdas)
   summary(lambda.sam)
    
   par(mfrow=c(1,2))
   plot(density(beta.sam, from=0, to=5)); abline(v=mean(betas)); rug(betas)
   plot(density(lambda.sam, from=0, to=20)); abline(v=mean(lambdas)); rug(lambdas)

1. *Atividade 1* (individual ou duplas) Buscar algum pacote do R ou outro programa que permita obter os resultados (analíticos) vistos até aqui no curso. Evitar coincidẽncias entre os escolhidos
2. *Atividade 2* (individual ou duplas) Buscar algum pacote do R ou outro programa que permita obter por simulação resultados pera os exemplos vistos até aqui no curso. Evitar coincidẽncias entre os escolhidos
3. *Atividade 3* (individual ou duplas) Utilizar o recurso visto na Atividade 2 para analizar algum modelo/exemplo não visto no curso. Evitar coincidẽncias entre os escolhidos

1. Script R/JAGS para análise dos dados do Cap 8 (changepoint Poisson)

1. Coeficiente de correlação intraclasse

   ## Dados simulados do modelo:
   ## Y_{ij} \sim N(\mu_{i}, \sigma^2_y)
   ##     mu_{i} = theta + b_{i}
   ##     b_{i} \sim N(0, \sigma^2_b)
   ## que, por ser normal (com ligação identidade)
   ## pode ser escrito por:
   ## Y_{ij} = \beta_0 + b_{i} + \epsilon_{ij} 
   ##
   ## simulando dados:
   Ngr <-  25
   Nobs <- 10
   set.seed(12)
   sim <- data.frame(id  = Ngr*Nobs,
                     gr  = rep(1:Ngr, each=Nobs),
                     bs  = rep(rnorm(Ngr, m=0, sd=10), each=Nobs),
                     eps = rnorm(Ngr*Nobs, m=0, sd=4)
                     )
   sim <- transform(sim, y = 100 + bs + eps)
   sim
    
   ## estimativas "naive"
   resumo <- function(x) c(media=mean(x), var=var(x), sd=sd(x), CV=100*sd(x)/mean(x))
   (sim.res <- aggregate(y~gr, FUN=resumo, data=sim))
   var(sim.res$y[,1])
   mean(sim.res$y[,2])
   mean(sim$y)
    
   ## A seguir serão obtidas inferências de três formas diferentes:
   ## - ajuste modelo de efeito aleatório (não bayesiano)
   ## - ajuste via JAGS (inferência por simulação da posteriori)
   ## - ajuste via INLA (inferência por aproximação da posteriori)
    
   ##
   ## Modelo de efeitos aleatórios
   ##
   require(lme4)
   fit.lme <- lmer(y ~ 1|gr, data=sim)
   summary(fit.lme)
   ranef(fit.lme)
   coef(fit.lme)$gr - fixef(fit.lme)
   print(VarCorr(fit.lme), comp="Variance")
    
   ## JAGS
   require(rjags)
    
   sim.lst <- as.list(sim[c("gr","y")])
   sim.lst$N <- nrow(sim)
   sim.lst$Ngr <- length(unique(sim$gr))
   mean(sim.lst$y)
    
   cat("model{
       for(j in 1:N){
           y[j] ~ dnorm(mu[gr[j]], tau.e)
        }
       for(i in 1:Ngr){
           mu[i] ~ dnorm(theta, tau.b)
       }
       theta ~ dnorm(0, 1.0E-6)
       tau.b ~ dgamma(0.001, 0.001)
       sigma2.b <- 1/tau.b
       tau.e ~ dgamma(0.001, 0.001)
       sigma2.e <- 1/tau.e
       cci <- sigma2.e/(sigma2.e+sigma2.b)
   }", file="sim.jags")
    
   sim.jags <- jags.model(file="sim.jags", data=sim.lst, n.chains=3, n.adapt=1000)
   ## inits = ...
    
   fit.jags <- coda.samples(sim.jags, c("theta", "sigma2.b", "sigma2.e", "cci"), 10000, thin=10)
    
   summary(fit.jags)
   plot(fit.jags)
    
   ##
   require(INLA)
    
   fit.inla <- inla(y ~ f(gr) , family="gaussian", data=sim)
   summary(fit.inla)
   sqrt(1/fit.inla$summary.hyperpar[,1])

Atividades propostas:

1. Complementar as análise acima com exploração dos resultados, obtenção de gráficos e resultados de interesse
2. Ajustar o modelo acima aos dados de:
   Julio M. Singer, Carmen Diva Saldiva de André, Clóvis de Araújo Peres
   Confiabilidade e Precisão na Estimação de Médias
   Revista Brasileira de Estatística, v73, n. 236, jan./jun. 2012.
3. Identificar e ajustar modelos (não bayesianos, bayesianos por simulação ou aproximados) para dados simulados da seguinte forma:

   set.seed(123456L)
   n <- 50
   m <- 10
   w <- rnorm(n, sd=1/3)
   u <- rnorm(m, sd=1/4)
   b0 <- 0
   b1 <- 1
   idx <- sample(1:m, n, replace=TRUE)
   y <- rpois(n, lambda = exp(b0 + b1 * w + u[idx]