4
$\begingroup$

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: enter image description here

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];
}
$\endgroup$
7
$\begingroup$

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.

$\endgroup$
6
  • $\begingroup$ Obviously... But can you maybe elaborate for those it's not so obvious to them, like me :-) How do you show this? CDF? Forget the $\mu$ for a second, how do you show that $M+min(U_1,...U_Z) \sim Exp(1)$ ? $\endgroup$ May 1 at 19:43
  • $\begingroup$ ah ok, I didn't know about this bible, thanks! $\endgroup$ May 1 at 20:43
  • $\begingroup$ I can't find this in the link, and my university doesn't have access to full book. Also, maybe I got something wrong, but my calculations are: $\mathbb E^Z[1-(1-x)^Z)] = 1-\sum_{z=0}^\infty (1-x)^z \frac{e^{-1}}{z!} = 1-e^{-1}e^1e^{-x} = 1-e^{-x}$ $\endgroup$ May 1 at 21:06
  • $\begingroup$ Oh I got it, Z is from 1 up, without 0. $\endgroup$ May 1 at 21:29
  • $\begingroup$ Any chance of posting the Lemma IV 2? $\endgroup$ May 1 at 21:34

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.