disciplinas:ce227-2014-01:historico

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

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

Diferenças

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

--- disciplinas:ce227-2014-01:historico [2014/04/04 20:18]
paulojus
+++ disciplinas:ce227-2014-01:historico [2016/05/18 09:51] (atual)
paulojus
@@ Linha 14: / Linha 14: @@
 ===== Conteúdos das Aulas =====
-^ ^^ B & M ^^ Online ^
 ^ Data ^ Conteúdo ^ Leitura ^ Exercícios ^ Tópico ^
 | 12/02 Qua |Teorema de Bayes: revisão, interpretações e  generalização. Expressão probabilística de informação subjetiva, estimação baseada nos dados, estimação combinando informação prévia (subjetiva) e dados. Exemplo Binomial-Beta | [[#12/02|Ver abaixo]] | [[#12/02|Ver abaixo]] |  |
@@ Linha 31: / Linha 30: @@
 | 29/03 Sex |Discussão sobre a questão de 3a avaliação semanal | | |[[#26/03|Ver abaixo]] |
 | 02/03 Qua |Computação da questão de 3a avaliação semanal. Análise Bayesiana de dados: modelos linerares (Gaussianos). Construção do algorítmo de Gibbs | | |[[#02/04|Ver abaixo]] |
+| 04/04 Sex |Extensões e fundamentos do algorítmo de Gibbs. Gibbs metropolizado. Predição. Introdução do [[http:mcmc-jags.sourceforge.net|JAGS]] e rjags | | |[[#04/04|Ver abaixo]] |
+| 09/04 Qua |Inferência Bayesiana em modelos lineares generalizados (ex distribuição de Poisson). Introdução aos modelos de efeitos aleatórios - comparações com modelos de efeitos fixos e verossimilhança | | |[[#09/04|Ver abaixo]] |
+| 11/04 Sex |Modelos de efeitos aleatórios/hierárquicos e inferência Bayesiana. Algorítmos de inferência e INLA | | |[[#11/04|Ver abaixo]] |
+| 16/04 Qua |Avaliação semanal. Discussão da avaliação sobre modelos declarados em JAGS. Distribuições marginais <latex>[Y] = \int [Y|\theta][\theta] d\theta</latex> e interpretação. Exemplo: caso Poisson | |[[#16/04|Ver abaixo]] | |
+| 23/04 Qua |sem aula expositiva. Preparar e discutir [[#16/04|atividades propostas na última aula]]  para apresentação na próxima aula. | | | |
+| 25/04 Sex |Avaliação: apresentação das análises de dados correspondentes aos {{:disciplinas:ce227:ce227-av04.pdf|modelos da quarta avaliação semanal}}  | | | |
+| 30/04 Sex |continuação das apresentações e discussões  | | | |
+| 07/05 Qua |Modelos para efeitos aleatórios dependentes. Estrutura do modelo, comparação com outras estratégias/modelos. Exemplo: série temporal para dados binários.  | | |[[#07/05|Ver abaixo]] |
+| 09/05 Sex |Fundamentos do INLA | | |{{:disciplinas:ce227:apresentacao-inla.pdf|Apresentação}} |
+| 14/05 Sex |Fundamentos do INLA - II - Comentários sobre a (excelente!) apresentação de Gianluca Baio | | |[[http://www.statistica.it/gianluca/Talks/INLA.pdf|A apresentação]] |
@@ Linha 131: / Linha 141: @@
 plot(rlG2[,2], type="l")
 plot(rlG2[,3], type="l")
-plot(density(rlG[,1])); abline(v=coef(reg)[1])
+plot(density(rlG2[,1])); abline(v=coef(reg)[1])
-plot(density(rlG[,2])); abline(v=coef(reg)[2])
+plot(density(rlG2[,2])); abline(v=coef(reg)[2])
-plot(density(rlG[,3])); abline(v=summary(reg)$sigma^2)
+plot(density(rlG2[,3])); abline(v=summary(reg)$sigma^2)
 </code>
@@ Linha 139: / Linha 149: @@
 === 04/04 ===
   - Estender os códigos anteriores para incluir a predição de novos valores
+<code R>
+#### código agora incluindo predições
+reg.lin.GIBBS <- function(y, x, Nsim, iniBeta, x.pred = NULL){
+    require(MASS)
+    n <- length(y)
+    X <- cbind(1, x)
+    XX <- crossprod(X)
+    bhat <- solve(XX,crossprod(X,y))
+    ##
+    NC <- ncol(X) + 1
+    if(!is.null(x.pred)){
+        Xpred <- cbind(1, x.pred)
+        NC <- NC + nrow(Xpred)
+    }
+    sim <- matrix(0, nrow = Nsim, ncol=NC)
+    sim[1,1:2] <- iniBeta
+    sim[1, 3]  <- 1/rgamma(1, shape=(n-2)/2, scale=2/crossprod(y-X%*%sim[1,1:2]))
+    ##
+    for(i in 2:Nsim){
+        sim[i, 3]  <- 1/rgamma(1, shape=(n-2)/2, scale=2/crossprod(y-X%*%sim[(i-1),1:2]))
+        sim[i,1:2] <- mvrnorm(n=1, mu=bhat, Sigma = solve(XX)*sim[i, 3])
+        if(!is.null(x.pred))
+            sim[i,-(1:3)] <- rnorm(nrow(Xpred), m=Xpred%*%sim[i,1:2], sd=sqrt(sim[i,3]))
+    }
+    return(sim)
+}
+## agora com predição
+rlG0 <- reg.lin.GIBBS(y=y, x=x, Nsim=15000, iniBeta=c(0,0.6), x.pred=1:20)
+dim(rlG0)
+rlG1 <- rlG0[-(1:3000),]
+dim(rlG1)
+rlG2 <- rlG1[10*(1:1200),]
+dim(rlG2)
+rbind(dim(rlG0), dim(rlG1), dim(rlG2))
+reg
+summary(reg)$sigma^2
+par(mfrow=c(4,5))
+apply(rlG2[,-(1:3)], 2, function(x) plot(density(x),type="l"))
+par(mfrow=c(1,1))
+y.pred <- apply(rlG2[,-(1:3)], 2, mean)
+plot(y ~ x)
+abline(reg)
+lines(1:20, y.pred, col=2)
+#
+y.pred <- apply(rlG1[,-(1:3)], 2, mean)
+plot(y ~ x)
+abline(reg)
+lines(1:20, y.pred, col=2)
+</code>
+  * Incluir bandas de predição usuais e bayesianas no gráfico
+  * Generalizar código para família de priori Normal-GammaInversa
+  * Verificar o efeito/sensitividade à especificação das prioris
+**Exemplos JAGS/rjags**
+  - Média e variância (precisão) da distribuição normal<code R>
+n <- 20
+x <- rnorm(n, 70, 5)
+#write.table(x,
+#            file = 'normal.data',
+#            row.names = FALSE,
+#            col.names = FALSE)
+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"
+)
+require(rjags)
+## 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))
+jags <- jags.model('normal.modelo',
+                   data = list('x' = x, 'n' = n),
+                   n.chains = 3,
+                   inits = inis,
+                   n.adapt = 100)
+#update(jags, 1000)
+#sam <- jags.samples(jags, c('mu', 'tau'), 1000)
+sam <- coda.samples(jags, c('mu', 'tau'), n.iter=10000, thin=10)
+par(mfrow=c(2,2))
+plot(sam)
+str(sam)
+summary(sam)
+HPDinterval(sam)
+</code>
+  - regressão linear simples<code R>
+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")
+## alternativa:
+## tau ~dgamma(0.001, 0.001)
+## sigma2 <- 1/tau
+n <- 20
+x <- sort(runif(n, 0, 20))
+epsilon <- rnorm(n, 0, 2.5)
+y <- 2 + 0.5*x + epsilon
+#write.table(data.frame(X = x, Y = y, Epsilon = epsilon),
+#            file = 'reglin.dados',
+#            row.names = FALSE,
+#            col.names = TRUE)
+require(rjags)
+# require(R2jags)
+# require(dclone)
+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))
+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)
+sam <- jags.samples(jags,
+             c('b0', 'b1', 'sigma'),
+)
+class(sam)
+sam <- coda.samples(jags,
+             c('b0', 'b1', 'sigma'),
+)
+class(sam)
+str(sam)
+plot(sam)
+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>
+=== 09/04 ===
+  - Fazer uma análise Bayesiana de uma regressão de Poisson e comparar com resultados de um MLG (modelo linear generalizado) usual. Utilize um conjunto de dados qualquer da literatura (sugestão: conjunto ''gala'' do pacote "faraway")
+  - Use os dados do artigo //Confiabilidade e Precisão na Estimação de Médias// de Singer et al na [[http://www.rbes.ibge.gov.br/c/document_library/get_file?uuid=007b1342-7715-4539-9a26-e53cd9dc220b&groupId=37690208|Revista Brasileira de Estatística, v73, n. 236, jan./jun 2012]] e:
+    - ajuste um modelo "ingênuo" de média e resíduo
+    - ajuste um modelo de média, efeito de local e resíduo (medidas repetidas com efeitos fixos)
+    - ajuste um modelo de efeitos aleatórias
+    - ajuste um modelo Bayesiano.
+Compare e discuta os resultados.
+=== 11/04 ===
+<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}
+##
+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)
+## 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>
+=== 16/04 ===
+  - (Individual ou em duplas) Encontre um conjunto de dados para cada um dos casos na avaliação semanal de 16/04. Proceda análises não Bayesianas e Bayesianas e discuta os resultados. Trazer //scripts// para discussão na aula de 23/04 e apresentação em 25/04.
+  - Considere o modelo de verossimilhança <latex>[Y|\mu, \sigma^2] \sim N(\theta, \sigma^2)</latex> e a priori <latex>\tau = 1/\sigma^2 \sim Ga(a, b)</latex>. Mostre como obter a densidade: \\ <latex>[Y|\theta, a, b] = \frac{\Gamma((n/2)+a)}{\pi^{n/2} \Gamma(a) (\sum_i (x_i - \theta)^2 + 2b)^{(n/2)+a}}</latex>. \\ Como este resultado pode ser interpretado?
+=== 07/05 ===
+<code R>
+require(INLA)
+##
+## Visualizado dados
+##
+data(Tokyo)
+head(Tokyo)
+plot(y ~ time, data=Tokyo)
+## colocando na forma de proporção de dias com chuva
+plot(y/2 ~ time, data=Tokyo)
+##
+## 1. Modelo "Nulo": só intercepto
+## estimando a probabilidade de chuva como uma constante:
+fit.glm <- glm(cbind(y, n-y) ~ 1, family=binomial, data=Tokyo)
+abline(h=exp(coef(fit.glm))/(1+exp(coef(fit.glm))), col=2, lty=3, lwd=3)
+## como, como este modelo, todos os valores preditos são iguais bastaria fazer:
+abline(h=fitted(fit.glm)[1], col=2, lty=3, lwd=3)
+##
+## 2. Modelo com probabilidades variando no tempo (variável/processo latente)
+## modelando o logito(probabilidade) como um efeito aleatório correlacionado no tempo
+## segundo um "random walk" cíclico de ordm 2
+modelo = y ~ 0 + f(time, model="rw2", cyclic=T, param=c(1, 0.0001))
+fit <- inla(modelo, data=Tokyo, family="binomial", Ntrials=n, control.predictor=list(compute=TRUE))
+##
+names(fit)
+head(fit$summary.fitted.values)
+## sobrepondo ao gráfico dos dados (moda, mediana e médi são praticamente indistinguíveis
+with(fit, matlines(summary.fitted.values[,c(1,3:6)], lty=c(1,2,2,2,3), col=1))
+##
+## 3. Agora o mesmo modelo nulo anterior porém jusado pelo INLA
+##
+modelo0 = y ~ 1
+fit0 <- inla(modelo0, data=Tokyo, family="binomial", Ntrials=n, control.predictor=list(compute=TRUE))
+summary(fit0)
+fit0$summary.fitted.values
+with(fit0, matlines(summary.fitted.values[,c(1,3:6)], lty=c(1,2,2,2,3), col=2))
+##
+## 4. Modelando usando GAM (generalised additive model)
+##
+require(mgcv)
+fit.gam <- gam(cbind(y, n-y) ~ s(time), family=binomial, data=Tokyo)
+names(fit.gam)
+fitted(fit.gam, se=T)
+pred.gam <- predict(fit.gam, type="response", se=T, newdata=Tokyo["time"])
+names(pred.gam)
+with(pred.gam, matlines(cbind(fit, fit+2*se.fit, fit-2*se.fit), lty=c(1,2,2), col=4))
+</code>

Diferenças

Navegação

Busca

Ferramentas

QR Code