RNG for truncated distributions


#1

So I have a truncated normal in my model. But how do I do a posterior predictive check with such a model? I know there are two ways to force a truncated distribution in the model block, either by parameter bounds or by explicit truncation wtih the T[a,b] modifier. But the explicit truncation approach does not work for the _rng functions. So the question is: how do I draw samples from truncated distributions using _rng functions in the generated quantities block?

Or is there a simple trick to transform draws from non-truncated distribution to a truncated one, similar to the desugaring of the T[a,b] syntax for sampling statements shown in the docs? The best I have so far is rejection sampling which luckily is not slowing my model down a lot, but it feels hacky.

Thanks for any ideas.


#2

I would use brute force. You sample from the rng and if it violates your truncation, then you sample again until you are happy. Not nice, but should work.


#3

Yes, that’s what I am doing (maybe I should’ve clarified this is what I mean by “rejection sampling”), I just wondered whether there is a better way.


#4

If the PRNG uses the inverse CDF method of drawing randomly from a standard uniform and then transforming that by the inverse CDF of the desired distribution, then it is often relatively straightforward to draw uniformally from some subset of the [0,1] interval and transform that by the inverse CDF of the distribution in order to obtain a truncated draw.


#5

But I guess Stan does not support that right now - is it correct?


#6

There is no syntax for it in the Stan language, but you can draw from a uniform distribution and there are some inverse CDFs implemented.


#7

For concreteness, I’ll give the example I worked on this morning. I want to draw from a left-truncated Weibull to impute right-censored event times. Let’s say that I want a Weibull(alpha, sigma) truncated so that my random draws are never below some lower boundary t.

Then I should:
(1) Figure out the quantile corresponding to lower cutoff t. Let’s say it’s 0.2, or the 20th percentile. Call that p.
(2) Draw from a continuous uniform(0.2, 1). Call that draw u.
(3) Apply the inverse transform for the cdf, or “percent point” function. I did not feel like doing algebra on a Monday morning so I found mine at https://www.wolframalpha.com/input/?i=inverse+cdf+weibull. This will return a draw from the correct truncated weibull.

(The unstated step 4 is to post on the Stan discourse to check my logic. So far the draws look good when I plot them, but I have not checked it much beyond that.)

For a double truncation, you would add a step to calculate the cdf at the upper truncation point and draw from a uniform(p1, p2) instead of uniform(p, 1).

functions {
  real tweibull_rng(real alpha, real sigma, real t) {
    real p;
    real u;
    real x;
    p = weibull_cdf(t, alpha, sigma);
    u = uniform_rng(p, 1);
    x = sigma * (-log1m(u))^(1/alpha);
    return x;
  }
}


Generated quantites block
#8

Thanks, that’s exactly what I need to do and with a lot of detail.


#9

This is very nice. I’m going to add this example to the manual after a bit of renaming and collapsing:

real weibull_lb_rng(real alpha, real sigma, real lb) {
    real p = weibull_cdf(lb, alpha, sigma);  // cdf for bounds
    real u = uniform_rng(p, 1);
    return sigma * (-log1m(u))^(1 / alpha);  // inverse cdf for value
  }

Here’s the issue comment


#10

I’m working on the same problem (truncated normal in RNG block). The Weibull example makes sense, but would you be able to share your code for the truncated normal? I’d like to check my work.


#11

It’d be something like this for a normal RNG with truncation at a lower bound lb:

real normal_lb_rng(real mu, real sigma, real lb) {
    real p = normal_cdf(lb, mu, sigma);  // cdf for bounds
    real u = uniform_rng(p, 1);
    return (sigma * inv_Phi(u)) + mu;  // inverse cdf for value
}

The inv_Phi function gives the inverse of a Normal(0,1) cdf, and you can apply the scale parameter and location shift to get back to a Normal(mu, sigma).