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 | ||
disciplinas:ce227-2016-01:historico [2016/04/18 19:22] paulojus |
disciplinas:ce227-2016-01:historico [2016/06/02 14:47] (atual) paulojus |
||
---|---|---|---|
Linha 28: | Linha 28: | ||
| 06/04 Qua |desenvolver análises análogas às vistas na última aula para algum outro modelo com 1 parâmetro (excluindo da binomial ou algum dos parâmetros da normal) | | | | | | 06/04 Qua |desenvolver análises análogas às vistas na última aula para algum outro modelo com 1 parâmetro (excluindo da binomial ou algum dos parâmetros da normal) | | | | | ||
| 11/04 Seg |Discussão das análises feitas pelos participantes do curso. Modelos com mais de um parâmetro - ideais fundamentais. Distribuições posterioris marginais, conjuntas e condicionais. |Cap 4 do material do curso | | | | | 11/04 Seg |Discussão das análises feitas pelos participantes do curso. Modelos com mais de um parâmetro - ideais fundamentais. Distribuições posterioris marginais, conjuntas e condicionais. |Cap 4 do material do curso | | | | ||
- | | 13/04 Qua |Resumos da posteriori |Cap 5 do material do curso | | | | + | | 13/04 Qua |Resumos da posteriori |Cap 5 do material do curso | |Preparar material para discussão sobre FBST | |
| 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 | | | | | ||
+ | | 25/04 Seg |Testes FBST - parte 2/2 e revisão/dúvidas para 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]] | | ||
+ | | 18/05 Qua |Inferência Bayesiana utilizando o JAGS/INLA - mais exemplos | | |[[#18/05|Ver abaixo]] | | ||
+ | | 23/05 Seg |Estudos (prof. em congresso) | | | | | ||
+ | | 25/05 Qua |Estudos (prof. em congresso) | | | | | ||
+ | | 31/05 Seg |Aplicação de inferência Bayesiana - erros e incertezas em estimação de vazão de uma bacia - Apres. Alana | | | | | ||
+ | | 01/06 Qua |Fundamentos do INLA | | |[[#01/06|Ver abaixo]] | | ||
+ | |||
=== 29/02 === | === 29/02 === | ||
Linha 49: | Linha 63: | ||
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> | ||
+ | |||
+ | <fs large>**Atividades propostas:**</fs> | ||
+ | - 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: <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> | ||
+ | |||
+ | === 18/05 === | ||
+ | - {{:disciplinas:ce227:changepointjags.r|Script R/JAGS para análise dos dados do Cap 8}} (changepoint Poisson) | ||
+ | - {{:disciplinas:ce227:ce227-av04.pdf|Arquivo com exemplos de modelos visto em aula}} | ||
+ | - Mais um exemplo de análise com efeitos aleatórios (serialmente) correlacionados<code R> | ||
+ | ## | ||
+ | ## Análise de conjunto de dados com INLA com efeitos aleatórios temporalmente correlacionados | ||
+ | ## | ||
+ | require(INLA) | ||
+ | ## | ||
+ | ## Visualização dos 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) | ||
+ | ## ou então, como neste modelo todos os valores preditos são iguais bastaria fazer: | ||
+ | abline(h=fitted(fit.glm)[1], col=2, lty=3, lwd=3) | ||
+ | ## | ||
+ | ## 2. Agora o mesmo modelo nulo anterior porém ajustado 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)) | ||
+ | ## | ||
+ | ## 3. Modelo com probabilidades variando no tempo | ||
+ | ## através da inclusão de variável/processo latente | ||
+ | ## modelando o logito(probabilidade) como um efeito aleatório correlacionado no tempo | ||
+ | ## segundo um "random walk" cíclico de ordem 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édia são praticamente indistinguíveis) | ||
+ | with(fit, matlines(summary.fitted.values[,c(1,3:6)], lty=c(1,2,2,2,3), col=1)) | ||
+ | ## | ||
+ | ## 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> | ||
+ | |||
+ | === 01/06 === | ||
+ | - **INLA** | ||
+ | - {{:disciplinas:ce227:apresentacao-inla.pdf|INLA - idéias básicas}} | ||
+ | - [[https://www.math.ntnu.no/~ingelins/INLAmai09/Pres/talkHavard.pdf|Apresentação de H. Rue]] | ||
+ | - [[http://www.statistica.it/gianluca/Talks/INLA.pdf|Apresentação de Gianluca Baio]] | ||
+ | - [[https://arxiv.org/pdf/1604.00860v1.pdf|Artigo recente de revisão do INLA]] pelos seus desenvolvedores | ||