Unidade IV - Rjags & MCMCpack
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, ou ler direto do site
bombas=read.table("https://raw.githubusercontent.com/liamorita/estatisticabayesiana/master/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.69673 0.27029 0.0030219 0.0035031
## b 0.92565 0.53391 0.0059693 0.0067438
## lambda[1] 0.05955 0.02556 0.0002857 0.0002904
## lambda[2] 0.10184 0.07847 0.0008774 0.0008923
## lambda[3] 0.08905 0.03701 0.0004138 0.0004247
## lambda[4] 0.11529 0.03015 0.0003371 0.0003371
## lambda[5] 0.60291 0.31811 0.0035566 0.0034810
## lambda[6] 0.60843 0.13656 0.0015268 0.0015268
## lambda[7] 0.90854 0.75177 0.0084050 0.0084969
## lambda[8] 0.89770 0.71957 0.0080450 0.0080448
## lambda[9] 1.59475 0.76269 0.0085272 0.0085276
## lambda[10] 1.98082 0.42006 0.0046964 0.0047495
## x_hat[1] 5.58375 3.39200 0.0379237 0.0379238
## x_hat[2] 1.59350 1.79009 0.0200138 0.0203086
## x_hat[3] 5.60912 3.33099 0.0372416 0.0372997
## x_hat[4] 14.52937 5.38786 0.0602381 0.0602396
## x_hat[5] 3.18012 2.47006 0.0276161 0.0282577
## x_hat[6] 19.07888 6.07542 0.0679253 0.0679267
## x_hat[7] 0.92675 1.23342 0.0137901 0.0132700
## x_hat[8] 0.95900 1.25290 0.0140079 0.0140085
## x_hat[9] 3.36512 2.43255 0.0271967 0.0271984
## x_hat[10] 20.75425 6.38090 0.0713407 0.0713434
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## a 0.28828 0.50130 0.65471 0.85103 1.3251
## b 0.19138 0.54135 0.83068 1.20029 2.2488
## lambda[1] 0.02045 0.04110 0.05592 0.07393 0.1199
## lambda[2] 0.00780 0.04447 0.08249 0.13918 0.2987
## lambda[3] 0.03168 0.06224 0.08398 0.10999 0.1756
## lambda[4] 0.06435 0.09381 0.11255 0.13425 0.1814
## lambda[5] 0.14975 0.36984 0.54623 0.77661 1.3589
## lambda[6] 0.37226 0.51224 0.59751 0.69216 0.9019
## lambda[7] 0.07303 0.37840 0.71112 1.22543 2.8919
## lambda[8] 0.07327 0.38147 0.71585 1.20647 2.7588
## lambda[9] 0.48843 1.03676 1.47179 2.02104 3.3793
## lambda[10] 1.24013 1.68404 1.94832 2.24824 2.8863
## x_hat[1] 1.00000 3.00000 5.00000 8.00000 14.0000
## x_hat[2] 0.00000 0.00000 1.00000 2.00000 6.0000
## x_hat[3] 1.00000 3.00000 5.00000 7.00000 13.0000
## x_hat[4] 5.00000 11.00000 14.00000 18.00000 27.0000
## x_hat[5] 0.00000 1.00000 3.00000 4.00000 9.0000
## x_hat[6] 9.00000 15.00000 19.00000 23.00000 32.0000
## x_hat[7] 0.00000 0.00000 1.00000 1.00000 4.0000
## x_hat[8] 0.00000 0.00000 1.00000 1.00000 4.0000
## x_hat[9] 0.00000 2.00000 3.00000 5.00000 9.0000
## x_hat[10] 10.00000 16.00000 20.00000 25.00000 35.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.1249 Min. :0.02186 Min. :0.006367 Min. :0.0001776
## 1st Qu.:0.5013 1st Qu.:0.54135 1st Qu.:0.041099 1st Qu.:0.0444707
## Median :0.6547 Median :0.83068 Median :0.055923 Median :0.0824929
## Mean :0.6967 Mean :0.92565 Mean :0.059554 Mean :0.1018409
## 3rd Qu.:0.8510 3rd Qu.:1.20029 3rd Qu.:0.073927 3rd Qu.:0.1391801
## Max. :2.3291 Max. :4.31260 Max. :0.222945 Max. :0.6331831
## lambda[3] lambda[4] lambda[5] lambda[6]
## Min. :0.00520 Min. :0.03875 Min. :0.02635 Min. :0.2334
## 1st Qu.:0.06224 1st Qu.:0.09381 1st Qu.:0.36984 1st Qu.:0.5122
## Median :0.08398 Median :0.11255 Median :0.54623 Median :0.5975
## Mean :0.08905 Mean :0.11529 Mean :0.60291 Mean :0.6084
## 3rd Qu.:0.10999 3rd Qu.:0.13425 3rd Qu.:0.77661 3rd Qu.:0.6922
## Max. :0.29548 Max. :0.27471 Max. :2.67005 Max. :1.2863
## lambda[7] lambda[8] lambda[9] lambda[10]
## Min. :0.000462 Min. :0.003306 Min. :0.1036 Min. :0.7371
## 1st Qu.:0.378404 1st Qu.:0.381471 1st Qu.:1.0368 1st Qu.:1.6840
## Median :0.711124 Median :0.715852 Median :1.4718 Median :1.9483
## Mean :0.908545 Mean :0.897705 Mean :1.5948 Mean :1.9808
## 3rd Qu.:1.225430 3rd Qu.:1.206469 3rd Qu.:2.0210 3rd Qu.:2.2482
## Max. :7.411214 Max. :5.906046 Max. :5.8221 Max. :4.7173
## 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.584 Mean : 1.593 Mean : 5.609 Mean :14.53
## 3rd Qu.: 8.000 3rd Qu.: 2.000 3rd Qu.: 7.000 3rd Qu.:18.00
## Max. :24.000 Max. :13.000 Max. :22.000 Max. :40.00
## x_hat[5] x_hat[6] x_hat[7] x_hat[8]
## Min. : 0.00 Min. : 3.00 Min. : 0.0000 Min. : 0.000
## 1st Qu.: 1.00 1st Qu.:15.00 1st Qu.: 0.0000 1st Qu.: 0.000
## Median : 3.00 Median :19.00 Median : 1.0000 Median : 1.000
## Mean : 3.18 Mean :19.08 Mean : 0.9267 Mean : 0.959
## 3rd Qu.: 4.00 3rd Qu.:23.00 3rd Qu.: 1.0000 3rd Qu.: 1.000
## Max. :20.00 Max. :55.00 Max. :11.0000 Max. :10.000
## x_hat[9] x_hat[10]
## Min. : 0.000 Min. : 2.00
## 1st Qu.: 2.000 1st Qu.:16.00
## Median : 3.000 Median :20.00
## Mean : 3.365 Mean :20.75
## 3rd Qu.: 5.000 3rd Qu.:25.00
## Max. :23.000 Max. :58.00summary(samps2[,1])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1249 0.5013 0.6547 0.6967 0.8510 2.3291Diagnó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 1.00
## b 1 1.00
## lambda[1] 1 1.00
## lambda[2] 1 1.00
## lambda[3] 1 1.00
## lambda[4] 1 1.00
## lambda[5] 1 1.01
## lambda[6] 1 1.00
## lambda[7] 1 1.00
## lambda[8] 1 1.00
## lambda[9] 1 1.00
## lambda[10] 1 1.00
## x_hat[1] 1 1.00
## x_hat[2] 1 1.00
## x_hat[3] 1 1.00
## x_hat[4] 1 1.00
## x_hat[5] 1 1.01
## x_hat[6] 1 1.00
## x_hat[7] 1 1.00
## x_hat[8] 1 1.00
## x_hat[9] 1 1.00
## x_hat[10] 1 1.00
##
## Multivariate psrf
##
## 1gelman.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.