====== 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. ===== Conteúdos das Aulas ===== ^ Data ^ Conteúdo ^ Leitura ^ Exercícios ^ Tópico ^ | 19/02 Seg |Informações sobre o curso. Motivação inicial pela discussão de 3 situações |Cap 1 da "apostila" do curso |[[#19/02|Ver abaixo e sugerido em aula]] | | | 21/02 Qua |Discussão dos problemas propostos na aula anterior. Teorema de Bayes para o casos discretos e contínuos. Exemplo verossimilhança binomial e priori Beta. Obtenção de priori a partir da "opinião" sobre a proporção. |Cap 2 da "apostila" do curso, até o exemplo 2.2 |[[#21/02|Ver abaixo e sugerido em aula]] | | 26/02 Seg |Revisão de conceitos e exemplos discutidos até aqui. Discussão sobre as versões discreta e continua do Teorema de Bayes. Relações entre problemas e notação unificada. Especificação de priori Beta a partir da informação subjetiva |Cap 2 da "apostila" do curso |[[#26/02|Ver abaixo e sugerido em aula]] | | 28/02 Qua |1a avaliação periódica. Discussão da avaliação. Características da abordagem bayesiana: atualização sequencial, suficiência, princípio da verossimilhança |Cap. 2 (todo) |[[http://www.leg.ufpr.br/~paulojus/CE227/ce227-av01-2018.pdf|1a avaliação periódica]] | | 05/03 Seg |Elicitação de priori para o exemplo de verossimilhança binomial com priori beta - apresentação de implementação. Distribuição marginal e relações com modelos de efeitos aleatórios. Especificação de prioris: conjugadas, impróprias e representações de ignorância |Cap. 3 até 3.4 | | | 07/03 Qua |Discussão sobre o Capitulo 3 do material |Cap. 3 |[[#07/03|Ver abaixo e sugerido em aula]] | | 12/03 Seg |Implementação de problemas de aula e materiais (Alcides, Bruna e Hektor). Apresentação inicial de exercício de inferência sobre a variância. | |{{ :disciplinas:ce227:ce227-av01.pdf |Exercício para estudo e discussão }} | | 14/03 Qua |Inferência para posterioris de forma desconhecida: aproximação normal, discretização e amostragem por MCMC. | |[[#14/03|ver abaixo]] | | 19/03 Seg |Revisão e dúvidas sobe materia da aula anterior. Implementação da discretização. Amostragem e "jitter" das amostras da discreta | |[[#14/03|ver abaixo]] | | 21/03 Qua |Não haverá aula expositiva |Revisar conteúdos até aqui. | | | 26/03 Seg |Resumindo posteriori: Decisão (espaço discreto), estimação pontual (espaço contínuo), intervalos e testes |Cap 5, sec 5.1 |[[#26/03|ver abaixo]] | | 28/03 Qua |2a avaliação | | | | 02/04 Seg |Discussão da 2a avaliação | | | | 04/04 Qua |Resumindo posteriori: estimação pontual (espaço contínuo), intervalos e testes, Predição Bayesiana |Cap 5 e Cap 6 |[[#04/04|ver abaixo]] | | 09/04 Seg |Inferência em problemas com mais de um parâmetro |Cap 4: ler, estudar e refazer exemplos |[[#09/04|ver abaixo]] | | 11/04 Qua |Sem aula expositiva. Fazer atividades recomendadas da aula anterior ([[#09/04|ver abaixo]]) | | | | 16/04 Seg |Resolução e discussão do exercício 6.1. Obtendo a preditiva: (i) analiticamente, (ii) por aproximação normal (iii) por simulação | |[[#16/04|ver abaixo]] | | 18/04 Qua |Algoritmo amostrador de Gibbs (Gibbs sampler). Exemplo na inferência para distribuição normal | |[[#18/04|ver abaixo]] | | 23/04 Seg |Revisão Gibbs sampler. Modelo Poisson com priori Gamma e hiperpriori InvGamma. Derivação da posteriori, condicionais completas e implementação do algoritmo de Gibbs. Regressão linear: expressões para amostragem exata e via Gibbs| |[[#23/04|ver abaixo]] | | 25/04 Qua |Gibbs sampler compasso Metrópolis. Modelo Poisson com priori Normal. Derivação da posteriori, condicionais completas e implementação do algoritmo de Gibbs com um passo metrópolis. | |[[#23/04|ver abaixo]] | | 30/04 Seg |Feriado | | | | 02/05 Qua |1a prova | | | | 07/05 Seg |Discussão das questões da 1a prova. Definição de atividades para sequencia do curso | |[[#07/05|ver abaixo]] | | 09/05 Qua |Recursos computacionais para inferência Bayesiana - Atividades indicadas na aula anterior. Sem aula expositiva. | | | | 14/05 Seg |Gibbs sampler: exemplo das minas de carvão. Programação e utilização do JAGS | |[[#14/05|ver abaixo]] | | 16/05 Qua |Modelo de componentes de variância e correlação intraclasse, análise Bayesiana e não Bayesiana. Análise via JAGS | |[[#16/05|ver abaixo]] | | 21/05 Seg |Sem aula expositiva: semana dedicada às atividades do RDay e RBras | | | | 23/05 Qua |Sem aula expositiva: semana dedicada às atividades do RDay e RBras | | | | 28/05 Seg |Sem aula expositiva: interrupção de aulas na UFPR | | | | 30/05 Qua |Sem aula expositiva: interrupção de aulas na UFPR | | | | 04/06 Seg |Reorganização do curso. Dúvidas e perguntas dos alunos | | | | 06/06 Qua |2a prova. Toda matéria discutida até aqui | | | | 11/06 Qua |Bayesiano empírico - modelo Poisson-Gamma para taxas | | | | 13/06 Qua |Apresentações | | | | 18/06 Qua |Apresentações | | | | 20/06 Qua |3a prova. | | | === 19/02 === - **Problema 1:** catapora ou varíola? Formalizar e responder pergunta de interesse - 90% dos que tem varíola apresentam os sintomas reportados - 80% dos que tem catapora apresentam os sintomas reportados - prevalência de varíola na população: 1/1000 - prevalência de catapora na população: 1/1000 - **Problema 2:** - Assistir a parte inicial do [[https://www.youtube.com/watch?v=Ozpek_FrOPs|//sketch// four candles/fork handles]] e notar o problema do reconhecimento de voz - Relacionar elementos com o problema anterior - Formalizar notações, atribuir probabilidades e calcular quantidades de interesse - **Problema 3:** (está na apostila mas **tentar resolver de alguma forma** antes de ler o material) - 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? - Ler e resolver exercícios do Capítulo 1 da apostila === 21/02 === - Revisitar **Problema 2:** acima - 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. - 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) ^ Estimativa ^ Intervalo ^ Probabilidade ^ | 0,63 | (0,40 ; 0,75) | 90% | | 0,42 | (0,32 ; 0,52) | 80% | | 0,20 | (0,05 ; 0,35) | 80% | | 0,50 | (0 ; 1) | 100% | | 0,30 | (0,20 ; 0,40) | 50% | === 26/02 === - Completar problemas propostas nas aulas anteriores após as discussões em aula - Escrever um código para o Exemplo da Poisson (2.3 do material), que permita desenhar as funções e avaliar efeitos de prioris e dados - Ler e resolver exercícios do Capítulo 2 da apostila === 07/03 === - Exercícios do Cap 3 - Escrever um código que receba: modelo, dados, priori conjugada e retorne posteriori (parâmetros e gráficos) === 14/03 === - Reler o material do {{ :disciplinas:ce227:ce227-av01.pdf |Exercício para estudo e discussão }} - Montar um algoritmo para aproximação por discretização do exemplo - Tome algum outro modelo de um parâmetro e desenvolva os resultados análogos aos vistos em aula. === 19/03 === - Visualizar, experimentar e comentar o [[https://brunaw.shinyapps.io/bayes_app/|aplicativo shiny]] construído por Bruna, Hektor e Alcides === 26/03 === - Refazer exemplos e fazer Exercício 5.1 do Cap 5 - 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)) === 02/04 === - Na questão 1 verificar como a mudança na priori (proporção de motoristas acima do limite) afeta os resultados - Na questão 2 refazer com a parametrização alternativa da exponencial e gama - Ainda na questão 2 fazer utilizando o resultado genérico de prioris conjugadas para família exponencial - 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 === 04/04 === - Refazer exemplos e fazer Exercício 5.2 a 5.4 do Cap 5 - Refazer exemplos e fazer Exercícios cap 6 - Escrever funções mostrando como média, mediana e quartis podem ser calculados a partir de minimização de função perda: - para um conjunto de dados - para uma distribuição discreta - para uma distribuição contínua === 09/04 === - 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. - 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. - No Exemplo 3 obter a marginal [\sigma^2|y] e a posteriori condicional [\mu|\sigma^2,y] - Ainda no exemplo 3 definir os hiperparâmetros de obter uma simulação de dados do modelo - Com os dados simulados obter as expressões da posteriori conjunta, marginais (do material) e condicional conjunta (item anterior) - Obter uma simulação da posteriori. Comparar a conjunta e marginais teórica e simulada. === 16/04 === - 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) - Experimentar diferentes - 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) === 18/04 === 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 provou-se 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) === 23/04 === - Implementar modelo semelhante ao visto em aula porém com log(lambda ~Normal). (ver detalhes na versão revisada do Cap 8 do material do curso. - Implementar a regressão linear via algoritmo de Gibbs. Usar dados simulados de uma regressão linear simples. Incluir amostras da preditiva no algoritmo - 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) === 07/05 === - **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 - **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 - **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 === 14/05 === - {{:disciplinas:ce227:changepointjags.r|Script R/JAGS para análise dos dados do Cap 8}} (changepoint Poisson) === 16/05 === - 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:** - Complementar as análise acima com exploração dos resultados, obtenção de gráficos e resultados de interesse - 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**\\ [[http://www.rbes.ibge.gov.br/images/doc/rbe_236_jan_jun2012.pdf|Revista Brasileira de Estatística, v73]], n. 236, jan./jun. 2012. - 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]