Jags
########################################################################################
############# JAGS #############
########################################################################################
set.seed(07042019)
y = rnorm(1000,3,2)
y.censored <- ifelse(y>=0, y, NA)
y.ind <- as.numeric(y>=0)
N = length(y)
model_code= 'model {
for (j in 1:N){
y.ind[j] ~ dinterval(y.censored[j], 0)
y.censored[j] ~ dnorm(a,b)
}
a ~ dnorm(0, .0001)
b ~ dgamma(0.1, .01)
sigma2 <- 1/b
}'
tobit.inits <- function() {
list(a=0, b=1,
y.censored=ifelse(y.ind, NA, -abs(rnorm(N, 0, 1))))
}
m <- jags(model.file=textConnection(model_code),
parameters.to.save = c("a", "b","sigma2"),
data=list(y.censored=y.censored, y.ind=y.ind, N=N),
inits=tobit.inits)
## module glm loaded
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 1935
## Unobserved stochastic nodes: 67
## Total graph size: 2009
##
## Initializing model
hdi(m$BUGSoutput$sims.matrix)
## a b deviance sigma2
## lower 2.940574 0.2274579 3732.848 3.646531
## upper 3.196898 0.2742332 3781.122 4.396418
## attr(,"credMass")
## [1] 0.95
## mean sd 2.5% 25% 50%
## a 3.0640807 0.06413510 2.9350411 3.0227062 3.0634980
## b 0.2504887 0.01192089 0.2267155 0.2423206 0.2505965
## deviance 3754.2112538 12.71349514 3733.0808869 3744.8567090 3753.2475789
## sigma2 4.0013044 0.19186555 3.6544933 3.8669060 3.9904785
## 75% 97.5% Rhat n.eff
## a 3.1054023 3.1931684 1.000596 3000
## b 0.2586047 0.2736357 1.000752 3000
## deviance 3762.6436467 3781.8276063 1.000625 3000
## sigma2 4.1267635 4.4108144 1.000752 3000
Stan
########################################################################################
############# STAN #############
########################################################################################
set.seed(07042019)
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> sigma;
}
model {
y_obs ~ normal(mu, sigma);
y_cens ~ normal(0, 10);
}"
data <- list(y_obs=y, y_cens=y.censored, N_obs = N, N_cens = N, U = 0)
fit <- stan(model_code = model,data = data, iter = 1000, warmup = 100, chains = 1)
##
## SAMPLING FOR MODEL '0a72551b5979edf68dcd95d9ae1d8601' 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: WARNING: There aren't enough warmup iterations to fit the
## Chain 1: three stages of adaptation as currently configured.
## Chain 1: Reducing each adaptation stage to 15%/75%/10% of
## Chain 1: the given number of warmup iterations:
## Chain 1: init_buffer = 15
## Chain 1: adapt_window = 75
## Chain 1: term_buffer = 10
## Chain 1:
## Chain 1: Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 1: Iteration: 100 / 1000 [ 10%] (Warmup)
## Chain 1: Iteration: 101 / 1000 [ 10%] (Sampling)
## Chain 1: Iteration: 200 / 1000 [ 20%] (Sampling)
## Chain 1: Iteration: 300 / 1000 [ 30%] (Sampling)
## Chain 1: Iteration: 400 / 1000 [ 40%] (Sampling)
## Chain 1: Iteration: 500 / 1000 [ 50%] (Sampling)
## Chain 1: Iteration: 600 / 1000 [ 60%] (Sampling)
## Chain 1: Iteration: 700 / 1000 [ 70%] (Sampling)
## Chain 1: Iteration: 800 / 1000 [ 80%] (Sampling)
## Chain 1: Iteration: 900 / 1000 [ 90%] (Sampling)
## Chain 1: Iteration: 1000 / 1000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 1.833 seconds (Warm-up)
## Chain 1: 10.899 seconds (Sampling)
## Chain 1: 12.732 seconds (Total)
## Chain 1:
summary(fit,probs = c(0.025, 0.975),pars=c("mu","sigma"))
## $summary
## mean se_mean sd 2.5% 97.5% n_eff Rhat
## mu 3.046645 0.002128089 0.05951281 2.936888 3.154791 782.0622 0.9991069
## sigma 2.032143 0.001482839 0.04423419 1.949623 2.121730 889.8728 0.9989364
##
## $c_summary
## , , chains = chain:1
##
## stats
## parameter mean sd 2.5% 97.5%
## mu 3.046645 0.05951281 2.936888 3.154791
## sigma 2.032143 0.04423419 1.949623 2.121730
Nimble
########################################################################################
############# NIMBLE #############
########################################################################################
set.seed(07042019)
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, .0001)
b ~ dgamma(0.1, .01)
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)))))
## defining model...
## Warning in nimbleModel(simpleCode1, data = list(y.censored = y.censored, : BUGSmodel: found the same variable(s) in both 'data' and 'constants'; using variable(s) from 'data'.
## building model...
## setting data and initial values...
## Warning in model$setData(data): data not used in model: N
## 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.
mcmc.out <- nimbleMCMC(code = simpleModel1,
monitors=c("a", "b","sigma2"),
nchains = 1,
niter = 1000,
summary = TRUE,
progressBar = TRUE)
## compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compiler details.
## compilation finished.
## running chain 1...
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## Mean Median St.Dev. 95%CI_low 95%CI_upp
## a 3.064161 3.0637196 0.06354095 2.9384860 3.1910948
## b 0.251025 0.2505204 0.01215968 0.2279741 0.2758376
## sigma2 3.992986 3.9916910 0.19301280 3.6253213 4.3864634
Laplace
########################################################################################
############# LAPLACE #############
########################################################################################
d = ifelse(y < 0, 0, 1)
model <- function(p,y,d) {
log_lik <- pnorm(0, p["mu"], p["sigma"], log = T)*sum(1-d) + sum(d * dnorm(y, p["mu"], p["sigma"], log = T))
log_post <- log_lik + dnorm(p["mu"], 0, 10, log = T) + dunif(p["sigma"], 0, 100, log = T)
log_post
}
inits <- c(mu = 0, sigma = 10)
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)
## mu sigma
## 3.06 2.00
## mu sigma
## mu 4.037578e-03 -9.679021e-05
## sigma -9.679021e-05 2.216696e-03