Notas de aula de Estatística Bayesiana
Lia Hanna Martins Morita
2020
Aplicações com o pacote rjags, coda & MCMCpack
require(rjags)
require(coda)
require(DT)
require(tidyverse)
require(car)
#install.packages("MCMCpack")
require(MCMCpack)
#install.packages("bayesplot")
require(bayesplot)Exemplo 1:
Utilização de um modelo hierárquico do tipo Gamma-Poisson no exemplo das bombas de usina referência na página 20 .
Temos 10 bombas de usina de força e a variável número de falhas por bomba (X) é dada por um modelo de Poisson: \(X | \lambda_i \sim Poisson (\lambda_i t_i), i=1, \ldots,10\) onde \(\lambda_i\) é a taxa de falha por bomba e \(t_i\) é o tempo de operação da bomba (em 1000 segundos). Abaixo estão os dados observados:
Copiar os dados em um txt no seu computador
bombas=read.table("exemplo1_bombas.txt",sep=";",head=TRUE)
datatable(bombas)- Deseja-se estimar a taxa média de falhas para cada bomba. Assim, como a conjugada da Poisson é a Distribuição Gamma, assumiu-se uma Distribuição Gamma para \(\lambda_i\):
\[ \lambda_i \sim \text{Gamma}(a,b), \text{ onde a e b são os hiperparâmetros:}\\ a \sim \text{Exponencial}(1)\\ b \sim \text{Gamma}(0.1,1) \]
Entrada dos dados e valores iniciais
dados=list("n"=nrow(bombas),"t"=bombas$t,"x"=bombas$x)
iniciais=list(
list(a=1,b=1),
list(a=2,b=2)
)Sintaxe e inicialização do modelo
cat("
model{
for (i in 1 : n) {
lambda[i] ~ dgamma(a, b)
x[i] ~ dpois(lambda[i]*t[i])
x_hat[i] ~ dpois(lambda[i]*t[i])
}
a ~ dexp(1)
b ~ dgamma(0.1, 1.0)
}
",file="modelo3.txt")
modelo=jags.model(file="modelo3.txt",data=dados,inits=iniciais,n.chains=2, n.adapt=1000)## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 10
## Unobserved stochastic nodes: 22
## Total graph size: 55
##
## Initializing modelEspecifica as quantidades a serem monitorados
algumas vêm direto dos parâmetros que foram amostrados, e outras são funções dos parâmetros
- Implementa o modelo e armazena as amostras da distribuição a posteriori
params = c("a","b","lambda","x_hat")
samps = coda.samples( modelo, params, n.iter=20000,thin=5)
summary(samps)##
## Iterations = 1005:21000
## Thinning interval = 5
## Number of chains = 2
## Sample size per chain = 4000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## a 0.69614 0.27400 0.0030635 0.0035668
## b 0.92223 0.54850 0.0061324 0.0070102
## lambda[1] 0.05957 0.02505 0.0002800 0.0002820
## lambda[2] 0.10300 0.08000 0.0008944 0.0009070
## lambda[3] 0.08984 0.03798 0.0004247 0.0004246
## lambda[4] 0.11583 0.03032 0.0003389 0.0003435
## lambda[5] 0.60296 0.31815 0.0035571 0.0036349
## lambda[6] 0.60663 0.13589 0.0015193 0.0015193
## lambda[7] 0.89008 0.72091 0.0080601 0.0078042
## lambda[8] 0.89365 0.72351 0.0080891 0.0084512
## lambda[9] 1.58794 0.77606 0.0086766 0.0086766
## lambda[10] 1.99368 0.41967 0.0046921 0.0046923
## x_hat[1] 5.63400 3.34387 0.0373856 0.0373862
## x_hat[2] 1.64637 1.82537 0.0204083 0.0204353
## x_hat[3] 5.63112 3.38275 0.0378203 0.0388128
## x_hat[4] 14.58900 5.41357 0.0605256 0.0598271
## x_hat[5] 3.11650 2.42449 0.0271066 0.0271070
## x_hat[6] 19.00600 6.12568 0.0684871 0.0684894
## x_hat[7] 0.92000 1.22231 0.0136658 0.0131729
## x_hat[8] 0.95163 1.23226 0.0137771 0.0137737
## x_hat[9] 3.31400 2.47945 0.0277211 0.0277222
## x_hat[10] 20.91112 6.25263 0.0699065 0.0699096
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## a 0.286898 0.49766 0.65423 0.84662 1.3351
## b 0.191624 0.52840 0.81518 1.19162 2.2635
## lambda[1] 0.021323 0.04121 0.05625 0.07388 0.1182
## lambda[2] 0.008709 0.04386 0.08317 0.14025 0.3086
## lambda[3] 0.031114 0.06210 0.08477 0.11159 0.1796
## lambda[4] 0.064075 0.09405 0.11315 0.13474 0.1814
## lambda[5] 0.148797 0.36661 0.54591 0.78208 1.3631
## lambda[6] 0.369026 0.50934 0.59833 0.69390 0.8947
## lambda[7] 0.068607 0.37901 0.70216 1.20184 2.7984
## lambda[8] 0.071766 0.37443 0.70953 1.21241 2.7634
## lambda[9] 0.488402 1.03096 1.45441 1.99195 3.4787
## lambda[10] 1.257767 1.69542 1.96024 2.26072 2.8909
## x_hat[1] 1.000000 3.00000 5.00000 8.00000 13.0250
## x_hat[2] 0.000000 0.00000 1.00000 2.00000 6.0000
## x_hat[3] 1.000000 3.00000 5.00000 8.00000 14.0000
## x_hat[4] 6.000000 11.00000 14.00000 18.00000 26.0000
## x_hat[5] 0.000000 1.00000 3.00000 4.00000 9.0000
## x_hat[6] 8.000000 15.00000 18.00000 23.00000 32.0000
## x_hat[7] 0.000000 0.00000 1.00000 1.00000 4.0000
## x_hat[8] 0.000000 0.00000 1.00000 1.00000 4.0000
## x_hat[9] 0.000000 2.00000 3.00000 5.00000 9.0000
## x_hat[10] 10.000000 17.00000 21.00000 25.00000 34.00001 Outro modo de armazenar
utilizando a função as.mcmc.list do pacote coda
- desta forma o parâmetro “a” fica guardado no objeto samps2[[1]] da lista samps2
samps2=rbind(as.mcmc.list(samps)[[1]],as.mcmc.list(samps)[[2]])
summary(samps2)## a b lambda[1] lambda[2]
## Min. :0.1006 Min. :0.009067 Min. :0.007352 Min. :0.000222
## 1st Qu.:0.4977 1st Qu.:0.528396 1st Qu.:0.041211 1st Qu.:0.043862
## Median :0.6542 Median :0.815181 Median :0.056246 Median :0.083174
## Mean :0.6961 Mean :0.922231 Mean :0.059573 Mean :0.103001
## 3rd Qu.:0.8466 3rd Qu.:1.191616 3rd Qu.:0.073881 3rd Qu.:0.140246
## Max. :2.3333 Max. :4.598982 Max. :0.180295 Max. :0.640348
## lambda[3] lambda[4] lambda[5] lambda[6]
## Min. :0.01069 Min. :0.04079 Min. :0.0324 Min. :0.2217
## 1st Qu.:0.06210 1st Qu.:0.09405 1st Qu.:0.3666 1st Qu.:0.5093
## Median :0.08477 Median :0.11315 Median :0.5459 Median :0.5983
## Mean :0.08984 Mean :0.11583 Mean :0.6030 Mean :0.6066
## 3rd Qu.:0.11159 3rd Qu.:0.13474 3rd Qu.:0.7821 3rd Qu.:0.6939
## Max. :0.28711 Max. :0.25396 Max. :2.4768 Max. :1.2323
## lambda[7] lambda[8] lambda[9] lambda[10]
## Min. :0.000297 Min. :0.002948 Min. :0.08235 Min. :0.7608
## 1st Qu.:0.379015 1st Qu.:0.374433 1st Qu.:1.03096 1st Qu.:1.6954
## Median :0.702160 Median :0.709535 Median :1.45441 Median :1.9602
## Mean :0.890077 Mean :0.893653 Mean :1.58794 Mean :1.9937
## 3rd Qu.:1.201835 3rd Qu.:1.212406 3rd Qu.:1.99195 3rd Qu.:2.2607
## Max. :6.854103 Max. :7.225609 Max. :6.73311 Max. :3.7544
## x_hat[1] x_hat[2] x_hat[3] x_hat[4]
## Min. : 0.000 Min. : 0.000 Min. : 0.000 Min. : 1.00
## 1st Qu.: 3.000 1st Qu.: 0.000 1st Qu.: 3.000 1st Qu.:11.00
## Median : 5.000 Median : 1.000 Median : 5.000 Median :14.00
## Mean : 5.634 Mean : 1.646 Mean : 5.631 Mean :14.59
## 3rd Qu.: 8.000 3rd Qu.: 2.000 3rd Qu.: 8.000 3rd Qu.:18.00
## Max. :24.000 Max. :14.000 Max. :22.000 Max. :43.00
## x_hat[5] x_hat[6] x_hat[7] x_hat[8]
## Min. : 0.000 Min. : 3.00 Min. : 0.00 Min. :0.0000
## 1st Qu.: 1.000 1st Qu.:15.00 1st Qu.: 0.00 1st Qu.:0.0000
## Median : 3.000 Median :18.00 Median : 1.00 Median :1.0000
## Mean : 3.116 Mean :19.01 Mean : 0.92 Mean :0.9516
## 3rd Qu.: 4.000 3rd Qu.:23.00 3rd Qu.: 1.00 3rd Qu.:1.0000
## Max. :17.000 Max. :45.00 Max. :10.00 Max. :9.0000
## x_hat[9] x_hat[10]
## Min. : 0.000 Min. : 4.00
## 1st Qu.: 2.000 1st Qu.:17.00
## Median : 3.000 Median :21.00
## Mean : 3.314 Mean :20.91
## 3rd Qu.: 5.000 3rd Qu.:25.00
## Max. :21.000 Max. :52.00summary(samps2[[1]])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.5306 0.5306 0.5306 0.5306 0.5306 0.5306Diagnóstico de Gelman & Rubin: através do pacote coda
- Relaciona as variâncias entre as cadeias e dentro das cadeias,
- O valor da estatística deve ser entre 1.0 e 1.2 para indicar convergência.
- Saiba mais sobre diagnósticos de Gelman & Rubin .
- Saiba mais sobre Métodos MCMC .
- Verificamos que houve convergência das cadeias
Referências na literatura:
[1] Gelman, A and Rubin, DB (1992) Inference from iterative simulation using multiple sequences, Statistical Science, 7, 457-511.
[2] Brooks, SP. and Gelman, A. (1998) General methods for monitoring convergence of iterative simulations. Journal of Computational and Graphical Statistics, 7, 434-455.
gelman.diag(samps, confidence = 0.95, transform=FALSE)## Potential scale reduction factors:
##
## Point est. Upper C.I.
## a 1.00 1.00
## b 1.00 1.00
## lambda[1] 1.00 1.00
## lambda[2] 1.00 1.00
## lambda[3] 1.00 1.00
## lambda[4] 1.00 1.00
## lambda[5] 1.00 1.00
## lambda[6] 1.01 1.01
## lambda[7] 1.00 1.00
## lambda[8] 1.00 1.00
## lambda[9] 1.00 1.01
## lambda[10] 1.00 1.00
## x_hat[1] 1.00 1.00
## x_hat[2] 1.00 1.00
## x_hat[3] 1.00 1.01
## x_hat[4] 1.00 1.00
## x_hat[5] 1.00 1.00
## x_hat[6] 1.00 1.00
## x_hat[7] 1.00 1.01
## x_hat[8] 1.00 1.01
## x_hat[9] 1.00 1.00
## x_hat[10] 1.00 1.00
##
## Multivariate psrf
##
## 1.01gelman.plot(samps)
Traços e densidades através do pacote bayesplot
- Desconsiderar as densidades a posteriori dos valores preditos “x_hat”, pois como estas quantidades são discretas (distribuição de Poisson), os gráficos são inapropriados.
mcmc_trace(samps)
mcmc_dens(samps)
Gráficos de autocorrelação do pacote coda
- os gráficos em preto são referentes a primeira cadeia
- os gráficos em vermelho são referentes a segunda cadeia
autocorr.plot(samps[[1]],col=1)
autocorr.plot(samps[[2]],col=2)
Conclusões: Com a inferência Bayesiana e o pacote Rjags, é possível estabelecer a taxa média de falha para cada bomba. Verificou-se convergência do algoritmo.
Pendências:
- Verificar se existe uma forma de estabelecer semente - igual set.seed(), testei nesta análise com rjags mas não funciona.
- Segundo este fórum precisa “setar” a semente de outra forma https://sourceforge.net/p/mcmc-jags/discussion/610037/thread/c6b2ee2b/
- Seria interessante agregar os sumários a posteriori na tabela dos dados das bombas.