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-2018-01:historico [2018/04/10 13:47] paulojus |
disciplinas:ce227-2018-01:historico [2018/05/28 15:44] paulojus |
||
|---|---|---|---|
| Linha 31: | Linha 31: | ||
| | 09/04 Seg |Inferência em problemas com mais de um parâmetro |Cap 4: ler, estudar e refazer exemplos |[[#09/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]]) | | | | | 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 | | | | ||
| + | |||
| === 19/02 === | === 19/02 === | ||
| Linha 156: | Linha 171: | ||
| - Obter uma simulação da posteriori. Comparar a conjunta e marginais teórica e simulada. | - 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<code R> | ||
| + | ## 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) | ||
| + | </code> | ||
| + | |||
| + | === 18/04 === | ||
| + | Código visto em aula<code R> | ||
| + | ## | ||
| + | ## 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 provouse 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) | ||
| + | </code> | ||
| + | |||
| + | === 23/04 === | ||
| + | - Implementar modelo semelhante ao visto em aula porém com <math>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:<code R> | ||
| + | ## 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) | ||
| + | </code> | ||
| + | |||
| + | === 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 <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> | ||
