Não foi possível enviar o arquivo. Será algum problema com as permissões?
Diferenças
Aqui você vê as diferenças entre duas revisões dessa página.
Ambos lados da revisão anterior Revisão anterior Próxima revisão | Revisão anterior Próxima revisão Ambos lados da revisão seguinte | ||
disciplinas:ce227-2016-01:historico [2016/04/28 16:48] paulojus |
disciplinas:ce227-2016-01:historico [2016/05/18 09:44] paulojus |
||
---|---|---|---|
Linha 31: | Linha 31: | ||
| 18/04 Seg |Predição Bayesiana |Cap 6 do material do curso | |[[#18/04|Ver abaixo]] | | | 18/04 Seg |Predição Bayesiana |Cap 6 do material do curso | |[[#18/04|Ver abaixo]] | | ||
| 20/04 Qua |Testes FBST - parte 1/2 | | | | | | 20/04 Qua |Testes FBST - parte 1/2 | | | | | ||
- | | 25/04 Qua |Testes FBST - parte 2/2 e revisão/dúvidas para prova | | | | | + | | 25/04 Seg |Testes FBST - parte 2/2 e revisão/dúvidas para prova | | | | |
| 27/04 Qua |1a prova | | | | | | 27/04 Qua |1a prova | | | | | ||
+ | | 02/05 Seg |Discussão da 1a prova | | | | | ||
+ | | 04/05 Qua |Atividades dos alunos - revisão da prova | | | | | ||
+ | | 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 === | === 29/02 === | ||
Linha 52: | Linha 58: | ||
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? | 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? | ||
+ | === 09/05 === | ||
+ | - Obter os resultados analíticos possíveis para o problema da questão 5 da prova (posteriori, constante de integração, aproximação quadrática, etc) | ||
+ | - Implementar os diferentes métodos para inferência baseada na posteriori (exata, aproximação normal, discretização, amostragem) | ||
+ | |||
+ | === 11/05 === | ||
+ | - 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'), | ||
+ | 1000) | ||
+ | class(sam) | ||
+ | |||
+ | ## ... ou via coda | ||
+ | sam <- coda.samples(jags, | ||
+ | c('b0', 'b1', 'sigma'), | ||
+ | 1000) | ||
+ | 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"), | ||
+ | 1000) | ||
+ | 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. | ||
+ | - Identificar e ajustar modelos (não bayesianos, bayesianos por simulação ou aproximados) para dados simulados da seguinte forma: <code R> | ||
+ | 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] | ||
+ | </code> | ||
+ |