# Efficiently sampling a thresholded Beta distribution

How should I efficiently sample from the following distribution?

$$x \sim B(\alpha, \beta),\space x > k$$

If $k$ is not too big then rejection sampling may be the best approach, but I am not sure how to proceed when $k$ is large. Perhaps there is some asymptotic approximation which can be applied?

• It's not unambiguously clear what you intend there by "$x \sim B(\alpha, \beta),\space x > k$". Do you mean a truncated beta distribution (truncated on the left at $k$)? – Glen_b -Reinstate Monica Apr 19 '17 at 9:33
• @Glen_b exactly. – user1502040 Apr 19 '17 at 9:38
• For both shape parameters greater than 1 the beta distribution is log-concave, so exponential envelopes can be used for rejection sampling. As to generate untruncated beta variates you're already sampling from truncated exponential distributions (which is easy to do) it should be straightforward to adapt this method. – Scortchi - Reinstate Monica Apr 19 '17 at 10:40

The simplest way, and most general way, that applies to any truncated distribution (it can be also generalized to truncation on both sides), is to use inverse transform sampling. If $F$ is the cumulative distribution of interest, then set $p_0 = F(k)$ and take

$$U \sim \mathcal{U}(p_0, 1) \\ X = F^{-1}(U)$$

where $X$ is a sample from $F$ left-truncated at $k$. The quantile function $F^{-1}$ will map probabilities to samples from $F$. Since we take values of $U$ only from the "area" that matches the values of beta distribution from the non-truncated region, you will sample only those values.

This method is illustrated on the image below where the truncated area is marked by a gray rectangle, points in red are drawn from $\mathcal{U}(p_0, 1)$ distribution and then transformed to $\mathcal{B}(2, 8)$ samples.

• (+1) It's worth remarking that the quantile function isn't so easily evaluated. – Scortchi - Reinstate Monica Apr 19 '17 at 10:56
• @Scortchi If either a or b are 1 or at least an integer, there is a not so bad form (see wikipedia). And in Python there is scipy.special.betainc for the inverse and in R there is pbeta. – Graipher Apr 19 '17 at 13:06
• @Graipher: I've should've said "cheaply, in general" - it'd be better to avoid Newton-Raphson or other iterative solutions if possible. (BTW it's qbeta for the quantile function in R.) – Scortchi - Reinstate Monica Apr 19 '17 at 13:35
• @Scortchi you are right, but in most cases, for modern-day computers this should not be a major problem. I also recommend this approach since it is directly available in most software and can be generalized to any truncated distribution, only if one has access to quantile function. – Tim Apr 19 '17 at 14:04
• Undoubtedly it's good to have an generally applicable, easily implemented method to hand whose run-time doesn't grow with $k$; & for distributions with closed-form quantile functions, e.g. the Weibull, it must be as good as it gets. Nevertheless I suspect $k$ will have to be set to chop off a rather large part of the beta distribution before it beats the efficient rejection-sampling algorithms that are also available in most software & that rely only on calculation of the beta's probability density. – Scortchi - Reinstate Monica Apr 19 '17 at 20:48

@Tim's answer shows how inverse transform sampling can be adapted for truncated distributions, freeing run-time of dependency on the threshold $k$. Further efficiencies can be got by avoiding costly numerical evaluation of the beta quantile function & using inverse transform sampling as part of rejection sampling.

The density function of a beta distribution with shape parameters $\alpha$ & $\beta$ doubly truncated at $k_1<k_2$ (for a little more generality) is

$$f(x) = \frac{x^{(\alpha-1)}(1-x)^{(\beta-1)}}{\operatorname{B}(k2, \alpha, \beta) - \operatorname{B}(k_1, \alpha, \beta)}$$

Take any monotonically increasing part of the density between $x_\mathrm{L}$ & $x_\mathrm{U}$: for $\alpha,\beta>1$ it's log-concave, so you can envelop it with an exponential function drawn at a tangent to any point along it:

$$g(x) = c \cdot \lambda \mathrm{e}^{-\lambda (x-x_\mathrm{L})}$$

Find $\lambda$ by setting the gradients of the log densities equal

$$-\lambda = \frac{a-1}{x} - \frac{b-1}{1-x}$$ & find $c$ by working out how much the exponential density needs to be scaled up to meet the density at that point $$c = \frac{f(x)}{\lambda\mathrm{e}^{-\lambda(x-x_\mathrm{L})}}$$

The best envelope for rejection sampling is the one with the smallest area below it. The area is $$A = c \cdot (1 - \mathrm{e}^{-\lambda(x_\mathrm{U}-x_\mathrm{L})})$$ Substituting expressions in $x$ for $\lambda$ & $c$, & dropping constant factors gives

\begin{align} Q(x)= & \frac{x^a (1-x)^b}{(a+b-2)x - a+1} \cdot\\ & \left[\exp\left(\frac{(b-1)(x-x_L)}{1-x} + \frac{x_L (a-1)}{x} - (a-1)\right) - \right.\\ & \left. \exp\left(\frac{(b-1)(x-x_U)}{1-x} + \frac{x_U(a-1)}{x} - (a-1)\right)\right]\\ \end{align}

Writing the derivative $\frac{\mathrm{d} Q}{\mathrm{d} x}$ is left as an exercise for readers or their computers. Any root-finding algorithm can then be used to find the $x$ for which $\frac{\mathrm{d} Q}{\mathrm{d} x} = 0$.

The same goes, mutatis mutandis, for monotonically decreasing parts of the density. So the truncated beta distribution can be quite neatly enveloped by two exponential functions, say, one from $k_1$ to the mode & one from the mode to $k_2$, ready for rejection sampling. (For a truncated uniform random variable $U$, $\frac{- \log(1-U)}{\lambda}$ has a truncated exponential distribution with rate parameter $\lambda$.)

The beauty of this approach is that all the hard work's in the set up. Once the envelope function's defined, the normalizing constant for the truncated beta density calculated, all that's left is to generate uniform random variates, & perform on them a few simple arithmetical operations, logs & powers, & comparisons. Tightening the envelope function—with horizontal lines or more exponential curves—can of course reduce the number of rejections.

• +1 Nice idea. Because Beta is approximately normal for modest to large values of its parameters, depending on how close to each other they are, using a Gaussian envelope might be quite a bit more efficient yet. – whuber Apr 26 '17 at 20:32
• @whuber: I wanted a distribution with a closed-form quantile distribution for the envelope, but I suppose you could generate the truncated Gaussian variates efficiently using one of those good approximations to the Gaussian quantile function. I'm still interested in what you'd do when $\alpha<1$ or $\beta<1$. – Scortchi - Reinstate Monica Apr 27 '17 at 9:00
• For such small $\alpha$ and $\beta$, you would want to switch to an exponential tail. I'm unsure what you mean by "closed-form," though: when you look hard at computer implementations of exponentials, Gaussians, hypergeometric functions, etc--that is, all the non-algebraic functions--you discover that none of them have a general "closed form": they are computed through successive approximations such as Taylor series, partial fractions, or asymptotic expansions. There's not much difference between computing an exponential and computing a Gaussian quantile. – whuber Apr 27 '17 at 13:10
• @whuber: (1) The approach I've taken here to constructing envelopes wouldn't work because the densities aren't log-concave. (2) (a) I meant certainly algebraic functions + logs & powers, trig. functions if I'd have been asked, & perhaps even gamma functions - I confess I didn't have a precise notion. (b) Point taken - speedy function evaluations aren't confined to those with closed forms. – Scortchi - Reinstate Monica Apr 27 '17 at 16:33
• Good point about failure of log-concavity. I suspect a power law distribution ought to make a good envelope for $\alpha\lt 1$ or $\beta \lt 1$. – whuber Apr 27 '17 at 17:42