RNG for truncated distributions

Modeling
,
23

I’m looking for a way to do the same for the beta distribution: to sample from the beta distribution with an upper bound. But it seems that stan does not have a built-in inverse cdf for the beta distribution. Any help is much appreciated!

real beta_ub_rng(real alpha, real beta, real ub) {
    real p_ub = beta_cdf(ub, alpha, beta);  // cdf for upper bound
    real u = uniform_rng(0, p_ub);
    return ????  // inverse cdf for value
}
24

There is one in Boost but you would have to call it from C++.

25

Thanks! I was also wondering for gamma: How would I sample from gamma with lower or upper bound?
This is what I have but I’m too new to this to be sure that I’m doing it right:

real gamma_lb_rng(real alpha, real beta, real lb) {
    real p_lb = gamma_cdf(lb, alpha, beta);
    real u = uniform_rng(p_lb, 1);
    real y  = gamma_cdf(u, alpha, beta);
    return y;
}
3 Likes
26

Unfortunately, that’s not correct, you need
real y = gamma_inverse_cdf(u, alpha, beta);

Where gamma_inverse_cdf is the inverse CDF (also called quantile functions and PPF) of the gamma function, which AFAIK is not available in Stan. But if you are willing to invoke C++, it should be easy to implement it using boost (e.g. as discussed at Quantile functions in boost (C++) - Stack Overflow).

Also note that there is a possibility for confusion as there is also “inverse gamma” distribution, so that “CDF of (inverse gamma distribution)” is something completely different than “(inverse CDF) of gamma distribution” (which is the one you want).

Best of luck!

3 Likes
27

In the Stan user manual Phi(u) is used rather than inv_Phi(u). Which is correct?
link: https://mc-stan.org/docs/2_18/stan-users-guide/truncated-random-number-generation.html

image
image1048×447 24.5 KB

28

inv_Phi note that the current version (https://mc-stan.org/docs/2_23/stan-users-guide/truncated-random-number-generation.html) has this fixed.

29

Thanks @martinmodrak! so inv_gamma_cdf() which exists in stan is the cdf of inverse gamma distribution, which is NOT what I want, correct? Thanks for pointing this out!

30

Correct. Also, instead of dipping into C++ it might just be easier to prepare the necessary data in Stan but call the inverse CDF from R or whatever language you invoke Stan from (in R it is called qgamma)