How R sample from exponential distribution?

Question

I managed to find the source code in sexp.c, and the algorithm (Ahrens & Dieter). I mostly understand the first half of the code - it seems like it finds the coarse location of the returned value, using the fact that $\int_{k\cdot\ln2}^{(k+1)\cdot\ln2}e^{-x}dx=2^{-k-1}$. So if we stopped there, more or less, we would get a distribution which is a step function of the exp:

I don't understand the 2nd half.

I guess I agree that $u$ (in this point) is again uniform distributed in $U(0, \ln 2)$ (the $\ln2$ comes from multiplying the result). I more or less understand that: "With probabilities ... $(\ln 2)^i/i!$ consider the minimum of $i$ uniform samples from $(0, \ln2)$". I just don't know how you show this is indeed an exponential distribution ... I know how to find the pdf of $\min u_i$ but here it has different probabilities for each $i$.

Can anyone show this?

Here is the code:

double exp_rand(void)
{
    /* q[k-1] = sum(log(2)^k / k!)  k=1,..,n, */
    /* The highest n (here 16) is determined by q[n-1] = 1.0 */
    /* within standard precision */
    const static double q[] =
    {
    0.6931471805599453,
    0.9333736875190459,
    0.9888777961838675,
    0.9984959252914960,
    0.9998292811061389,
    0.9999833164100727,
    0.9999985691438767,
    0.9999998906925558,
    0.9999999924734159,
    0.9999999995283275,
    0.9999999999728814,
    0.9999999999985598,
    0.9999999999999289,
    0.9999999999999968,
    0.9999999999999999,
    1.0000000000000000
       };
    double a = 0.;
    double u = unif_rand();    /* precaution if u = 0 is ever returned */
    while(u <= 0. || u >= 1.) u = unif_rand();     for (;;) {
        u += u;
        if (u > 1.)
            break;
        a += q[0];
        }
        u -= 1.;
    
        if (u <= q[0])
        return a + u;  # Up to here I understand
    
        int i = 0;
        double ustar = unif_rand(), umin = ustar;
        do {
        ustar = unif_rand();
        if (umin > ustar)
            umin = ustar;
        i++;
        } while (u > q[i]);
        return a + umin * q[0];
    }

Answer 1

The answer is (as usual!) provided in Devroyes' Non-uniform random variate generation (1986, p.396). The principle for the Ahrens-Dieter (1972) algorithm is a result due to George Marsaglia (1961):

Theorem IX.2.1$\ \ \ $If $U_1,U_2,\ldots$ is a series of iid $\mathcal U(0,1)$ random variables, if $Z$ is an independent positive Poisson $\mathcal P_+(\mu)$ random variable, and if $M$ is an independent Geometric $\mathcal G(1-e^{-\mu})$ random variable, then $$X = \mu(M+\min(U_1,\ldots,U_Z))\sim\mathcal E(1)$$

The proof proceeds as follows (reproducing from Devroye, p.395): \begin{align*} \mathbb P(\mu\min(U_1,\ldots,U_Z)\le x) &=\mathbb E^Z[\mathbb P(\mu\min(U_1,\ldots,U_z)\le x|Z=z)]\\ &=\mathbb E^Z[1-(1-x/\mu)^Z)]\\ &=1-\dfrac{e^{\mu-x}-1}{e^{\mu}-1}\\ &=\dfrac{1-e^{-x}}{1-e^{-\mu}} \end{align*} which is the cdf of the exponential distribution truncated to $(0,\mu)$ (see here). A preliminary result$^1$ due to Von Neumann (Lemma IV.2., p.125 & p.393) is that, when $$Z\sim \mathcal G(1-e^{-\mu})\qquad Y\sim \dfrac{e^{-y}}{1-e^{-\mu}}\mathbb I_{(0,1)}(y)$$then$$\mu(Z-1)+Y\sim\mathcal E(1)$$ This follows from considering the moment generating function (for $t<\min(\mu,1)$) $$\mathbb E[e^{tX}]=\mathbb E[e^{t\mu(Z-1)}]\mathbb E[e^{tY}]= \frac{1-e^{-\mu}}{1-e^{-\mu(1-t)}}\frac{1-e^{-\mu(1-t)}}{(1-e^{-\mu)}(1-t)}=\frac{1}{1-t}$$ which concludes the proof.

Relating to the graph included in the question, this means that $\mu M$ corresponds to the area under the step function, while $\mu\min(U_1,\ldots,U_Z))$ corresponds to the residual area between the step function and the Exponential $\mathcal E(1)$ density, which is therefore independent from $M$.

Marsaglia then derives his Exponential algorithm from Theorem IX.2.1:

1. Generate a Geometric $M\sim\mathcal G(1-e^{-\mu})$ variable
2. Generate two Uniform $U$ and $V$
3. Set $Y=V$ and $Z=1$
4. While $U>F_\mu(Z)$, increase $Z$ to $Z+1$ and decrease $Y$ to $\min(Y,W)$, where $W$ is Uniform
5. Return $\mu(M+Y)$ as an Exponential $\mathcal E(1)$ variate.

and Ahrens and Dieter (1972) consists in an optimisation of the above, for instance by choosing$^2$ $\mu=\log(2)$, refining the generation of $M$, and storing the cdf $F(\cdot)$.

A detailed explanation of Ahrens and Dieter (1972) version:

1. The prerecorded table corresponds to the first terms of the cdf of the Poisson $\mathcal P_+(\log 2)$ distribution (with q[0] equal to $\mu$)
2. The first loop returns a rescaled Geometric $\mathcal G(1-e^{-\mu})$ variate $\mu M$ as a (using a sequential search inversion method and the property that $F_M(i)=1-2^{-i}$ when $\mu=\log 2$)
3. The residual u-1 produces$^3$ an independent uniform variate $V$
4. The case u <= q[0] corresponds to $Z=1$ and avoids running the second loop
5. The second loop while (u > q[i]) produces the Poisson $\mathcal P_+(\log 2)$ variate $Z$ as i-1and the associated $Y=\min(U_1,\ldots,U_Z)$ as umin
6. The outcome a+umin*q[0] is indeed $X = \mu(M+\min(U_1,\ldots,U_Z))\sim\mathcal E(1)$

---

$^1$An ingenious algorithm for generating from the Exponential distribution is derived from this lemma and consists in only producing sequences of Uniform variates $U_0,U_1,\ldots$.

$^2$When $\mu=\log 2$, the geometric random variate corresponds to the number of $0$ before the first $1$ in the binary expansion of $U\sim\mathcal U(0,1)$. Sampling directly these bits proves to be much faster than the first loop until (u > 1.).

$^3$Devroye remarks that "Ahrens and Dieter squeeze the first uniform [0,1] random variate $U$ dry". The efficiency of the method is such that it requires on average $1+\log(2)$ uniforms.