Comparação entre os pacotes de MCMC

Pacotes

library("R2WinBUGS")
library("R2jags")
library("rstan")
library("nimble")
library("boa")
library("coda")
library("xtable")
library("kableExtra")
library("tictoc")
library("LaplacesDemon")

Dados

######################################################################################
#                                SIMULAÇÃO DOS DADOS                                 #
######################################################################################
set.seed(123)
#Parâmetros verdadeiros
beta0=-2
beta1=2
beta2=10
sigma2=5
n=100
tau=1/sigma2
#Geração
x1=rnorm(n,0,1)
x2=rnorm(n,2,10)
e=rnorm(n,0, sd = sqrt(sigma2))
y=beta0+beta1*x1+beta2*x2+e

BUGS

Código

######################################################################################
#                                          BUGS                                      #
######################################################################################
sink("mod1.txt")        
cat("
    MODEL LR1 {
    #Verossimilhança 
    for(i in 1:N) {
    y[i] ~ dnorm(mu[i], sigma.theta )
    mu[i] <- beta0 + beta1*x1[i] + beta2*x2[i]
    }
    #Proris
    beta0 ~ dnorm(0,0.1)
    beta1 ~ dnorm(0,0.1)
    beta2 ~ dnorm(0,0.1)
    xi ~ dnorm(0,1)
    eta ~ dgamma(1,0.1)
    # A half-cauchy prior distribution on sigma.theta is induced using
    # a normal prior on xi and an inverse-gamma on tau.eta
    tau.eta <- 1/eta
    sigma.theta <- abs(xi)/sqrt(tau.eta)
    
    # Transformando em desvio
    sd <- 1/sqrt(sigma.theta)
    
     
    } ", fill = TRUE)
sink()
#Obeservação: A dist normal no bugs está definido para a precisão. 
#Tamanho amostral
N = length(y)
#Conjunto de dados observados
data = list("N","y","x1","x2")
#Parâmetros estudados
params = c("beta0", "beta1","beta2","sd")  
#Valores iniciais
inits <- function () {list(beta0 = 0,
                           beta1 = 0,
                           beta2 = 0,
                           sigma.teta = 2)}
tic()
result <- bugs(data = data, inits = inits, parameters.to.save = params,
               model.file = "mod1.txt", n.chains = 3, n.iter = 10000,
               n.burnin = 1000, n.thin=9, bugs.directory = "C:\\Program Files\\WinBUGS14",
               debug = TRUE, save.history = TRUE, DIC = FALSE)

cannot create file 'C:\Program Files\WinBUGS14/System/Rsrc/Registry_Rsave.odc', reason 'Permission denied'cannot open file 'C:\Program Files\WinBUGS14/System/Rsrc/Registry.odc': Permission deniedError in file(con, "wb") : cannot open the connection

toc()

21.9 sec elapsed

Resultados

result$summary %>%
  kable() %>%
  kable_styling()

	mean	sd	2.5%	25%	50%	75%	97.5%	Rhat	n.eff
beta0	-1.694912	0.2130390	-2.114025	-1.83325	-1.6965	-1.549	-1.264000	1.000536	3000
beta1	1.690878	0.2363106	1.235925	1.53100	1.6900	1.849	2.150025	1.002342	1100
beta2	10.004714	0.0224442	9.961000	9.98900	10.0000	10.020	10.050000	1.001220	2800
sd	2.119934	0.1512555	1.855000	2.01100	2.1120	2.214	2.438050	1.002110	1200

bugs = result$sims.list
bugs_b0 = bugs$beta0
bugs_b1 = bugs$beta1
bugs_b2 = bugs$beta2
bugs_sd = bugs$sd
plot(bugs_b0,type = "l",main=expression(beta[0]),ylab="",
     xlab="",cex.main=2)
abline(h=-2,col="red")

plot(bugs_b1,type = "l",main=expression(beta[1]),ylab="",
     xlab="",cex.main=2)
abline(h=2,col="red")

plot(bugs_b2,type = "l",main=expression(beta[2]),ylab="",
     xlab="",cex.main=2)
abline(h=10,col="red")

plot(bugs_sd,type = "l",main=expression(sigma),ylab="",
     xlab="",cex.main=2)
abline(h=sqrt(5),col="red")

densplot(as.mcmc(bugs_b0),main=expression(beta[0]),xlab="",ylab="",cex.main=2)

densplot(as.mcmc(bugs_b1),main=expression(beta[1]),xlab="",ylab="",cex.main=2)

densplot(as.mcmc(bugs_b2),main=expression(beta[2]),xlab="",ylab="",cex.main=2)

densplot(as.mcmc(bugs_sd),main=expression(sigma),xlab="",ylab="",cex.main=2)

JAGS

Código

######################################################################################
#                                          JAGS                                      #
######################################################################################
model_code = 'model{
     #Verossimilhança 
    for(i in 1:N) {
    y[i] ~ dnorm(mu[i], sigma.theta )
    mu[i] <- beta0 + beta1*x1[i] + beta2*x2[i]
    }
    #Proris
    beta0 ~ dnorm(0,0.1)
    beta1 ~ dnorm(0,0.1)
    beta2 ~ dnorm(0,0.1)
    xi ~ dnorm(0,1)
    eta ~ dgamma(1,0.1)
    
    tau.eta <- 1/eta
    sigma.theta <- abs(xi)/sqrt(tau.eta)
    
    }'
  
N = length(y)
model_data = list(N = N, y = y, x1 = x1, x2 = x2)
model_parameters =  c("beta0", "beta1","beta2","sigma.theta")
tic()
model_run = jags(data = model_data,
                 parameters.to.save = model_parameters,
                 model.file=textConnection(model_code),
                 n.chains=3, # Number of different starting positions
                 n.iter=10000, # Number of iterations
                 n.burnin=1000, # Number of iterations to remove at start
                 n.thin=9,# Amount of thinning
                 DIC=FALSE) 

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 100
   Unobserved stochastic nodes: 5
   Total graph size: 613

Initializing model
Deleting model

Error in jags.model(model.file, data = data, inits = init.values, n.chains = n.chains,  : 
  Error in node y[1]
Invalid parent values

Resultados

STAN

Código

######################################################################################
#                                          STAN                                      #
######################################################################################
model = "data {
int<lower=1> N;
vector[N] y;
vector[N] x1;
vector[N] x2;
}
parameters {
real alpha;
real beta1;
real beta2;
real<lower=0> sigma2;
}
transformed parameters{
real<lower=0> sigma;
sigma = fabs(sigma2);
}
model {
alpha ~ normal(0, 10);   
beta1 ~ normal(0, 10); 
beta2 ~ normal(0, 10);
sigma2 ~ cauchy(0,10);
y ~ normal(alpha + beta1 * x1 + beta2*x2, sigma);
}
generated quantities {
vector[N] y_rep;
for(n in 1:N) {
y_rep[n] = normal_rng(alpha + beta1 * x1[n] + beta2*x2[n], sigma);
}
}"
data <- list(y = y, x1 = x1, x2 = x2 , N = 100)
tic()
fit <- stan(model_code = model,
            data = data, iter = 2000, warmup = 1000, chains = 3)

SAMPLING FOR MODEL 'b52d6740126bb2b9be591795cab4219a' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 0.316 seconds (Warm-up)
Chain 1:                0.268 seconds (Sampling)
Chain 1:                0.584 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'b52d6740126bb2b9be591795cab4219a' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 0 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 0.266 seconds (Warm-up)
Chain 2:                0.272 seconds (Sampling)
Chain 2:                0.538 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL 'b52d6740126bb2b9be591795cab4219a' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 0 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 0.499 seconds (Warm-up)
Chain 3:                0.659 seconds (Sampling)
Chain 3:                1.158 seconds (Total)
Chain 3: 

stan = extract(fit)
toc()

62.9 sec elapsed

Resultados

summary(fit,probs = c(0.025, 0.25, 0.50, 0.75, 0.975),pars=c("alpha","beta1","beta2","sigma")) %>%
  kable() %>%
  kable_styling()

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat alpha -1.705149 0.0039007 0.2203186 -2.128741 -1.851513 -1.706788 -1.558252 -1.266535 3190.227 0.9994377 beta1 1.704664 0.0042980 0.2346550 1.239475 1.547721 1.702359 1.869757 2.162556 2980.765 1.0002684 beta2 10.004990 0.0004031 0.0226243 9.962109 9.990043 10.004909 10.020552 10.048529 3150.105 0.9996364 sigma 2.159418 0.0031805 0.1589902 1.874101 2.049011 2.152554 2.260097 2.490462 2498.909 0.9995500

mean.chain:1 sd.chain:1 2.5%.chain:1 25%.chain:1 50%.chain:1 75%.chain:1 97.5%.chain:1 mean.chain:2 sd.chain:2 2.5%.chain:2 25%.chain:2 50%.chain:2 75%.chain:2 97.5%.chain:2 mean.chain:3 sd.chain:3 2.5%.chain:3 25%.chain:3 50%.chain:3 75%.chain:3 97.5%.chain:3 alpha -1.707477 0.2235333 -2.126845 -1.855703 -1.716247 -1.560096 -1.282388 -1.703702 0.2068422 -2.103886 -1.844171 -1.703053 -1.570987 -1.271763 -1.704267 0.2301279 -2.147936 -1.860284 -1.701615 -1.551234 -1.259488 beta1 1.698001 0.2377841 1.214058 1.535801 1.696142 1.866904 2.154469 1.710654 0.2302538 1.266526 1.553476 1.710989 1.860361 2.178951 1.705337 0.2359243 1.231920 1.544955 1.700726 1.875302 2.155878 beta2 10.004423 0.0237017 9.958859 9.988924 10.003503 10.021196 10.049177 10.005320 0.0216338 9.963577 9.991049 10.005570 10.020792 10.047768 10.005227 0.0225016 9.963601 9.989519 10.005487 10.019398 10.049736 sigma 2.157314 0.1589916 1.861853 2.046991 2.152927 2.264420 2.484893 2.156001 0.1518553 1.880637 2.049673 2.154927 2.245384 2.475617 2.164939 0.1658292 1.876279 2.051884 2.150899 2.269062 2.528738

stan_b0 = stan$alpha
stan_b1 = stan$beta1
stan_b2 = stan$beta2
stan_sd = stan$sigma
plot(stan_b0,type = "l",main=expression(beta[0]),ylab="",
     xlab="",cex.main=2)
abline(h=-2,col="red")

plot(stan_b1,type = "l",main=expression(beta[1]),ylab="",
     xlab="",cex.main=2)
abline(h=2,col="red")

plot(stan_b2,type = "l",main=expression(beta[2]),ylab="",
     xlab="",cex.main=2)
abline(h=10,col="red")

plot(stan_sd,type = "l",main=expression(sigma),ylab="",
     xlab="",cex.main=2)
abline(h=sqrt(5),col="red")

plot(density(stan_b0),main=expression(beta[0]),cex.main=2)

plot(density(stan_b1),main=expression(beta[1]),cex.main=2)

plot(density(stan_b2),main=expression(beta[2]),cex.main=2)

plot(density(stan_sd),main=expression(sigma),cex.main=2)

NIMBLE

Código

######################################################################################
#                                        NIMBLE                                     #
######################################################################################
simpleCode1 <- nimbleCode({
  
beta0 ~ dnorm(0,0.1)
beta1 ~ dnorm(0,0.1)
beta2 ~ dnorm(0,0.1)
xi ~ dt(0,0.01,1)
sigma  <- abs(xi)
sigma2 <- pow(xi,2)
tau <- pow(sigma2,-1)
  
for(i in 1:N) {
Ypred[i] <- beta0 + beta1 * x1[i] + beta2*x2[i]
Y[i] ~ dnorm(Ypred[i], tau)
}
  
  
})
simpleModel1 <- nimbleModel(simpleCode1,
                            data = list(Y = y, x1 = x1,x2=x2),
                            constants = list(N = N),
                            inits = list(beta0 = 0, beta1 = 0,beta2=0, xi = 2 ))

defining model...
building model...
setting data and initial values...
running calculate on model (any error reports that follow may simply reflect missing values in model variables) ... 
checking model sizes and dimensions...
model building finished.

tic()
mcmc.out <- nimbleMCMC(code = simpleModel1, 
                       data = list(Y = y, x1 = x1,x2=x2), 
                       inits = list(beta0 = 0, beta1 = 0,beta2=0, xi = 2),
                       monitors=c("beta0","beta1","beta2","sigma","sigma2"),
                       nchains = 3, 
                       niter = 10000,
                       thin = 9,
                       nburnin = 1000,
                       summary = TRUE, 
                       progressBar = TRUE)

compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compiler details.
compilation finished.
runMCMC's handling of nburnin changed in nimble version 0.6-11. Previously, nburnin samples were discarded *post-thinning*.  Now nburnin samples are discarded *pre-thinning*.  The number of samples returned will be floor((niter-nburnin)/thin).
running chain 1...

|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|

running chain 2...

|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|

running chain 3...

|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|

toc()

78.03 sec elapsed

Resultados

mcmc.out$summary %>%
  kable() %>%
  kable_styling()

Mean Median St.Dev. 95%CI_low 95%CI_upp beta0 -1.701853 -1.700324 0.2247996 -2.134344 -1.269477 beta1 1.694019 1.699318 0.2317892 1.237297 2.166930 beta2 10.006603 10.006732 0.0223889 9.962705 10.051267 sigma 2.151196 2.139595 0.1634963 1.852558 2.511840 sigma2 4.654348 4.577867 0.7162964 3.431973 6.309339

Mean Median St.Dev. 95%CI_low 95%CI_upp beta0 -1.690154 -1.683048 0.2156525 -2.118917 -1.284972 beta1 1.686202 1.683455 0.2413500 1.218901 2.165015 beta2 10.004045 10.004224 0.0220070 9.959289 10.049076 sigma 2.149450 2.140115 0.1534882 1.879498 2.459997 sigma2 4.643669 4.580092 0.6692461 3.532516 6.051586

Mean Median St.Dev. 95%CI_low 95%CI_upp beta0 -1.703272 -1.700163 0.2156711 -2.120387 -1.287375 beta1 1.691633 1.686491 0.2365682 1.233774 2.149293 beta2 10.004938 10.005058 0.0229612 9.959744 10.047435 sigma 2.153613 2.146234 0.1546927 1.847456 2.494988 sigma2 4.661956 4.606319 0.6731883 3.413092 6.224967

Mean Median St.Dev. 95%CI_low 95%CI_upp beta0 -1.698426 -1.692678 0.2187562 -2.123892 -1.284737 beta1 1.690618 1.690694 0.2365450 1.232427 2.158214 beta2 10.005195 10.005548 0.0224733 9.959804 10.049202 sigma 2.151420 2.142408 0.1572459 1.856467 2.484581 sigma2 4.653324 4.589912 0.6863865 3.446468 6.173146

nimble_b0  =  c(mcmc.out$samples$chain1[,1],mcmc.out$samples$chain2[,1],mcmc.out$samples$chain3[,1])
nimble_b1  =  c(mcmc.out$samples$chain1[,2],mcmc.out$samples$chain2[,2],mcmc.out$samples$chain3[,2])
nimble_b2  =  c(mcmc.out$samples$chain1[,3],mcmc.out$samples$chain2[,3],mcmc.out$samples$chain3[,3])
nimble_sigma  =  c(mcmc.out$samples$chain1[,4],mcmc.out$samples$chain2[,4],mcmc.out$samples$chain3[,4])
plot(nimble_b0,type = "l",main=expression(beta[0]),ylab="",
     xlab="",cex.main=2)
abline(h=-2,col="red")

plot(nimble_b1,type = "l",main=expression(beta[1]),ylab="",
     xlab="",cex.main=2)
abline(h=2,col="red")

plot(nimble_b2,type = "l",main=expression(beta[2]),ylab="",
     xlab="",cex.main=2)
abline(h=10,col="red")

plot(nimble_sigma,type = "l",main=expression(sigma),ylab="",
     xlab="",cex.main=2)
abline(h=sqrt(5),col="red")

plot(density(nimble_b0),main=expression(beta[0]),cex.main=2)

plot(density(nimble_b1),main=expression(beta[1]),cex.main=2)

plot(density(nimble_b2),main=expression(beta[2]),cex.main=2)

plot(density(nimble_sigma),main=expression(sigma),cex.main=2)