# Sampling from Truncated Distribution in STAN

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']])

traceplot(fit)


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']]) 