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?
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.
qbeta for the quantile function in R.)
$\endgroup$
Commented
Apr 19, 2017 at 13:35
@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.
