Inverse CDFs for truncated RNG distributions

Modeling
1

I’m trying to fit a model that involves several truncated distributions - at the moment, a bounded normal distribution and an upper-bounded poisson. I want to write my own truncated X_rng for different distributions - poisson, negative binomial, student’s-t, etc - similar to normal_lb_ub_rng from this post. Are there other inverse CDFs available in Stan, similar to inv_phi? If not - is there any other approach for generating truncated quantities for posterior predictive checks?

2

Boost has some, but you would have to write your own C++.

3

You can also do rejection sampling - draw from the full distribution, if the result is outside the bounds, draw a new one until it is. This is easy to implement and tends to be very efficient if your bounds don’t exclude a large fraction of the mass (even when rejecting say 2/3 of samples on average, it may easily still be faster to draw 3 times more samples than to compute the inverse CDF).

2 Likes
4

Yes, this is pretty much what I did. Something along these lines:

generated quantities{
    y_pred = upper_bound+1;
    while y_pred > upper_bound
         y_pred ~ poisson_rng(lam);
}

It works, slows things down but not by much… just wanted to make sure I’m not missing a more “natural”, Stan-based solution. Thanks!

1 Like
5

Currently the most efficient method to construct truncated RNGS for continuous densities is to use the algebraic solver to implicitly solve for the quantile function. For discrete densities you can sum up the atomic probabilities explicitly to compute the quantile function.

2 Likes
6

Thanks! Do you know of any tutorials/examples of this I can learn from? I’m new to Stan, and this seems like something I would love to learn!

7

Here’s code that samples from a Poisson truncated at U,

int y;

real sum_p = 0;
real u = uniform_rng(0, 1);

for (b in 0:U) {
  sum_p = sum_p + exp(poisson_lpmf(b | lambda) - poisson_lcdf(U | lambda));
  if (sum_p >= u) {
    y = b;
    break;
  }
}

I don’t have any Stan code for sampling from a continuous truncated immediately available, but the pseudo code would be

  • Compute p_low and p_high corresponding to the quantiles of the lower and upper truncated points. For a Gaussian that would be p_low = normal_lcdf(x_low | mu, sigma), etc.
  • Sample a uniform value between p_low and p_high with real u = uniform_rng(p_low, p_high).
  • Find the x whose quantile is equal to u. If the inverse CDF is available in closed form then this is easy, but if you only have the CDF available then you have to solve for u = target_cdf(x) using the algebraic solver.
2 Likes
8

Thanks!