Comparação

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


######################################################################################
#                                SIMULAÇÃO DOS DADOS                                 #
######################################################################################
set.seed(2019)
beta0=-2
beta1=2
beta2=10
sigma2=5
tau=1/sigma2

yt=list()

for (i in 1:100){
n=100
x=rnorm(n,0,1)
e=rnorm(n,0, sd = sqrt(sigma2))  
yt[[i]]=list(beta0+beta1*x+e,x)
}


amostra = list()
final = list()
b0  = list()
b1  = list()
tau = list()
media_b0 = vector()
media_b1 = vector()
media_tau = vector()
mediana_b0 = vector()
mediana_b1 = vector()
mediana_tau = vector()
var_b0 = vector()
var_b1 = vector()
var_tau = vector()
eqm_b0  = vector()
eqm_b1  = vector()
eqm_tau = vector()
hdi = list()

tic()
for (i in 1:100){
y =  yt[[i]][[1]]
x =  yt[[i]][[2]]
sink("mod1.txt")        
cat("
    MODEL LR1 {
    #Verossimilhança 
    for(i in 1:N) {
    y[i] ~ dnorm(mu[i], tau )
    mu[i] <- beta0 + beta1*x[i]
    }
    #Proris
    beta0 ~ dnorm(0,0.1)
    beta1 ~ dnorm(0,0.1)
    tau ~ dgamma(1,0.1)
    } ", fill = TRUE)
sink()

N = 100;data = list("N","y","x");params = c("beta0", "beta1","tau") ;inits <- function () {list(beta0 = 0,
                           beta1 = 0,
                           tau = 2)}

amostra[[i]] <- 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 = FALSE, save.history = FALSE, DIC = FALSE)

final[[i]] = amostra[[i]]$sims.matrix

b0[[i]]  = final[[i]][,1]
b1[[i]]  = final[[i]][,2]
tau[[i]] = final[[i]][,3]

media_b0[i]  = mean(b0[[i]]) 
media_b1[i]  = mean(b1[[i]]) 
media_tau[i] = mean(tau[[i]]) 

mediana_b0[i]  = median(b0[[i]]) 
mediana_b1[i]  = median(b1[[i]]) 
mediana_tau[i] = median(tau[[i]])  

var_b0[i]  = var(b0[[i]])
var_b1[i]  = var(b1[[i]])
var_tau[i] = var(tau[[i]])

eqm_b0[i]  = mean((b0[[i]]+2)^2)
eqm_b1[i]  = mean((b1[[i]]-2)^2)
eqm_tau[i] = mean((tau[[i]]-.2)^2)

hdi[[i]] = hdi(final[[i]])

}
toc()

media=c(mean(media_b0),mean(media_b1),mean(media_tau))
mediana=c(mean(mediana_b0),mean(mediana_b1),mean(mediana_tau))
var=c(mean(var_b0),mean(var_b1),mean(var_tau))
eqm=c(mean(eqm_b0),mean(eqm_b1),mean(eqm_tau))
medidas=cbind(media,mediana,var,eqm)
rownames(medidas) = c("b0","b1","tau");medidas

int = as.data.frame(cbind(t(as.data.frame((hdi))),rep(c(-2,2,0.2),100)))

write.csv2(int,"C:\\Users\\Rafael\\Desktop\\int.xls")



######################################################################################
#                                        JAGS                                        #
######################################################################################



amostra = list()
final = list()
b0  = list()
b1  = list()
tau = list()
media_b0 = vector()
media_b1 = vector()
media_tau = vector()
mediana_b0 = vector()
mediana_b1 = vector()
mediana_tau = vector()
var_b0 = vector()
var_b1 = vector()
var_tau = vector()
eqm_b0  = vector()
eqm_b1  = vector()
eqm_tau = vector()
hdi = list()

tic()
for (i in 1:100){
  y =  yt[[i]][[1]]
  x =  yt[[i]][[2]]
  

  model_code = '
model
{
  # Vero
  for (i in 1:n) {
  y[i] ~ dnorm(beta1 + beta2 * x[i], tau)
  }
  # Prioris
  beta1 ~ dnorm(0, 01)
  beta2 ~ dnorm(0,0.1)
  tau ~ dgamma(1,0.1)
}
'
  
model_data = list(n = n, y = y, x = x)
  
model_parameters =  c("beta1", "beta2", "tau")
  
amostra[[i]] = 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)

final[[i]] = amostra[[i]][["BUGSoutput"]][["sims.matrix"]]

b0[[i]]  = final[[i]][,1]
b1[[i]]  = final[[i]][,2]
tau[[i]] = final[[i]][,3]

media_b0[i]  = mean(b0[[i]]) 
media_b1[i]  = mean(b1[[i]]) 
media_tau[i] = mean(tau[[i]]) 

mediana_b0[i]  = median(b0[[i]]) 
mediana_b1[i]  = median(b1[[i]]) 
mediana_tau[i] = median(tau[[i]])  

var_b0[i]  = var(b0[[i]])
var_b1[i]  = var(b1[[i]])
var_tau[i] = var(tau[[i]])

eqm_b0[i]  = mean((b0[[i]]+2)^2)
eqm_b1[i]  = mean((b1[[i]]-2)^2)
eqm_tau[i] = mean((tau[[i]]-.2)^2)

hdi[[i]] = hdi(final[[i]])

}
toc()

media=c(mean(media_b0),mean(media_b1),mean(media_tau))
mediana=c(mean(mediana_b0),mean(mediana_b1),mean(mediana_tau))
var=c(mean(var_b0),mean(var_b1),mean(var_tau))
eqm=c(mean(eqm_b0),mean(eqm_b1),mean(eqm_tau))
medidas=cbind(media,mediana,var,eqm)
rownames(medidas) = c("b0","b1","tau");medidas

int = as.data.frame(cbind(t(as.data.frame((hdi))),rep(c(-2,2,0.2),100)))

write.csv2(int,"C:\\Users\\Rafael\\Desktop\\int.xls")


######################################################################################
#                                        STAN                                        #
######################################################################################


amostra = list()
final = list()
b0  = list()
b1  = list()
tau = list()
media_b0 = vector()
media_b1 = vector()
media_tau = vector()
mediana_b0 = vector()
mediana_b1 = vector()
mediana_tau = vector()
var_b0 = vector()
var_b1 = vector()
var_tau = vector()
eqm_b0  = vector()
eqm_b1  = vector()
eqm_tau = vector()
hdi = list()

tic()
for (i in 1:100){
  y =  yt[[i]][[1]]
  x =  yt[[i]][[2]]
  
model = "data {
int<lower=1> N;
vector[N] y;
vector[N] x;
}
parameters {
real beta1;
real beta2;
real<lower=0> sigma2;
}
model {
beta1 ~ normal(0, 10);   
beta2 ~ normal(0, 10); 
sigma2 ~ gamma(1,0.1);
y ~ normal(beta1 + beta2 * x, 1/sqrt(sigma2));
}

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

data <- list(y = y, x = x, N = 100)
fit <- stan(model_code = model,data = data, iter = 4000, warmup = 1000, chains = 1)

amostra[[i]] = extract(fit)

b0[[i]]  = amostra[[i]][["beta1"]]
b1[[i]]  = amostra[[i]][["beta2"]]
tau[[i]] = amostra[[i]][["sigma2"]]

final[[i]] = cbind(amostra[[i]][["beta1"]],amostra[[i]][["beta2"]],amostra[[i]][["sigma2"]])

media_b0[i]  = mean(b0[[i]]) 
media_b1[i]  = mean(b1[[i]]) 
media_tau[i] = mean(tau[[i]]) 

mediana_b0[i]  = median(b0[[i]]) 
mediana_b1[i]  = median(b1[[i]]) 
mediana_tau[i] = median(tau[[i]])  

var_b0[i]  = var(b0[[i]])
var_b1[i]  = var(b1[[i]])
var_tau[i] = var(tau[[i]])

eqm_b0[i]  = mean((b0[[i]]+2)^2)
eqm_b1[i]  = mean((b1[[i]]-2)^2)
eqm_tau[i] = mean((tau[[i]]-.2)^2)

hdi[[i]] = hdi(final[[i]])

}
toc()

media=c(mean(media_b0),mean(media_b1),mean(media_tau))
mediana=c(mean(mediana_b0),mean(mediana_b1),mean(mediana_tau))
var=c(mean(var_b0),mean(var_b1),mean(var_tau))
eqm=c(mean(eqm_b0),mean(eqm_b1),mean(eqm_tau))
medidas=cbind(media,mediana,var,eqm)
rownames(medidas) = c("b0","b1","tau");medidas
int = as.data.frame(cbind(t(as.data.frame((hdi))),rep(c(-2,2,0.2),100)))
write.csv2(int,"C:\\Users\\Rafael\\Desktop\\int.xls")


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


amostra = list()
final = list()
b0  = list()
b1  = list()
tau = list()
media_b0 = vector()
media_b1 = vector()
media_tau = vector()
mediana_b0 = vector()
mediana_b1 = vector()
mediana_tau = vector()
var_b0 = vector()
var_b1 = vector()
var_tau = vector()
eqm_b0  = vector()
eqm_b1  = vector()
eqm_tau = vector()
hdi = list()




tic()
for (i in 1:100){
  y =  yt[[i]][[1]]
  x =  yt[[i]][[2]]
  
  simpleCode1 <- nimbleCode({
    beta0 ~ dnorm(0,0.1)
    beta1 ~ dnorm(0,0.1)
    tau ~ dgamma(1, 0.1)
    for(i in 1:N) {
      Ypred[i] <- beta0 + beta1 * x[i]
      Y[i] ~ dnorm(Ypred[i], tau)
    }
  })
  
  
  


simpleModel1 <- nimbleModel(simpleCode1,
                            data = list(Y = y, x = x),
                            constants = list(N = 100),
                            inits = list(beta0 = 0, beta1 = 0, tau = 2 ))


amostra[[i]] <- nimbleMCMC(code = simpleModel1, 
                       data = list(Y = y, x1 = x1,x2=x2), 
                       inits = list(beta0 = 0, beta1 = 0, tau = 2),
                       monitors=c("beta0","beta1","tau"),
                       nchains = 1, 
                       niter = 10000,
                       thin = 9,
                       nburnin = 1000,
                       summary = TRUE, 
                       progressBar = FALSE)

final[[i]] = amostra[[i]][["samples"]]

b0[[i]]  = final[[i]][,1]
b1[[i]]  = final[[i]][,2]
tau[[i]] = final[[i]][,3]

media_b0[i]  = mean(b0[[i]]) 
media_b1[i]  = mean(b1[[i]]) 
media_tau[i] = mean(tau[[i]]) 

mediana_b0[i]  = median(b0[[i]]) 
mediana_b1[i]  = median(b1[[i]]) 
mediana_tau[i] = median(tau[[i]])  

var_b0[i]  = var(b0[[i]])
var_b1[i]  = var(b1[[i]])
var_tau[i] = var(tau[[i]])

eqm_b0[i]  = mean((b0[[i]]+2)^2)
eqm_b1[i]  = mean((b1[[i]]-2)^2)
eqm_tau[i] = mean((tau[[i]]-.2)^2)

hdi[[i]] = hdi(final[[i]])

}
toc()

media=c(mean(media_b0),mean(media_b1),mean(media_tau))
mediana=c(mean(mediana_b0),mean(mediana_b1),mean(mediana_tau))
var=c(mean(var_b0),mean(var_b1),mean(var_tau))
eqm=c(mean(eqm_b0),mean(eqm_b1),mean(eqm_tau))
medidas=cbind(media,mediana,var,eqm)
rownames(medidas) = c("b0","b1","tau");medidas
int = as.data.frame(cbind(t(as.data.frame((hdi))),rep(c(-2,2,0.2),100)))
write.csv2(int,"C:\\Users\\Rafael\\Desktop\\int.xls")

###########################################################

modelo = list()
amostra = list()
final = list()
b0  = list()
b1  = list()
tau = list()
media_b0 = vector()
media_b1 = vector()
media_tau = vector()
mediana_b0 = vector()
mediana_b1 = vector()
mediana_tau = vector()
var_b0 = vector()
var_b1 = vector()
var_tau = vector()
eqm_b0  = vector()
eqm_b1  = vector()
eqm_tau = vector()
hdi = list()

tic()
for (i in 1:100){
y =  yt[[i]][[1]]
x =  yt[[i]][[2]]
dado=as.data.frame(cbind(y,x))
modelo[[i]] = inla(y~x,family="gaussian",data=dado)

final[[i]] = modelo[[i]]$marginals.fixed

b0[[i]]  = final[[i]][[1]][,1]
b1[[i]]  = final[[i]][[2]][,1]
tau[[i]] = modelo[[i]]$marginals.hyperpar[[1]][,1]

media_b0[i]  = mean(b0[[i]]) 
media_b1[i]  = mean(b1[[i]]) 
media_tau[i] = mean(tau[[i]]) 

mediana_b0[i]  = median(b0[[i]]) 
mediana_b1[i]  = median(b1[[i]]) 
mediana_tau[i] = median(tau[[i]])  

var_b0[i]  = var(b0[[i]])
var_b1[i]  = var(b1[[i]])
var_tau[i] = var(tau[[i]])

eqm_b0[i]  = mean((b0[[i]]+2)^2)
eqm_b1[i]  = mean((b1[[i]]-2)^2)
eqm_tau[i] = mean((tau[[i]]-.2)^2)

hdi[[i]] = hdi(final[[i]])


}
toc()

media=c(mean(media_b0),mean(media_b1),mean(media_tau))
mediana=c(mean(mediana_b0),mean(mediana_b1),mean(mediana_tau))
var=c(mean(var_b0),mean(var_b1),mean(var_tau))
eqm=c(mean(eqm_b0),mean(eqm_b1),mean(eqm_tau))
medidas=cbind(media,mediana,var,eqm)
rownames(medidas) = c("b0","b1","tau");medidas
int = as.data.frame(t(as.data.frame((hdi))))
write.csv2(int,"C:\\Users\\Rafael\\Desktop\\int.xls")


int2 = as.data.frame(t(as.data.frame((hdi(tau)))))
write.csv2(int2,"C:\\Users\\Rafael\\Desktop\\int2.xls")