generate random number from exponential distribution with random right truncation

Question

I need to draw a random number from an exponential distribution (rate $mu$) that is right-truncated with the truncation value coming from a gamma distribution (shape $k$, rate $lambda$).

My naive approach would have been to first draw the truncation value, and then draw from the truncated exponential distribution via Inverse Transform Sampling. However, such draws do not follow the same distribution which I get from simulation results.

The following R code demonstrates that this approach is for some reason not working:

    # parameter settings
    k <- 2
    lambda <- 1
    mu <- 0.1
    n <- 1e7 # number of draws

    # draw right truncation values
    trunc_value <- rgamma(n=n, shape=k, rate=lambda)

    # draw from right-truncated exponential via inverse transform sampling
    itexp <- function(u, rate, trunc) { -log(1-u*(1-exp(-trunc*rate)))/rate }
    rtexp <- function(n, rate, trunc) { itexp(runif(n), rate, trunc) }
    x0 <- rtexp(n=n, rate=mu, trunc_value)

    # simulate by drawing from exponential and then discarding values larger than trunc
    x1 <- rexp(n, mu)
    x1 <- x1[x1 < trunc_value]

    # means differ significantly
    mean(x0) # 0.95
    mean(x1) # 1.34

Answer 1

trunc_value has a mean of about $2$ by construction.

But if you take those which were greater than the original x1 you will find those used for the the non-rejected x1 had a mean of about $2.86$.

So your rejection method has biased trunc_value upwards and so biased the non-rejected x1 values upwards.

Answer 2

I now found a way to draw samples so that they match x1 as stated in my question.

One can specify the likelihood of the truncated random variable by multiplying the non-truncated probability distribution with the probability of the truncation limit being larger: $$L(x | k,\lambda,\mu) = e^{-\mu x} \cdot P(\text{trunc}\geq x) \\= e^{-\mu x} \cdot (1-F_{\Gamma}(x|k,\lambda)),$$ with $F_\Gamma$ denoting the c.d.f. of the gamma distribution. And then one could generate draws e.g. via slice-sampling.

Here some additional R code, which relies on the slice-sampling method from the diversitree package:

    require("diversitree")
    loglik <- function(x) -mu*x + pgamma(x, k, lambda, lower=F, log.p=T)
    x2 <- mcmc(loglik, x.init = k/lambda, nsteps = n/100, 
               w = k/(2*lambda), print.every = n/1000, lower = 0)$pars
    mean(x2) # 1.34