Sampling from a truncated distribution

[Glen DePalma]

I posted this on Stack Exchange but thought it might be more useful here.

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)

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)

[Kevin Van Horn]

According to the data generation code you give

    nobs=10

    xtrue=runif(nobs,0,5)

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

xtrue is, in fact, not guaranteed to lie between xobs and xobs - 1.

Leaving that aside, if you want to constrain a variable to lie in a particular range, you have to make that explicit. For  example to constrain xtrue to lie between x and x - 1:

    parameters {

      vector<lower=0, upper=1>[N] pre_xtrue;

    }

    derived parameters {

      vector[N] xtrue;

      xtrue <- x - pre_xtrue;

[Glen DePalma]

Thanks but I'm not sure why what I was doing doesn't work..  And, yes, I am aware the statements do not guarantee the true value is above the observed.

[Bob Carpenter]

You'll also run into both statistical and computational issues
if you have interval constraints, but the values aren't truly
constrained, when all of the mass piles up at a boundary.

- Bob

P.S. I really wish one "stackexchange" would win the battle. I
hadn't signed up for this one because I wasn't expecting computational
questions on the stats thingy.



> On Mar 20, 2016, at 1:11 AM, Glen DePalma <glen.d...@gmail.com> wrote:
>
> Thanks this is what I was looking for. And, yes, I am aware the statements do not guarantee the true value is above the observed.

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Glen DePalma]

Thanks Bob, and for all your work.  Do you have a recommendation how to code latent xtrue so it is between xobs-1 and xobs?  

And further how to calculate the likelihood function?  I was thinking something like:

 increment_log_prob(normal_log(xobs,xtrue,xsig)-normal_log(xobs-1,xtrue,xsig));

[Bob Carpenter]

Yes, but you'll have to do it yourself with the Jacobians.

data {
int<lower=0> N;
ordered[N] xobs; // make sure these go up in order
...

parameters {
real xtrue_raw[N];
...

transformed parameters {
real xtrue[N];
for (n in 1:N) {
xtrue[n - 1] <- x_obs[n] + inv_logit(xtrue_raw[n]);
}

model {
// Jacobian adjustment for xtrue_raw -> xtrue transform
for (n in 1:N)
increment_log_prob(log_inv_logit(xtrue_raw[n])
+ log1m_inv_logit(x_true_raw[n]));


There's a discussion in the manual chapter on transforms that
shows this is what we do internally for upper and lower bounds.

- Bob