Comparação entre os pacotes de MCMC

Apresentação

Bugs

Bayesian inference Using Gibbs Sampler ou ‘Inferência Bayesina utilizando o amostrador de Gibbs’, dorovante BUGS, é um projeto concebido na criação de um software flexível para análise de modelos estatísticos de alta complexidade, utilizando Markov Chain Monte Carlo como método. O projeto começou em 1989 sob a tutela de Spiegelhalter et al, na MCR Biostatistics Unit, em Camdribge, resultando o programa inicial denominado de BUGS. Furthermore, em um projeto conjunto com a Imperial College, Londres, foi desenvolvido o WinBugs, uma versão compatível com o R.

Posteriormente, foi apresentado ao público o OpenBUgs, uma versão acessível ao público geral e igualmente compatível com o R, assim como a versão anterior. Howerver, a produção de seu conteúdo foi encerrada em 2011, embora haja versões disponíveis até o presente dia. Atualmente, conta com mais de 30 mil downloas de usuários e uma interface gráfica única.

Em 2018, Goudie el at desenvolveram o projeto MultiBugs, que oferece um algoritmo em computação parelala, objetivando a compilação rápida de modelos mais robustos.

Jags

Just Another Gibbs Sampler (JAGS), ou Apenas Outro amostra de Gibbs, é um software ao moldes do GIBS, apresentado anteriormente, criado sob as mesmas premissas e com o mesmo objetivo, Inferência Bayesiana via MCMC. Para além das similaridades mencionadas, o JAGS é mais flexível, se comparado com o GIBS, no que tange a interconectividade com outros softwares estatísticos, em especial, o R.

Desenvolvido por Martyn Plummer, é consideravelmente mais recente que o projeto anterior e sua última atualização data de 2017. Embora, assim como o BUGS, é necessário que o compilador esteja atuando em concomitância com o R, ou seja, não há independência. O JAGS é de uso livre e licensiado sob a GNU General Public License.

Stan

Stan é pacote estatístico com o objetivo de realizar inferência bayesiana sob a ótica do MCMC. Desenvolvido por Gelman, Carpenter e um time de 34 pesquisados, o programa é uma homenagem ao pinoeiro do método de Monte Carlo e físico responsável pela descoberta da fissão atômica, Stanislaw Ulam, um dos membros do projeto Manhattam. É licensiado pela New BDS License, sendo de origem recente, produzido em 2012, conta com atualização e feedback de sua equipe até o presente dia.

Em suma, as diferenças para o JAGS e o BUGS consta na compilação do modelo, os primeiros, mais simples, seguem o mesmo padrão. O Stan, por sua vez, adotou uma linguagem C++ para a escrita da modelagem. O seu diferencial está em seus algoritmos, sendo o Hamiltonian Monte Carlo (HMC) (default) e o No U turner sampler (NUTS), alternativo ao primeiro.

Nimble

O nimble é sistema integrado ao R, em formato de pacote, desenvolvido por uma pesquisa colaborativa entre P de Valpine, C Paciorek e D Temple Lang. Tem como objetivo principal a difusão da análise e cômputo de modelos de alta dimensionalidade, ou computacionalmente intensivos, em especial os modelos hierárquicos. O compilador utiliza-se da linguagem do BUGS, já bem definida pela comunidade, ademias, seus rotinas internas são integradas em C++ para ganhos de velocidade.

Embora seja considerado um compilador User Friendly para MCMC, o Nimble possui outras ferramentas. Há a possiblidade também de escrever modelos Bayesianos Não Paramétricos (BNP) de modelos de mistura e estimação de densidade não paramétrica via processo de mistura Dirichlet.

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                                      #
######################################################################################

#Modelo do tipo BUGS
sink("mod1.txt")        
cat("
    MODEL LR1 {

    #Verossimilhança 
    for(i in 1:N) {
    y[i] ~ dnorm(mu[i], tau)
    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)
    tau ~ dgamma(1,0.1)
    
    #Transformação para variância
    sigma2 <- 1/tau
    } ", 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","tau","sigma2")  

#Valores iniciais
inits <- function () {list(beta0 = rnorm(1),
                           beta1 = rnorm(1),
                           beta2 = rnorm(1),
                           tau = 1)
                           sigma2 = 1}

#Modelo final
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)
toc()

Resultados

	mean	sd	2.5%	25%	50%	75%	97.5%	Rhat	n.eff
beta0	-1.7017525	0.2174842	-2.1550000	-1.84525	-1.70400	-1.554000	-1.279975	1.002923	830
beta1	1.6939550	0.2304004	1.2440000	1.53975	1.69900	1.845000	2.144000	1.001317	2400
beta2	10.0042037	0.0220753	9.9620000	9.98900	10.00000	10.020000	10.050000	1.000688	3000
tau	0.2250701	0.0327324	0.1658975	0.20160	0.22385	0.246025	0.293710	1.000817	3000
sigma2	4.5387863	0.6724222	3.4048749	4.06400	4.46750	4.959250	6.026150	1.000817	3000

JAGS

Código

######################################################################################
#                                          JAGS                                      #
######################################################################################
model_code = 'model{
    #Verossimilhança 
    for(i in 1:N) {
    y[i] ~ dnorm(mu[i], tau)
    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)
    tau ~ dgamma(1,0.1)
    
    #Transformação para variância
    sigma2 <- 1/tau}'
  
N = length(y)

model_data = list(N = N, y = y, x1 = x1, x2 = x2)

model_parameters =  c("beta0", "beta1","beta2", "sigma2","tau")

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) 
toc()

Resultados

	mean	sd	2.5%	25%	50%	75%	97.5%	Rhat	n.eff
beta0	-1.6860944	0.2105909	-2.0879681	-1.8296104	-1.687723	-1.5483238	-1.266608	1.000977	3000
beta1	1.6909286	0.2311049	1.2313686	1.5387469	1.688274	1.8493560	2.144819	1.000647	3000
beta2	10.0052898	0.0219353	9.9620366	9.9905397	10.005205	10.0198417	10.048766	1.001766	1600
sigma2	4.5186170	0.6587993	3.4324159	4.0422921	4.452478	4.9056101	5.975496	1.000686	3000
tau	0.2258834	0.0320833	0.1673501	0.2038482	0.224594	0.2473844	0.291340	1.000686	3000

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;
}
model {
alpha ~ normal(0, 10);   
beta1 ~ normal(0, 10); 
beta2 ~ normal(0, 10);
sigma2 ~ inv_gamma(3,20);
y ~ normal(alpha + beta1 * x1 + beta2*x2, sqrt(sigma2));
}

generated quantities {
vector[N] y_rep;
for(n in 1:N) {
y_rep[n] = normal_rng(alpha + beta1 * x1[n] + beta2*x2[n], sqrt(sigma2));
}
}"


data <- list(y = y, x1 = x1, x2 = x2 , N = 100)

tic()
fit <- stan(model_code = model,
            data = data, iter = 10000, warmup = 1000, chains = 3)

stan = extract(fit)
toc()

Resultados

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat alpha -1.705771 0.0013234 0.2187635 -2.138415 -1.851237 -1.704914 -1.559763 -1.276514 27325.03 0.9999686 beta1 1.699761 0.0014114 0.2395953 1.228784 1.540529 1.698885 1.861845 2.168289 28818.32 0.9999350 beta2 10.005224 0.0001332 0.0225257 9.961374 9.990201 10.005070 10.020402 10.049088 28582.34 0.9999602 sigma2 4.745150 0.0041508 0.6696007 3.614317 4.275004 4.679747 5.153950 6.223051 26024.15 1.0000428

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.704459 0.2182247 -2.134405 -1.845417 -1.702805 -1.560221 -1.272177 -1.708991 0.2205426 -2.143654 -1.856214 -1.707929 -1.561689 -1.281975 -1.703863 0.2175001 -2.132866 -1.852140 -1.704202 -1.557887 -1.279405 beta1 1.698571 0.2366761 1.234933 1.541722 1.696640 1.857513 2.164830 1.701441 0.2444538 1.217817 1.537478 1.701037 1.868434 2.175044 1.699270 0.2375979 1.229826 1.541493 1.699162 1.860382 2.159939 beta2 10.005271 0.0223819 9.961223 9.990564 10.005212 10.020494 10.048512 10.005315 0.0224908 9.961959 9.990279 10.004997 10.020309 10.049125 10.005086 0.0227050 9.960872 9.989856 10.004952 10.020399 10.049404 sigma2 4.739498 0.6709983 3.616892 4.267374 4.665881 5.148719 6.239155 4.746080 0.6703198 3.616933 4.276826 4.680797 5.153938 6.220253 4.749871 0.6675122 3.603574 4.279996 4.692062 5.159586 6.208506

NIMBLE

Código

######################################################################################
#                                        NIMBLE                                     #
######################################################################################

simpleCode1 <- nimbleCode({
  beta0 ~ dnorm(0,0.1)
  beta1 ~ dnorm(0,0.1)
  beta2 ~ dnorm(0,0.1)
  tau ~ dgamma(1, 0.1)
  sigma2 <- 1/tau
  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, tau = 2 ))

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

Resultados

Mean Median St.Dev. 95%CI_low 95%CI_upp beta0 -1.7041573 -1.7083052 0.2091917 -2.1043706 -1.2794896 beta1 1.6974375 1.7005161 0.2272487 1.2548992 2.1237661 beta2 10.0050927 10.0049245 0.0223281 9.9638410 10.0516948 sigma2 4.5113108 4.4523352 0.6241347 3.4605394 5.9630216 tau 0.2258069 0.2246013 0.0304566 0.1677002 0.2889724

Mean Median St.Dev. 95%CI_low 95%CI_upp beta0 -1.6968087 -1.6944266 0.2076834 -2.0947722 -1.2962685 beta1 1.6877346 1.6932318 0.2321149 1.2386992 2.1394175 beta2 10.0043223 10.0043528 0.0223773 9.9611573 10.0475828 sigma2 4.4939257 4.4264768 0.6632215 3.4387254 5.9387300 tau 0.2272029 0.2259133 0.0323686 0.1683862 0.2908054

Mean Median St.Dev. 95%CI_low 95%CI_upp beta0 -1.6928879 -1.6885632 0.2167004 -2.102350 -1.270194 beta1 1.6852116 1.6877019 0.2381134 1.216738 2.150423 beta2 10.0055725 10.0063789 0.0221097 9.962371 10.049552 sigma2 4.5020546 4.4327260 0.6709563 3.398402 6.061527 tau 0.2268949 0.2255948 0.0327657 0.164975 0.294256

Mean Median St.Dev. 95%CI_low 95%CI_upp beta0 -1.6979513 -1.6964828 0.2112099 -2.1021654 -1.2796220 beta1 1.6901279 1.6935095 0.2325170 1.2372259 2.1414415 beta2 10.0049958 10.0049545 0.0222705 9.9622690 10.0492859 sigma2 4.5024303 4.4367283 0.6529133 3.4376043 5.9920623 tau 0.2266349 0.2253913 0.0318745 0.1668875 0.2909003

Comparativo entre os tempos

Pacote utilizado: tictoc

	Tempo em Segundos (Arredondado)
Bugs	6
Jags	2
Stan	8
Nimble	30