# generate random number from exponential distribution with random right truncation

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


## 2 Answers

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  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. • Could you elaborate further? How could I then generate random numbers efficiently so that they do follow the distribution of$x1$? – mplatzer Nov 2 '14 at 17:20 • or were you saying that the simulation$x1$is flawed, and the$x0\$ draws are correct? – mplatzer Nov 2 '14 at 17:49
• I am saying that your method for x1 is incorrect. If you want to reject some of them, then to avoid bias you need to replace the with others for the same trunc_value. – Henry Nov 2 '14 at 21:47
• I now found a solution how to draw in order to match x1. see the answer which I just posted. – mplatzer Nov 3 '14 at 13:10