# How r sample from exponential distribution?

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];
}


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.

• 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)$ ? May 1 at 19:43
• ah ok, I didn't know about this bible, thanks! May 1 at 20:43
• 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}$ May 1 at 21:06
• Oh I got it, Z is from 1 up, without 0. May 1 at 21:29
• Any chance of posting the Lemma IV 2? May 1 at 21:34