1
$\begingroup$

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
$\endgroup$
1
$\begingroup$

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
$\endgroup$
0
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Could you elaborate further? How could I then generate random numbers efficiently so that they do follow the distribution of $x1$? $\endgroup$ – mplatzer Nov 2 '14 at 17:20
  • $\begingroup$ or were you saying that the simulation $x1$ is flawed, and the $x0$ draws are correct? $\endgroup$ – mplatzer Nov 2 '14 at 17:49
  • $\begingroup$ 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. $\endgroup$ – Henry Nov 2 '14 at 21:47
  • $\begingroup$ I now found a solution how to draw in order to match x1. see the answer which I just posted. $\endgroup$ – mplatzer Nov 3 '14 at 13:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.