...And Then Three Came Along At Once!

| 0 Comments

In the last few posts[1][2] we have been looking at the properties of memoryless processes, being those processes that trigger events at intervals that are independent of how long we have been waiting for one to happen. First we answered the question of what is the probability that we must wait some given time that an event will occur with the exponential distribution, and then the question of what is the probability that we must wait some given time for the \(k^{\mathrm{th}}\) event to occur with the gamma distribution, which allows the somewhat counter-intuitive case of non-integer \(k\).

This time I'd like to ask another question; what is the probability that we'll observe \(k\) events, occurring at a rate \(\lambda\), in a period of one unit of time?

The Poisson Distribution

It turns out that the easiest way to answer this question is to use the gamma CDF
\[ \Gamma_{k,\lambda}(x) = \begin{cases} \frac{\gamma(k, \lambda x)}{\Gamma(k)} & x \geqslant 0\\ 0 & x < 0 \end{cases} \]
to determine the probability that the \(k+1^\mathrm{th}\) event arrives before the period ends
\[ \Gamma_{k+1,\lambda}(1) = \frac{\gamma(k+1, \lambda)}{\Gamma(k+1)} \]
where \(\gamma(k, \lambda)\) is the lower incomplete gamma function and \(\Gamma(k)\) is the gamma function. Note that since \(k\) is an integer
\[ \Gamma(k+1) = k! \]
where the exclamation mark stands for the factorial. If we denote the number of observed events with \(N\) then this means that
\[ \begin{align*} \Pr(N > k) &= \frac{\gamma(k+1, \lambda)}{\Gamma(k+1)}\\ \Pr(N \leqslant k) &= 1 - \frac{\gamma(k+1, \lambda)}{\Gamma(k+1)} = \frac{\Gamma(k+1) - \gamma(k+1, \lambda)}{\Gamma(k+1)} = \frac{\Gamma(k+1, \lambda)}{\Gamma(k+1)} = P_\lambda(k) \end{align*} \]
where \(\Gamma(k, \lambda)\) is the upper incomplete gamma function. \(P_\lambda(k)\) is the CDF of the distribution that governs the number of events, the Poisson distribution \(P(\lambda)\), and is simply the regularised upper incomplete gamma function which we implemented with ak.gammaQ[3].

The probability that the number of events equals \(k\) is given by the probability mass function, which we can easily derive from the CDF
\[ p_\lambda(k) = P_\lambda(k) - P_\lambda(k-1) = \frac{\Gamma(k+1, \lambda)}{\Gamma(k+1)} - \frac{\Gamma(k, \lambda)}{\Gamma(k)} \]
Recall that for the gamma function and the upper incomplete gamma function we have
\[ \begin{align*} \Gamma(k+1) &= k \Gamma(k)\\ \Gamma(k+1, \lambda) &= k \Gamma(k, \lambda) + \lambda^k e^{-\lambda} \end{align*} \]
and so
\[ p_\lambda(k) = \frac{k \Gamma(k, \lambda) + \lambda^k e^{-\lambda}}{k \Gamma(k)} - \frac{\Gamma(k, \lambda)}{\Gamma(k)} = \frac{\lambda^k e^{-\lambda}}{k \Gamma(k)} = \frac{\lambda^k e^{-\lambda}}{k!} \]
Processes that trigger memoryless events are known as Poisson processes since the number of arrivals in a given period is an extremely important quantity in many of the real world waiting time problems that we might model with such events.

Another Day, Another k

The next obvious question is what is the distribution of the number of events that might occur in \(n\) units of time. Given that \(\lambda\) is the rate per unit time, we should rather hope that it is \(P(n\lambda)\); an average of \(\lambda\) events per second should imply an average of \(60\lambda\) events per minute, after all.

To show that the Poisson distribution satisfies this requirement, we first note that the number of events in the first unit of time is independent of the number of events in the second is independent of the number in the third, etcetera, etcetera. The number of events that occur in \(n\) units of time is therefore simply the sum of the numbers of events that occurred in each of them.
Derivation 1: The Poisson CF
Rearranging the terms yields
\[ \begin{align*} \hat{p}_\lambda(t) &= e^{-\lambda} \sum_{k=0}^\infty e^{itk} \frac{e^{k \ln \lambda}}{k!}\\ &= e^{-\lambda} \sum_{k=0}^\infty \frac{1}{k!} e^{(it + \ln \lambda)k}\\ &= e^{-\lambda} \sum_{k=0}^\infty \frac{1}{k!} \left(e^{it + \ln \lambda}\right)^k \end{align*} \]
Now the sum is simply the series expansion of an exponential; specifically
\[ \sum_{k=0}^\infty \frac{1}{k!} \left(e^{it + \ln \lambda}\right)^k = \exp\left(e^{it + \ln \lambda} \right) \]
Rearranging again we have
\[ \begin{align*} \hat{p}_\lambda(t) &= e^{-\lambda} \exp\left(e^{it} e^{\ln \lambda}\right)\\ &= \exp\left(\lambda e^{it} - \lambda\right)\\ &= \exp\left(\lambda \left(e^{it} - 1\right)\right) \end{align*} \]
As usual when reasoning about sums of random variables, we shall turn to the characteristic function since the CF of the sum of a set of random variables is simply the product of their individual CFs. Now since Poisson random variables are discrete, rather than continuous, we use sums, rather than integrals, to calculate expectations and so the CF is given by
\[ \hat{p}_\lambda(t) = \mathrm{E}\left[e^{itN}\right] = \sum_{k=0}^\infty e^{itk} \frac{\lambda^k e^{-\lambda}}{k!} \]
where \(\sum\) is the summation sign. Whilst this may look like it might be rather time consuming to calculate, derivation 1 shows that it can be simplified to
\[ \hat{p}_\lambda(t) = \exp\left(\lambda \left(e^{it} - 1\right)\right) \]
where \(\exp(x)\) is another way of writing \(e^x\).
The CF of the sum of \(n\) observations of a Poisson random variate is therefore
\[ \begin{align*} \hat{p}_\lambda(t)^n &= \exp\left(\lambda \left(e^{it} - 1\right)\right)^n\\ &= \exp\left(n\lambda \left(e^{it} - 1\right)\right)\\ &= \hat{p}_{n\lambda}(t) \end{align*} \]
and so the number of events in \(n\) units of time are distributed as \(P(n\lambda)\), as hoped for!
Note that we can trivially extend this to fractional units of time by supposing that the number of events in \(\frac{1}{n}\) units of time is distributed as \(P(\lambda^\prime)\) and so, considering the total number of events in \(n\) of these smaller periods, we have
\[ \begin{align*} P(n\lambda^\prime) &= P(\lambda)\\ n\lambda^\prime &= \lambda\\ \lambda^\prime &= \frac{\lambda}{n} \end{align*} \]
Finally, extending the rule to periods of irrational length simply requires a limit argument on sequences of periods of rational lengths and we therefore have a general rule that if the number of events in unit time is distributed as \(P(\lambda)\) then the number of events in time \(T\) is distributed as \(P(T\lambda)\).

Implementing The PMF, The CDF And The CF

Given that the PMF only has mass for non-negative integer \(k\), it's tempting to round its argument to an integer to compensate for possible rounding errors. This could, however, lead to errors in which fractional numbers of events are treated as if they have non-zero probabilities; a nonsensical result that should be avoided at all costs. We shall therefore return zero for such arguments, as shown in listing 1, and leave it to the user to explicitly deal with rounding errors should they be an issue.

Listing 1: ak.poissonPMF
function pmf(l, ll, k) {
 if(k>=0 && k===ak.floor(k)) return Math.exp(k*ll - l - ak.logGamma(k+1));
 return isNaN(k) ? ak.NaN : 0;
}

ak.poissonPMF = function() {
 var state = {lambda: 1};
 var arg0  = arguments[0];
 var f, ll;

 constructors[ak.nativeType(arg0)](state, arg0, arguments);

 ll = Math.log(state.lambda);
 f = function(k){return pmf(state.lambda, ll, k);};
 f.lambda = function(){return state.lambda;};
 return Object.freeze(f);
};

var constructors = {};

Note that we're exponentiating the logarithm of the PMF to stave off overflow in the numerator and the denominator in the pmf helper function, using ak.logGamma to compute the logarithm of the factorial of \(k\).
We again defer the initialisation of the state object to a constructors object whose relevant members are given in listing 2.

Listing 2: The constructors Object
constructors[ak.UNDEFINED_T] = function() {
};

constructors[ak.FUNCTION_T] = function(state, rnd) {
 state.rnd = rnd;
};

constructors[ak.NUMBER_T] = function(state, lambda, args) {
 var arg1 = args[1];
 constructors[ak.NUMBER_T][ak.nativeType(arg1)](state, arg1, args);

 state.lambda = Number(lambda);
 if(state.lambda<=0 || !isFinite(state.lambda)) {
  throw new Error('invalid lambda in ak.poisson distribution');
 }
};

constructors[ak.NUMBER_T][ak.UNDEFINED_T] = function(state) {
};

Program 1 illustrates how the shape of the Poisson PMF changes with \(\lambda\).

Program 1: The Poisson PMF

Now unlike the PMF, we can safely round the argument of the CDF down to an integer, as is done in listing 3, since the number of events cannot be fractional and so the probability that there are less than or equal to a non-integer \(k\) of them is trivially equal to the probability that they number less than or equal to its integer part.

Listing 3: ak.poissonCDF
ak.poissonCDF = function() {
 var state = {lambda: 1};
 var arg0  = arguments[0];
 var f;

 constructors[ak.nativeType(arg0)](state, arg0, arguments);

 f = function(k){return k<0 ? 0 : ak.gammaQ(ak.floor(k+1), state.lambda);};
 f.lambda = function(){return state.lambda;};
 return Object.freeze(f);
};

Next we have the implementation of the CF, given in listing 4.

Listing 4: ak.poissonCF
ak.poissonCF = function() {
 var state = {lambda: 1};
 var arg0  = arguments[0];
 var f;

 constructors[ak.nativeType(arg0)](state, arg0, arguments);

 f = function(t) {
  return ak.exp(ak.complex(state.lambda*(Math.cos(t)-1), state.lambda*Math.sin(t)));
 };
 f.lambda = function(){return state.lambda;};
 return Object.freeze(f);
};

Note that we're slightly optimising its calculation by exploiting the fact that the exponent can be expressed as
\[ \lambda \left(e^{it}-1\right) = \lambda \left(\cos t - 1\right) + i \lambda \sin t \]

Random Numbers

The easiest way to generate Poisson distributed random numbers is to note that they are numbers of events with exponentially distributed waiting times that occur in one unit of time, so we can simply count how many observations of an exponential variate we need to make before their sum exceeds one, as illustrated by listing 5.

Listing 5: A Simple Poisson Random Number Generator
function rnd(lambda) {
 var expRnd = ak.exponentialRnd(lambda);
 var k = 0;
 var t = expRnd();

 while(t<=1) {
  t += expRnd();
  ++k;
 }
 return k;
}

The problem with this approach is that, for large \(\lambda\), we'll need to generate a lot of exponentially distributed random numbers before the loop terminates, seriously impacting its efficiency.
A better scheme is to divide the interval \([0,1)\) into an infinite set of non-overlapping intervals \(I_k\) with widths equal to the probability of observing \(k\) events, \(p_\lambda(k)\). If an observation of a standard uniform distributed random number \(U\) falls within \(I_k\) then \(k\) must therefore be an observation of a Poisson distributed random number with rate \(\lambda\) since
\[ 0 \leqslant a \leqslant b \leqslant 1 \implies \Pr(U \in [a,b]) = b-a \]
where \(\implies\) means implies and \(\in\) means within, and so
\[ \Pr\left(U \in I_k\right) = p_\lambda(k) \]
The easiest way to find the relevant interval is to add up values of \(p_\lambda(k)\) until their sum first equals or exceeds \(U\), implying that it falls somewhere within \(I_k\), as demonstrated by listing 6.

Listing 6: A Better Poisson Random Number Generator
function rnd(lambda) {
 var u = Math.random();
 var pmf = ak.poissonPMF(lambda);
 var k = 0;
 var sum = pmf(k);

 while(sum<u) sum += pmf(++k);
 return k;
}

This brings the number of random numbers we need per observation of a Poisson random variable down to just one, which is definitely a step in the right direction. Unfortunately, for large \(\lambda\), we'll instead need lots of calls to ak.poissonPMF, which definitely isn't.

We can improve matters still further by noting that finding the value for which the running sum of the PMF first equals or exceeds a given value is essentially inverting the CDF of a discrete distribution. We have already seen how we can reasonably efficiently invert a continuous function with the bisection method[4], repeatedly splitting an interval containing the sought after value into two parts until we have one whose width is smaller than some specified precision.
To apply this to a discrete CDF, we must round the splitting point down to an integer and keep the interval \((a_i, b_i]\) that satisfies
\[ P_\lambda(a_i) < U \leqslant P_\lambda(b_i) \]
As soon as we have \(b_i-a_i \leqslant 1\), the value we seek must be \(b_i\).

Since we're going to want the inverse of the CDF anyway, we shall first use this algorithm in ak.poissonInvCDF, as defined in listing 7.

Listing 7: ak.poissonInvCDF
ak.poissonInvCDF = function() {
 var state = {lambda: 1};
 var arg0  = arguments[0];
 var b1, g0, g1, f;

 constructors[ak.nativeType(arg0)](state, arg0, arguments);

 b1 = ak.ceil(state.lambda);
 g0 = ak.gammaQ(1, state.lambda);
 g1 = ak.gammaQ(b1+1, state.lambda);

 f = function(c){return invCDF(state.lambda, b1, g0, g1, c);};
 f.lambda = function(){return state.lambda;};
 return Object.freeze(f);
};

After initialising the state with the constructors object, we assume an initial guess for the bracket of \((0, b_1]\) whose CDF values are stored in g0 and g1 respectively and are passed to the helper function invCDF, given in listing 8.

Listing 8: invCDF
function invCDF(l, b1, g0, g1, c) {
 var b0 = 0;
 var m;

 if(isNaN(c)) return ak.NaN;
 if(c>=1)     return ak.INFINITY;
 if(c<=g0)    return b0;

 while(g1<c) {
  b0 = b1;
  b1 *= 2;
  g1 = ak.gammaQ(b1+1,l);
 }
 m = ak.floor(b0/2 + b1/2);

 while(m!==b0 && m!==b1) {
  if(ak.gammaQ(m+1,l)<c) b0 = m;
  else                   b1 = m;
  m = ak.floor(b0/2 + b1/2);
 }
 return b1;
}

Note that we handle the corner cases of a NaN argument and infinite or zero results first since they would violate the assumption that the result is greater than the lower bound of the initial interval and less than or equal to its upper bound.
To find that interval, we repeatedly double the upper bound of our first guess until it contains the required value of the CDF, moving the lower bound each time to rule out arguments that we already know can't yield that value. The termination condition of the subsequent loop is that the midpoint of the interval is equal to one of the end points, rather than its endpoints being separated by less than one, so that we don't get bitten by rounding errors for very large values of b0 and b1.

Now, given that we're seeking an integer result, this will terminate much faster than the numerical inverse of a continuous CDF, which we should typically want to converge to an interval of width much smaller than one, but it's still rather computationally inefficient for very large \(\lambda\). What we really need is a rejection method like the ones we used last time for the gamma distribution. Efficient schemes exist[5], but they are rather complicated and it strikes me that it would be instructive to derive a simple, if somewhat sub-optimal, rejection algorithm from scratch.

(I Want To Be) Rejected

Recall that for rejection sampling we draw a random number, \(x\), from a distribution for which we can easily generate them and only keep it if
\[ u \leqslant \frac{p(x)}{cq(x)} \]
where \(u\) is a uniform random number, \(p(x)\) is the PDF of the distribution we want to draw from, \(q(x)\) is the PDF of the distribution we can easily draw from and \(c\) satisfies
\[ \frac{p(x)}{q(x)} \leqslant c \]
for all \(x\) at which either PDF is non-zero.

Specifically, we shall seek to generate Poisson distributed random numbers by rejecting Cauchy distributed random numbers, having the PDF
\[ p_{\mu, \sigma}(x) = \left(\pi\sigma \left(1 + \left(\frac{x-\mu}{\sigma}\right)^2\right)\right)^{-1} \]
This has a number of convenient properties for rejection sampling. Firstly, it falls away to zero much more slowly than most distributions and so the ratio of the PDF of the distribution we want to draw from and its will usually tend to zero for \(x\) far from \(\mu\). Secondly, the CDF is given by
\[ P_{\mu, \sigma}(x) = \frac{1}{\pi}\arctan\left(\frac{x-\mu}{\sigma}\right)+\frac{1}{2} \]
which is trivial to invert and so we can easily generate random numbers from it with
\[ x = \mu + \sigma \tan\left(\left(u-\frac{1}{2}\right)\pi\right) \]
Since the Poisson distribution is discrete rather than continuous, we must instead find \(c\) such that
\[ \frac{p_\lambda(k)}{P_{\mu,\sigma}(k+1)-P_{\mu,\sigma}(k)} \leqslant c \]
for all non-negative integer \(k\), where \(p_\lambda\) is the Poisson PMF and \(P_{\mu,\sigma}\) is the Cauchy CDF, round down the Cauchy distributed random numbers to integer \(k\) and accept them if
\[ k \geqslant 0 \wedge u \leqslant \frac{p_\lambda(k)}{c\left(P_{\mu,\sigma}(k+1)-P_{\mu,\sigma}(k)\right)} \]
where \(\wedge\) means and.

Figure 1: The Ratio Of Masses For Each k
Now the Poisson PMF peaks near \(\lambda\) and the Cauchy PDF peaks at \(\mu\), so it makes sense to set \(\mu\) equal to \(\lambda\). Choosing an appropriate value for \(\sigma\) is rather more difficult, but extensive testing has strongly suggested to me that \(\sqrt\lambda\) is a reasonable one.
The mass of the Cauchy distribution initially falls towards zero faster than that of the Poisson distribution before levelling off as we get further from its peak. As a result, the ratio of their masses has two local maxima, located to the left and right of \(\lambda\), as illustrated by figure 1.

To my regret, I have not been able to prove that these maxima are located near \(\lambda-\sqrt\lambda\) and \(\lambda+\sqrt\lambda\), nor that the rightmost is consistently greater than the leftmost, as is the case in this example, but the evidence from my tests has been sufficient to convince me that it is true and hopefully program 2 will convince you too.

Program 2: Maxima Of The Rejection Constant

Do try changing lambda to both integer and non-integer values of varying magnitudes to see the effect on values of the ratio in the neighbourhoods of k1 and k2. Note that, on average, we'll need \(c\) Cauchy distributed random numbers per Poisson distributed random number, so this will give you a sense of how efficient this scheme will be.

Now despite my confidence in the claim that the rightmost peak is always the greater, we shall err on the side of caution and numerically search for the maximum starting from both of these points which, given my confidence that they are both always very close to the local maxima, shouldn't prove too costly.

Putting It All Together

Derivation 2: From A Sum To A Product
Noting that \(x_i = -\frac{\ln \left(1-u_i\right)}{\lambda}\) and that \(1-u_i\) is uniformly distributed, we have
\[ \begin{align*} \sum_{i=1}^k x_i \leqslant 1 &\implies \sum_{i=1}^k -\frac{\ln (1-u_i)}{\lambda} \leqslant 1\\ &\implies \sum_{i=1}^k \ln u_i \geqslant -\lambda\\ &\implies \ln \prod_{i=1}^k u_i \geqslant -\lambda\\ &\implies \prod_{i=1}^k u_i \geqslant e^{-\lambda}\\ \end{align*} \]
We can significantly improve the efficiency of random number generation for small \(\lambda\) by noting that since the waiting time for events is memoryless we can simply keep track of how much time there is to wait for the next one.
Furthermore, we can rearrange the test against the sum of exponentially distributed random numbers, \(x_i\), into a test against a product of uniformly distributed random numbers, \(u_i\)
\[ \begin{align*} \sum_{i=1}^k x_i &\leqslant 1\\ \prod_{i=1}^k u_i &\geqslant e^{-\lambda} \end{align*} \]
where \(\prod\) is the product sign, as shown in derivation 2.

Both of these observations are exploited by expRnd in listing 9, which we shall use when \(\lambda\) is less than or equal to two.

Listing 9: expRnd
function expRnd(state, el) {
 var k = 0;

 if(state.exp===0) state.exp = state.rnd();

 while(state.exp>=el) {
  state.exp *= state.rnd();
  ++k;
 }
 state.exp /= el;

 return k;
}

Note that el will be equal to \(e^{-\lambda}\) and state.exp is the rearranged sum of exponentials; its division by el after the loop being equivalent to subtracting one from the original sum to account for a unit of time having passed.

For \(\lambda\) up to eight we shall simply apply the inverse CDF to a uniform random number, so all that remains to do is to implement the rejection scheme. Listing 10 gives the implementation of the ratio of masses that we shall use to determine when to reject a rounded down Cauchy distributed random number, \(k\)

Listing 10: The Rejection Ratio
function reject(l, ll, rl, k) {
 var pp = pmf(l, ll, k);
 var pc = (Math.atan((k+1-l)/rl)-Math.atan((k-l)/rl))/ak.PI;
 return pp/pc;
}

Here l will equal \(\lambda\), ll will equal \(\ln \lambda\) and rl will equal \(\sqrt\lambda\). The rejection algorithm is easy to implement using this function, as shown by listing 11.

Listing 11: rejectRnd
function rejectRnd(l, ll, rl, c, rnd) {
 var k;

 do k = ak.floor(l + rl*Math.tan((rnd()-0.5)*ak.PI));
 while(k<0 || c*rnd()>reject(l, ll, rl, k));

 return k;
}

Note that the loop simply keeps assigning a rounded down Cauchy distributed random number to k until it passes the acceptance test. To find a candidate for the rejection constant we simply iterate over \(k\) from an initial guess until we find the local maximum, as shown by listing 12.

Listing 12: Finding The Maximum Rejection Ratio
function maxReject(l, ll, rl, k0) {
 var k1 = k0+1;
 var c0 = reject(l, ll, rl, k0);
 var c1 = reject(l, ll, rl, k1);
 var c;

 if(c0>c1) {
  while(c0>c1) {
   --k1; c1 = c0;
   --k0; c0 = reject(l, ll, rl, k0);
  }
  c = c1;
 }
 else {
  while(c0<c1) {
   ++k0; c0 = c1;
   ++k1; c1 = reject(l, ll, rl, k1);
  }
  c = c0;
 }
 return c;
}

The direction in which we iterate is decided by which of the rejection ratios, c0 and c1, is the greater.

Finally, we select which of these schemes we shall use according to the value of \(\lambda\) in ak.poissonRnd, given in listing 13.

Listing 13: ak.poissonRnd
ak.poissonRnd = function() {
 var state = {lambda: 1, rnd: Math.random};
 var arg0  = arguments[0];
 var b1, g0, g1, el, ll, rl, c, f;

 constructors[ak.nativeType(arg0)](state, arg0, arguments);

 if(state.lambda<=2) {
  state.exp = 0;
  el = Math.exp(-state.lambda);

  f = function(){return expRnd(state, el);};
 }
 else if(state.lambda<=8) {
  b1 = ak.ceil(state.lambda);
  g0 = ak.gammaQ(1, state.lambda);
  g1 = ak.gammaQ(b1+1, state.lambda);

  f = function(){return invCDF(state.lambda, b1, g0, g1, state.rnd());};
 }
 else {
  ll = Math.log(state.lambda);
  rl = Math.sqrt(state.lambda);
  c  = maxReject(state.lambda, ll, rl, ak.floor(state.lambda-rl));
  c  = Math.max(c, maxReject(state.lambda, ll, rl, ak.floor(state.lambda+rl)));

  f = function(){return rejectRnd(state.lambda, ll, rl, c, state.rnd);};
 }

 f.lambda = function(){return state.lambda;};
 f.rnd = function(){return state.rnd;};
 return Object.freeze(f);
};

constructors[ak.FUNCTION_T] = function(state, rnd) {
 state.rnd = rnd;
};

constructors[ak.NUMBER_T][ak.FUNCTION_T] = function(state, rnd) {
 state.rnd = rnd;
};

Note that we set up the constants that we need for each approach and pass them to the implementations of f in the bodies of the if statements that select which one to use. As usual, we provide additional constructors to allow the user to use their own uniform random number generator should they not wish to use the default, Math.random.

Program 3 compares a bar chart of the results of calls to ak.poissonRnd with a graph of the Poisson PMF, with \(\lambda\) chosen to use our own rejection algorithm.

Program 3: ak.poissonRnd

Yet another good match, I'm sure you'll agree, but you should of course run this program with other values of \(\lambda\) to test the inversion and summation random number generators.

We're almost through with the properties of memoryless processes, but until next time all of these functions can be found in PoissonDistribution.js.

References

[1] How Long Must I Wait?, www.thusspakeak.com, 2015.

[2] kth Time's The Charm, www.thusspakeak.com, 2015.

[3] Gamma Gamma Hey!, www.thusspakeak.com, 2015.

[4] Rooting Around For Answers, www.thusspakeak.com, 2015.

[5] Devroye, L., Non-Uniform Random Variate Generation, Springer Verlag, 1986.

Leave a comment