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 |
| - | disciplinas:ce227-2018-01:historico [2018/02/26 21:12] paulojus |
|---|---|
| + | Revisão anterior |
| + | disciplinas:ce227-2018-01:historico [2019/11/04 11:52] (atual) paulojus |
| @@ Linha -1,3 +1,3 @@ | |
| ====== | |
| No quadro abaixo será anotado o conteúdo dado em cada aula do curso. \\ | |
| @@ Linha -17,6 +17,38 @@ | |
| | 19/02 Seg |Informações sobre o curso. Motivação inicial pela discussão de 3 situações |Cap 1 da "apostila" do curso |[[#19/02|Ver abaixo e sugerido em aula]] | | | |
| | 21/02 Qua |Discussão dos problemas propostos na aula anterior. Teorema de Bayes para o casos discretos e contínuos. Exemplo verossimilhança binomial e priori Beta. Obtenção de priori a partir da "opinião" sobre a proporção. |Cap 2 da "apostila" do curso, até o exemplo 2.2 |[[#21/02|Ver abaixo e sugerido em aula]] | | |
| | 26/02 Seg |Revisão de conceitos e exemplos discutidos até aqui. Discussão sobre as versões discreta e continua do Teorema de Bayes. Relações entre | |
| | 28/02 Qua |1a avaliação periódica. Discussão da avaliação. Características da abordagem bayesiana: atualização sequencial, suficiência, princípio da verossimilhança |Cap. 2 (todo) |[[http://www.leg.ufpr.br/~paulojus/CE227/ce227-av01-2018.pdf|1a avaliação periódica]] | | |
| | 05/03 Seg |Elicitação de priori para o exemplo de verossimilhança binomial com priori beta - apresentação de implementação. Distribuição marginal e relações com modelos de efeitos aleatórios. Especificação de prioris: conjugadas, impróprias e representações de ignorância |Cap. 3 até 3.4 | | | |
| | 07/03 Qua |Discussão sobre o Capitulo 3 do material |Cap. 3 |[[#07/03|Ver abaixo e sugerido em aula]] | | |
| | 12/03 Seg |Implementação de problemas de aula e materiais (Alcides, Bruna e Hektor). Apresentação inicial de exercício de inferência sobre a variância. | |{{ :disciplinas:ce227:ce227-av01.pdf |Exercício para estudo e discussão }} | | |
| | 14/03 Qua |Inferência para posterioris de forma desconhecida: aproximação normal, discretização e amostragem por MCMC. | |[[#14/03|ver abaixo]] | | |
| | 19/03 Seg |Revisão e dúvidas sobe materia da aula anterior. Implementação da discretização. Amostragem e "jitter" das amostras da discreta | |[[#14/03|ver abaixo]] | | |
| | 21/03 Qua |Não haverá aula expositiva |Revisar conteúdos até aqui. | | | |
| | 26/03 Seg |Resumindo posteriori: Decisão (espaço discreto), estimação pontual (espaço contínuo), intervalos e testes |Cap 5, sec 5.1 |[[#26/03|ver abaixo]] | | |
| | 28/03 Qua |2a avaliação | | | | |
| | 02/04 Seg |Discussão da 2a avaliação | | | | |
| | 04/04 Qua |Resumindo posteriori: estimação pontual (espaço contínuo), intervalos e testes, Predição Bayesiana |Cap 5 e Cap 6 |[[#04/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]]) | | | | |
| | 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 | | | | |
| | 04/06 Seg |Reorganização do curso. Dúvidas e perguntas dos alunos | | | | |
| | 06/06 Qua |2a prova. Toda matéria discutida até aqui | | | | |
| | 11/06 Qua |Bayesiano empírico - modelo Poisson-Gamma para taxas | | | | |
| | 13/06 Qua |Apresentações | | | | |
| | 18/06 Qua |Apresentações | | | | |
| | 20/06 Qua |3a prova. | | | | |
| === 19/02 === | |
| @@ Linha -36,10 +68,438 @@ | |
| === 21/02 === | |
| - Revisitar **Problema 2:** acima | |
| - Considerar o exemplo de aula com dados n=200, y=75 e as prioris escolhidas. Obter gráficos das prioris e posterioris. Verificar o efeito do tamanho de amostra aumentando n e y na mesma proporção e repetindo os gráficos. | |
| - Escrever algum código para obtenção a priori a partir da opinião sobre a proporção cf discutido em aula. Verificar seu código com a sua opinião e as ilustradas em aula (ver tabela a seguir) | |
| ^ Estimativa ^ Intervalo ^ Probabilidade ^ | |
| | 0,63 | (0,40 ; 0,75) | 90% | | |
| | 0,42 | (0,32 ; 0,52) | 80% | | |
| | 0,20 | (0,05 ; 0,35) | 80% | | |
| | 0,50 | (0 ; 1) | 100% | | |
| | 0,30 | (0,20 ; 0,40) | 50% | | |
| === 26/02 === | |
| - Completar problemas propostas nas aulas anteriores após as discussões em aula | |
| - Escrever um código para o Exemplo da Poisson (2.3 do material), que permita | |
| - Ler e resolver exercícios do Capítulo 2 da apostila | |
| === 07/03 === | |
| - Exercícios do Cap 3 | |
| - Escrever um código que receba: modelo, dados, priori conjugada e retorne posteriori (parâmetros e gráficos) | |
| === 14/03 === | |
| - Reler o material do {{ :disciplinas:ce227:ce227-av01.pdf |Exercício para estudo e discussão }} | |
| - Montar um algoritmo para aproximação por discretização do exemplo | |
| - Tome algum outro modelo de um parâmetro e desenvolva os resultados análogos aos vistos em aula. | |
| === 19/03 === | |
| - Visualizar, experimentar e comentar o [[https://brunaw.shinyapps.io/bayes_app/|aplicativo shiny]] construído por Bruna, Hektor e Alcides | |
| === 26/03 === | |
| - Refazer exemplos e fazer Exercício 5.1 do Cap 5 | |
| - Código para o Exemplo 5.1: <code R> | |
| # Priori [\theta] | |
| (th <- c(th1=0.6, th2=0.4)) | |
| # Verossimilhança [Y|\theta] | |
| y.th1 <- c(y1=0.35, y2=0.30, y3=0.21, y4=0.14) | |
| y.th2 <- c(y1=0.09, y2=0.17, y3= 0.25, y4=0.49) | |
| (y.th <- rbind(y.th1, y.th2)) | |
| # Conjunta [Y, \theta] = [\theta] \cdot [Y|\theta] | |
| (yth <- th * y.th) | |
| rownames(yth) <- c("yth1", "yth2") | |
| yth | |
| # Marginal [Y] | |
| (y <- drop(crossprod(th, y.th))) | |
| # ou ... | |
| colSums(th * y.th) | |
| # Posteriori | |
| (th.y <- t(t(yth)/drop(y))) | |
| rownames(th.y) <- c("th1.y", "th2.y") | |
| th.y | |
| ## Fc Perda | |
| L <- diag(c(8,20)) | |
| rownames(L) <- paste("a", 1:2, sep="") | |
| L | |
| (L.th.y <- L %*% th.y) | |
| ## Decisão | |
| D.f <- function(x) ifelse(x[1] < x[2], "a1:Vacina", "a2:Não Vacina") | |
| apply(L.th.y, 2, D.f) | |
| ## Perda baseada na regra de decisão para cada resultado de exame | |
| apply(L.th.y, 2, min) | |
| ## Risco de Bayes (associado a uma determinada regra adotada aqui) | |
| ## -- perda média esperada | |
| sum(apply(L.th.y, 2, min) * y) | |
| ## Outra regra: Vacina todo mundo | |
| L <- diag(c(8,0)) | |
| (LT.th.y <- LT %*% th.y) | |
| sum(apply(LT.th.y, 2, sum) * y) | |
| ## ou simplesmente... | |
| sum(th * c(8,0)) | |
| ## Ainda outra regra: não vacina ninguém | |
| LN <- diag(c(0,20)) | |
| (LN.th.y <- LN %*% th.y) | |
| sum(apply(LN.th.y, 2, sum) * y) | |
| ## ou simplesmente... | |
| sum(th * c(0,20)) | |
| </code> | |
| === 02/04 === | |
| - Na questão 1 verificar como a mudança na priori (proporção de motoristas acima do limite) afeta os resultados | |
| - Na questão 2 refazer com a parametrização alternativa da exponencial e gama | |
| - Ainda na questão 2 fazer utilizando o resultado genérico de prioris conjugadas para família exponencial | |
| - Na questão 4 supor uma amostra de 5 valores (7, 3, 4, 5, 2), obter a posteriori e fazer os gráficos de priori, verossimilhança e posteriori | |
| === 04/04 === | |
| - Refazer exemplos e fazer Exercício 5.2 a 5.4 do Cap 5 | |
| - Refazer exemplos e fazer Exercícios cap 6 | |
| - Escrever funções mostrando como média, mediana e quartis podem ser calculados a partir de minimização de função perda: | |
| - para um conjunto de dados | |
| - para uma distribuição discreta | |
| - para uma distribuição contínua | |
| === 09/04 === | |
| - Fazer um código (com operações matriciais) para os cálculos do Exemplo 1. O código deve permitir definir diferentes prioris e verossimilhanças. Experimentar com valores diferentes do exemplo. | |
| - Especificar valores para os hiperparâmetros //p// e //q// no Exemplo 2 e simular um conjunto de dados. Obter a posteriori e maginais. Fazer gráficos conjuntos e marginais da priori e posteriori. | |
| - No Exemplo 3 obter a marginal <latex>[\sigma^2|y]</latex> e a posteriori condicional <latex>[\mu|\sigma^2,y]</latex> | |
| - Ainda no exemplo 3 definir os hiperparâmetros de obter uma simulação de dados do modelo | |
| - Com os dados simulados obter as expressões da posteriori conjunta, marginais (do material) e condicional conjunta (item anterior) | |
| - 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 provou-se 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)</math>. (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> | |
