R Notebook

library("R2OpenBUGS")
library("R2jags")
library("R2stan")
library("dplyr")
library("rstan")
library("nimble")
library("HDInterval")
library("tictoc")
library("mvtnorm")


########################################################################################
#############                          DADOS                               ############# 
########################################################################################

set.seed(6)

y  = list()
y.censt = list()
y.obs = list()
y.ind = list()
censura = vector()

for (i in 1:100){
  n=100
  y[[i]] <- rnorm(100,3,2)
  y.censt[[i]] <- ifelse(y[[i]]>0,y[[i]],NA)}

for (i in 1:100){y.ind[[i]] <- as.numeric(y[[i]]>=0)}

for (i in 1:100){censura[i] = as.vector(table(y.censt[[i]]>0)/100)};mean(censura)


N = length(y)


########################################################################################
#############                          BUGS                                ############# 
########################################################################################

amostra_bugs = list()
final_bugs = list()
mu_bugs  = list()
tau_bugs = list()
media_mu_bugs = vector()
media_tau_bugs = vector()
sd_mu_bugs = vector()
sd_tau_bugs = vector()
eqm_mu_bugs  = vector()
eqm_tau_bugs  = vector()
hdi_bugs = list()


tic()
for (i in 1:100){
  
y.cens = y.censt[[i]]

sink("regression.txt")
cat(" 
  model
  {
  # priors
  a ~ dnorm(0,0.01)
  b ~ dgamma(1,10)
  
  # likelihood
  for(i in 1:N)
  {
  y[i] ~ dnorm(a, b)C(0,1)}
  }
  ")
sink()

my.data <- list(y=y.cens, N=100)

params <- c("a", "b")

inits <- function(){list(beta0=1,beta1=1)}

amostra_bugs[[i]] <- bugs(data=my.data, inits=inits, parameters.to.save=params,
            n.iter=10000, n.chains=1, n.burnin=1000,
              model.file="regression.txt",DIC=FALSE)  
  

final_bugs[[i]] = amostra_bugs[[i]]$sims.matrix

mu_bugs[[i]]  = final_bugs[[i]][,1]
tau_bugs[[i]]  = final_bugs[[i]][,2]

media_mu_bugs[i]  = mean(mu_bugs[[i]]) 
media_tau_bugs[i]  = mean(tau_bugs[[i]]) 

sd_mu_bugs[i]  = sd(mu_bugs[[i]])
sd_tau_bugs[i]  = sd(tau_bugs[[i]])

eqm_mu_bugs[i]  = mean((mu_bugs[[i]]-3)^2)
eqm_tau_bugs[i]  = mean((tau_bugs[[i]]-0.25)^2)

hdi_bugs[[i]] = hdi(final_bugs[[i]])

}
toc()

t_bugs = 162 

media_bugs=c(mean(media_mu_bugs),mean(media_tau_bugs))
sd_bugs=c(mean(sd_mu_bugs),mean(sd_tau_bugs))
eqm_bugs=c(mean(eqm_mu_bugs),mean(eqm_tau_bugs))
medidas_bugs=cbind(media_bugs,sd_bugs,eqm_bugs)
rownames(medidas_bugs) = c("mu","tau")

hdp_bugs=t(as.data.frame(hdi(final_bugs)))
aux = as.character(rownames(hdp_bugs))
real=rep(c(3,0.25),10)
hdp_bugs=as.data.frame(cbind(aux,hdp_bugs,real))

hdp_bugs$lower = as.numeric(as.character(hdp_bugs$lower))
hdp_bugs$upper = as.numeric(as.character(hdp_bugs$upper))
hdp_bugs$real = as.numeric(as.character(hdp_bugs$real))

hdp_bugs = arrange(hdp_bugs,aux)

prop_bugs=ifelse(hdp_bugs$lower<hdp_bugs$real&hdp_bugs$real<hdp_bugs$upper,1,0)

hdp_bugs=as.data.frame(cbind(hdp_bugs,prop_bugs))

mu_prop_bugs=sum(hdp_bugs$prop_bugs[1:100])/100;tau_prop_bugs=sum(hdp_bugs$prop_bugs[101:200])/100
cred_bugs = rbind(mu_prop_bugs,tau_prop_bugs);colnames(cred_bugs)=c("cred_bugs")

tabela_bugs = cbind(medidas_bugs,cred_bugs,t_bugs)

########################################################################################
#############                          jags                                ############# 
########################################################################################

amostra_jags = list()
final_jags = list()
mu_jags  = list()
tau_jags = list()
media_mu_jags = vector()
media_tau_jags = vector()
sd_mu_jags = vector()
sd_tau_jags = vector()
eqm_mu_jags  = vector()
eqm_tau_jags  = vector()
hdi_jags = list()

tic()
for (i in 1:100){
  
  y.cens = y.censt[[i]]
  y.indx  = y.ind[[i]]

  model_code=  'model {
  for (j in 1:N){
  y.indx[j] ~ dinterval(y.cens[j], 0)
  y.cens[j] ~ dnorm(a,b)
  }
  a ~ dnorm(0,0.01)
  b ~ dgamma(1,10)
}'

tobit.inits <- function(){list(a=0, b=1,y.cens=ifelse(y.indx, NA, -abs(rnorm(N, 0, 1))))}
  
amostra_jags[[i]] <- jags(model.file=textConnection(model_code),
            parameters.to.save = c("a", "b"),
            data=list(y.cens=y.cens, y.indx=y.indx, N=N),
            inits=tobit.inits,DIC = FALSE,n.iter = 10000,n.burnin = 1000)

final_jags[[i]] = amostra_jags[[i]][["BUGSoutput"]][["sims.matrix"]]

mu_jags[[i]]  = final_jags[[i]][,1]
tau_jags[[i]]  = final_jags[[i]][,2]

media_mu_jags[i]  = mean(mu_jags[[i]]) 
media_tau_jags[i]  = mean(tau_jags[[i]]) 

sd_mu_jags[i]  = sd(mu_jags[[i]])
sd_tau_jags[i]  = sd(tau_jags[[i]])

eqm_mu_jags[i]  = mean((mu_jags[[i]]-3)^2)
eqm_tau_jags[i]  = mean((tau_jags[[i]]-0.25)^2)

hdi_jags[[i]] = hdi(final_jags[[i]])

}
toc()

t_jags = 13

media_jags=c(mean(media_mu_jags),mean(media_tau_jags))
sd_jags=c(mean(sd_mu_jags),mean(sd_tau_jags))
eqm_jags=c(mean(eqm_mu_jags),mean(eqm_tau_jags))
medidas_jags=cbind(media_jags,sd_jags,eqm_jags)
rownames(medidas_jags) = c("mu","tau")

hdp_jags=t(as.data.frame(hdi(final_jags)))
aux = as.character(rownames(hdp_jags))
real=rep(c(3,0.25),10)
hdp_jags=as.data.frame(cbind(aux,hdp_jags,real))

hdp_jags$lower = as.numeric(as.character(hdp_jags$lower))
hdp_jags$upper = as.numeric(as.character(hdp_jags$upper))
hdp_jags$real = as.numeric(as.character(hdp_jags$real))

hdp_jags = arrange(hdp_jags,aux)

prop_jags=ifelse(hdp_jags$lower<hdp_jags$real&hdp_jags$real<hdp_jags$upper,1,0)

hdp_jags=as.data.frame(cbind(hdp_jags,prop_jags))

mu_prop_jags=sum(hdp_jags$prop_jags[1:100])/100;tau_prop_jags=sum(hdp_jags$prop_jags[101:200])/100
cred_jags = rbind(mu_prop_jags,tau_prop_jags);colnames(cred_jags)=c("cred_jags")

tabela_jags = cbind(medidas_jags,cred_jags,t_jags)


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

fit = list()
amostra_stan = list()
final_stan = list()
mu_stan  = list()
tau_stan = list()
media_mu_stan = vector()
media_tau_stan = vector()
sd_mu_stan = vector()
sd_tau_stan = vector()
eqm_mu_stan  = vector()
eqm_tau_stan  = vector()
hdi_stan = list()

tic()
for (i in 1:100){
  y = rnorm(1000,3,2)
  y.censored <- ifelse(y>=0, y, NA)
  y.ind <- as.numeric(y>=0)
  N = length(y)
  
  model = "data {
int<lower=0> N_obs;
int<lower=0> N_cens;
real y_obs[N_obs];
real<lower=min(y_obs)> U;
}
parameters {
real<lower=U> y_cens[N_cens];
real mu;
real<lower=0> sigma2;
}
transformed parameters {
real<lower=0> tau;
tau = 1 / sigma2; 
}
model {
y_obs ~ normal(mu, sqrt(sigma2));
y_cens ~ normal(0, 10);
mu ~ normal(0,10);
sigma2 ~ gamma(1,10);
}"


data <- list(y_obs=y, y_cens=y.censored, N_obs = N, N_cens = N, U = 0)

fit[[i]] <- stan(model_code = model,data = data, iter = 10000, warmup = 1000, chains = 1)

amostra_stan[[i]] <- extract(fit[[i]])

mu_stan[[i]]  = amostra_stan[[i]][["mu"]]
tau_stan[[i]]  = amostra_stan[[i]][["tau"]]

final_stan[[i]] = cbind(amostra_stan[[i]][["mu"]],amostra_stan[[i]][["tau"]])

mu_stan[[i]]  = final_stan[[i]][,1]
tau_stan[[i]]  = final_stan[[i]][,2]

media_mu_stan[i]  = mean(mu_stan[[i]]) 
media_tau_stan[i]  = mean(tau_stan[[i]]) 

sd_mu_stan[i]  = sd(mu_stan[[i]])
sd_tau_stan[i]  = sd(tau_stan[[i]])

eqm_mu_stan[i]  = mean((mu_stan[[i]]-3)^2)
eqm_tau_stan[i]  = mean((tau_stan[[i]]-0.25)^2)

hdi_stan[[i]] = hdi(final_stan[[i]])}
toc()

t_stan = 3803 

media_stan=c(mean(media_mu_stan),mean(media_tau_stan))
sd_stan=c(mean(sd_mu_stan),mean(sd_tau_stan))
eqm_stan=c(mean(eqm_mu_stan),mean(eqm_tau_stan))
medidas_stan=cbind(media_stan,sd_stan,eqm_stan)
rownames(medidas_stan) = c("mu","tau")

hdp_stan=t(as.data.frame(hdi(final_stan)))
aux = as.character(rownames(hdp_stan))
real=rep(c(3,0.25),10)
hdp_stan=as.data.frame(cbind(aux,hdp_stan,real))


hdp_stan$lower = as.numeric(as.character(hdp_stan$lower))
hdp_stan$upper = as.numeric(as.character(hdp_stan$upper))
hdp_stan$real = as.numeric(as.character(hdp_stan$real))

hdp_stan = arrange(hdp_stan,aux)

prop_stan=ifelse(hdp_stan$lower<hdp_stan$real&hdp_stan$real<hdp_stan$upper,1,0)

hdp_stan=as.data.frame(cbind(hdp_stan,prop_stan))

mu_prop_stan=sum(hdp_stan$prop_stan[1:43])/44;tau_prop_stan=sum(hdp_stan$prop_stan[44:88])/44
cred_stan = rbind(mu_prop_stan,tau_prop_stan);colnames(cred_stan)=c("cred_stan")

tabela_stan = cbind(medidas_stan,cred_stan,t_stan)

rm(fit)

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

amostra_nimble = list()
final_nimble = list()
mu_nimble  = list()
tau_nimble = list()
media_mu_nimble = vector()
media_tau_nimble = vector()
sd_mu_nimble = vector()
sd_tau_nimble = vector()
eqm_mu_nimble  = vector()
eqm_tau_nimble  = vector()
hdi_nimble = list()

tic()
for (i in 1:100){
y = rnorm(1000,3,2)
y.censored <- ifelse(y>=0, y, NA)
y.ind <- as.numeric(y>=0)
N = length(y)

simpleCode1 <- nimbleCode({
  for (j in 1:N){
    y.ind[j] ~ dinterval(y.censored[j], 0)
    y.censored[j] ~ dnorm(a,b)
  }
  a ~ dnorm(0, .01)
  b ~ dgamma(1,10)
  sigma2 <- 1/b})

tobit.inits <- function() {
  list(a=0, b=1, 
       y.censored=ifelse(y.ind, NA, -abs(rnorm(N, 0, 1))))
}

simpleModel1 <- nimbleModel(simpleCode1,
                            data=list(y.censored=y.censored, y.ind=y.ind, N=N),
                            constants = list(N = N),
                            inits = list(a=0, b=1,y.censored=ifelse(y.ind, NA, -abs(rnorm(N, 0, 1)))))

amostra_nimble[[i]] <- nimbleMCMC(code = simpleModel1, 
                       monitors=c("a", "b"),
                       nchains = 1, 
                       niter = 1000,
                       summary = TRUE, 
                       progressBar = TRUE)

final_nimble[[i]] = amostra_nimble[[i]][["samples"]]

mu_nimble[[i]]  = final_nimble[[i]][,1]
tau_nimble[[i]]  = final_nimble[[i]][,2]

media_mu_nimble[i]  = mean(mu_nimble[[i]]) 
media_tau_nimble[i]  = mean(tau_nimble[[i]]) 

sd_mu_nimble[i]  = sd(mu_nimble[[i]])
sd_tau_nimble[i]  = sd(tau_nimble[[i]])

eqm_mu_nimble[i]  = mean((mu_nimble[[i]]-3)^2)
eqm_tau_nimble[i]  = mean((tau_nimble[[i]]-0.25)^2)

hdi_nimble[[i]] = hdi(final_nimble[[i]])

}
toc()

t_nimble = 5787

media_nimble=c(mean(media_mu_nimble),mean(media_tau_nimble))
sd_nimble=c(mean(sd_mu_nimble),mean(sd_tau_nimble))
eqm_nimble=c(mean(eqm_mu_nimble),mean(eqm_tau_nimble))
medidas_nimble=cbind(media_nimble,sd_nimble,eqm_nimble)
rownames(medidas_nimble) = c("mu","tau")

hdp_nimble=t(as.data.frame(hdi(final_nimble)))
aux = as.character(rownames(hdp_nimble))
real=rep(c(3,0.25),10)
hdp_nimble=as.data.frame(cbind(aux,hdp_nimble,real))

hdp_nimble$lower = as.numeric(as.character(hdp_nimble$lower))
hdp_nimble$upper = as.numeric(as.character(hdp_nimble$upper))
hdp_nimble$real = as.numeric(as.character(hdp_nimble$real))

hdp_nimble = arrange(hdp_nimble,aux)

prop_nimble=ifelse(hdp_nimble$lower<hdp_nimble$real&hdp_nimble$real<hdp_nimble$upper,1,0)

hdp_nimble=as.data.frame(cbind(hdp_nimble,prop_nimble))

mu_prop_nimble=sum(hdp_nimble$prop_nimble[1:100])/100;tau_prop_nimble=sum(hdp_nimble$prop_nimble[101:200])/100
cred_nimble = rbind(mu_prop_nimble,tau_prop_nimble);colnames(cred_nimble)=c("cred_nimble")

tabela_nimble = cbind(medidas_nimble,cred_nimble,t_nimble)



########################################################################################
#############                          LAPLACE                             ############# 
########################################################################################

amostra_laplace = list()
final_laplace = list()
mu_laplace  = list()
tau_laplace = list()
media_mu_laplace = vector()
media_tau_laplace = vector()
sd_mu_laplace = vector()
sd_tau_laplace = vector()
eqm_mu_laplace  = vector()
eqm_tau_laplace  = vector()
hdi_laplace = list()


for (i in 1:100){
y <- rnorm(100,3,2)
y.cens <- ifelse(y>0,y,NA)
d = ifelse(y < 0, 0, 1)

model <- function(p,y,d) {
  log_lik <- pnorm(0, p["mu"], 1/sqrt(p["tau"]), log = T)*sum(1-d) + sum(d * dnorm(y, p["mu"], 1/sqrt(p["tau"]), log = T))
  log_post <- log_lik + dnorm(p["mu"], 0, 10, log = T) + dgamma(p["tau"], 1, 10, log = T)
  log_post
}

inits <- c(mu = 0, tau = 1)
fit <- optim(inits, model, control = list(fnscale = -1), hessian = TRUE, y = y, d = d)
param_mean <- fit$par
param_cov_mat <- solve(-fit$hessian)
round(param_mean, 2)

final_laplace[[i]] <- rmvnorm(10000, param_mean, param_cov_mat)

mu_laplace[[i]]  = final_laplace[[i]][,1]
tau_laplace[[i]]  = final_laplace[[i]][,2]

media_mu_laplace[i]  = mean(mu_laplace[[i]]) 
media_tau_laplace[i]  = mean(tau_laplace[[i]]) 

sd_mu_laplace[i]  = sd(mu_laplace[[i]])
sd_tau_laplace[i]  = sd(tau_laplace[[i]])

eqm_mu_laplace[i]  = mean((mu_laplace[[i]]-3)^2)
eqm_tau_laplace[i]  = mean((tau_laplace[[i]]-0.25)^2)

hdi_laplace[[i]] = hdi(final_laplace[[i]])
}


t_laplace = 1

media_laplace=c(mean(media_mu_laplace),mean(media_tau_laplace))
sd_laplace=c(mean(sd_mu_laplace),mean(sd_tau_laplace))
eqm_laplace=c(mean(eqm_mu_laplace),mean(eqm_tau_laplace))
medidas_laplace=cbind(media_laplace,sd_laplace,eqm_laplace)
rownames(medidas_laplace) = c("mu","tau")

hdp_laplace=t(as.data.frame(hdi(final_laplace)))
aux = as.character(rownames(hdp_laplace))
real=rep(c(3,0.25),10)
hdp_laplace=as.data.frame(cbind(aux,hdp_laplace,real))

hdp_laplace$lower = as.numeric(as.character(hdp_laplace$lower))
hdp_laplace$upper = as.numeric(as.character(hdp_laplace$upper))
hdp_laplace$real = as.numeric(as.character(hdp_laplace$real))

hdp_laplace = arrange(hdp_laplace,aux)

prop_laplace=ifelse(hdp_laplace$lower<hdp_laplace$real&hdp_laplace$real<hdp_laplace$upper,1,0)

hdp_laplace=as.data.frame(cbind(hdp_laplace,prop_laplace))

mu_prop_laplace=sum(hdp_laplace$prop_laplace[1:100])/100;tau_prop_laplace=sum(hdp_laplace$prop_laplace[101:200])/100
cred_laplace = rbind(mu_prop_laplace,tau_prop_laplace);colnames(cred_laplace)=c("cred_laplace")

tabela_laplace = cbind(medidas_laplace,cred_laplace,t_laplace)