Generating random values from truncated Poisson?

As part of a prior predictive simulation, I would like to use Stan to generate values from a hurdle model, i.e. a Bernoulli trial to predict whether the outcome is 0 and a truncated Poisson for outcomes >= 1.

Does anyone have advice on how I can invert the Poisson CDF similar to R’s qpois function? The examples in the manual are for continuous functions where there’s a sensible analytical solution.

functions {
  int poisson_trunc_rng(real lambda) {
    real p = poisson_cdf(0, lambda);
    real u = uniform_rng(p, 1);
    // In R, the following line would be
    int y = qpois(u, lambda);
    return y;
  }
}```

To generate 0 truncated poisson rv the wikipedia article says you can do:

  init:
         Let k ← 1, t ← e−λ / (1 - e−λ) * λ, s ← t.
         Generate uniform random number u in [0,1].
    while s < u do:
         k ← k + 1.
         t ← t * λ / k.
         s ← s + t.
    return k.

In terms of Stan this would look like

functions {
  int zerotrunc_poisson_rng (real lambda) {
    int k = 1;
    real t = lambda * exp(-lambda) / (1 - exp(-lambda));
    real s = t;
    real u = uniform_rng(0, 1);
    
    while (s < u) {
      k += 1;
      t *= lambda / k;
      s += t;
    }
   
   return k; 
  }

}
1 Like