disciplinas:ce227-2016-01:historico

Não foi possível enviar o arquivo. Será algum problema com as permissões?

Você está aqui: start » disciplinas » ce227-2016-01 » historico

Diferenças

Aqui você vê as diferenças entre duas revisões dessa página.

--- disciplinas:ce227-2016-01:historico [2016/05/12 17:40]
paulojus
+++ disciplinas:ce227-2016-01:historico [2016/05/18 09:38]
paulojus
@@ Linha 37: / Linha 37: @@
 | 09/05 Seg |Discussão da prova e detalhamento do problema da questão 5  | | |[[#09/05|Ver abaixo]] |
 | 11/05 Qua |Discussão Caps 7 e 8 do material | | |[[#11/05|Ver abaixo]] |
+| 16/05 Seg |Inferência Bayesiana utilizando o JAGS - instalação e exemplos | | |[[#16/05|Ver abaixo]] |
 === 29/02 ===
 Manifestar uma opinião subjetiva sobre o parâmetro de uma distribuição binomial. (basear-se no contexto  de intenção de voto discutido em aula)
@@ Linha 62: / Linha 65: @@
   - Derivar os expressões das condicionais completas no problema do ponto de mudança da Poisson (ex. do capitulo 8)
   - Implementar o algorítmo de Gibbs para este exemplo.
+=== 16/05 ===
+Exemplos discutidos utilizando JAGS/rjags:
+  - Amostragem da normal <code R>
+## Simulando um conjunto de dados
+n <- 20
+x <- rnorm(n, 70, 5)
+## Exportar os dados (não é necessário) se utilizando o rjags
+#write.table(x,
+#            file = 'normal.data',
+#            row.names = FALSE,
+#            col.names = FALSE)
+## Especificação do modelo (deve ser exportada para um arquivo)
+cat( "model {
+	for (i in 1:n){
+		x[i] ~ dnorm(mu, tau)
+	}
+	mu ~ dnorm(0, 0.0001)
+	tau <- pow(sigma, -2)
+	sigma ~ dunif(0, 100)
+ }", file="normal.modelo"
+)
+## Carregando o pacotes rjags (pode-se ainda usar outros como runjags, R2jags etc)
+require(rjags)
+## Definindo valores iniciais. No caso três conjuntos porque iremos rodas 3 cadeias.
+## OBS: valores iniciais são dispensáveis neste exemplo
+inis <- list(list(mu=10, sigma=2),
+             list(mu=50, sigma=5),
+             list(mu=70, sigma=10))
+## O proximo comando prepara e " compila" o modelo e opções para o algorítmo
+jags <- jags.model('normal.modelo',
+                   data = list('x' = x, 'n' = n),
+                   n.chains = 3,
+                   inits = inis,
+                   n.adapt = 100)
+## Obtendo as amostras (diferentes opções, a última já prepara em formato para uso com o
+##  pacote ´coda´)
+#update(jags, 1000)
+#sam <- jags.samples(jags, c('mu', 'tau'), 1000)
+sam <- coda.samples(jags, c('mu', 'tau'), n.iter=10000, thin=10)
+## Visualizações e resultados
+par(mfrow=c(2,2))
+plot(sam)
+str(sam)
+summary(sam)
+HPDinterval(sam)
+</code>
+  - regressão linear simples<code R>
+## simulando dados
+n <- 20
+x <- sort(runif(n, 0, 20))
+epsilon <- rnorm(n, 0, 2.5)
+y <- 2 + 0.5*x + epsilon
+plot(y ~ x)
+lines(lm(y ~x))
+## especificando o modelo para o JAGS
+cat( "model {
+      	for (i in 1:n){
+		y[i] ~ dnorm(mu[i], tau)
+		mu[i] <- b0 + b1 * x[i]
+	}
+	b0 ~ dnorm(0, .0001)
+	b1 ~ dnorm(0, .0001)
+	tau <- pow(sigma, -2)
+	sigma ~ dunif(0, 100)
+}", file="reglin.modelo")
+## poderia-se redefinir o modelo acima com uma possível priori alternativa, por ex:
+## tau ~ dgamma(0.001, 0.001)
+## sigma2 <- 1/tau
+#write.table(data.frame(X = x, Y = y, Epsilon = epsilon),
+#            file = 'reglin.dados',
+#            row.names = FALSE,
+#            col.names = TRUE)
+require(rjags)
+## Valores iniciais (vamos rodar só duas cadeias neste exemplo)
+inis <- list(list(b0=0, b1=1, sigma=1),
+             list(b0=1, b1=0.5, sigma=2),
+             list(b0=2, b1=0.1, sigma=5))
+## Compilando modelo, dados e opções
+jags <- jags.model('reglin.modelo',
+                   data = list('x' = x,
+                               'y' = y,
+                               'n' = n),
+                   n.chains = 2,
+                   # inits=inits,
+                   n.adapt = 100)
+#update(jags, 1000)
+class(jags)
+## obtenção das amostras da posteriori ...
+## ... via rjags
+sam <- jags.samples(jags,
+                   c('b0', 'b1', 'sigma'),
+)
+class(sam)
+## ... ou via coda
+sam <- coda.samples(jags,
+             c('b0', 'b1', 'sigma'),
+)
+class(sam)
+str(sam)
+plot(sam)
+## Pode-se tb obter as distribuições preditivas correspondentes a cada observação
+sam <- coda.samples(jags,
+             c('b0', 'b1', 'sigma', "y"),
+)
+str(sam)
+int <- HPDinterval(sam)
+str(int)
+## complementar com gráficos, resumos, inferências de interesse, etc
+</code>
+  - Coeficiente de correlação  intraclasse <code R>
+## 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])
+</code>
+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.

Diferenças

Navegação

Busca

Ferramentas

QR Code