Sampling from Truncated Distribution in STAN

Question

Based on section 26.3 of the stan user guide, I'm trying to specify a model in which observed values are rounded and the true values are known to fall in a certain range (between observed and observed -1).

The STAN code is below. Looking at the traceplot for xtrue the sampled values are not constrained between (xobs-1,xobs). Any help is appreciated.

library(rstan)
rstan_options(auto_write = FALSE)
options(mc.cores = parallel::detectCores())



nobs=10
xtrue=runif(nobs,0,5)
xobs=ceiling(xtrue+rnorm(nobs,0,1))

dat=list(N=length(xobs),x=xobs)

init_fun <- function() {list(xtrue=xobs-.5) }

m="
data {
 int<lower = 1> N;
 real x[N];
}

parameters {
 real xtrue[N]; 
}

model{

    for(i in 1:N){
            increment_log_prob(normal_log(xtrue[i], x[i], 1));
            increment_log_prob(-log_diff_exp(normal_cdf_log(x[i],0,1),
                normal_cdf_log(x[i]-1,0,1)));

    }

}

"

fit=stan(model_code=m, data = dat,iter = 2000, chains = 1,thin=3,init=init_fun)

parms=extract(fit,c('xtrue'))
xtrue <- colMeans(parms[['xtrue']])

head(xobs)
head(xtrue)
traceplot(fit)

Answer 1

Given that xtrue[i]'s are constrained, Stan requires that these constraints are included in the variable declaration. To my knowledge, these constraints must be scalar quantities.

Below, I worked around this requirement by considering auxillary parameters, xraw[i], which have a truncated normal distribution.

m <- "
data {
  int<lower = 1> N;
  real x[N];
}

parameters {
  real<lower=-1, upper=0> xraw[N]; 
}

transformed parameters {
  real xtrue[N];
  for(i in 1:N)
  xtrue[i] = xraw[i] + x[i];
}

model{
  for(i in 1:N){
    target += normal_lpdf(xraw[i]| 0, 1);
    target += -log_diff_exp(normal_lcdf(0| 0, 1), normal_lcdf(-1| 0, 1));

  }
}
"

library(rstan)
rstan_options(auto_write = FALSE)
options(mc.cores = parallel::detectCores())

nobs=10

xtrue=runif(nobs,0,5)
xobs=ceiling(xtrue+rnorm(nobs,0,1))

dat=list(N=length(xobs),x=xobs)

init_fun <- function() {list(xtrue=xobs-.5) }

mod <- stan_model(model_code = m)
s <- sampling(mod, data = dat, iter = 2000, chains = 1, thin = 3, init = init_fun)
fit=stan(model_code=m, data = dat,iter = 2000, chains = 1,thin=3,init=init_fun)

parms=extract(s,c('xtrue'))
xtrue <- colMeans(parms[['xtrue']])

head(xobs)
[1] 4 2 5 6 4 2   
head(xtrue)
[1] 3.533775 1.507112 4.561159 5.545677 3.538002 1.520043
par(mfrow = c(2,5))
for(i in 1:10) {
  hist(samples$xtrue[,i], prob=T, main = paste(c("xtrue[",i,"]"), collapse=""), xlab=NULL)
  curve(dnorm(x, xobs[i], 1)/(.5 - pnorm(-1)), add=T, lty=2)
}

The posterior draws appear to follow the correct distributions: