Last time^{[1]} we took a first look at Bernoulli processes which are formed from a sequence of independent experiments, known as Bernoulli trials, each of which is governed by the Bernoulli distribution with a probability \(p\) of success. Since the outcome of one trial has no effect upon the next, such processes are memoryless meaning that the number of trials that we need to perform before getting a success is independent of how many we have already performed whilst waiting for one.
We have already seen^{[2]} that if waiting times for memoryless events with fixed average arrival rates are continuous then they must be exponentially distributed and in this post we shall be looking at the discrete analogue.
because of which the distribution of \(K\) is called the geometric distribution. The inverse of the CDF is given by
To make observations of \(K\) we can simply evaluate the inverse at the result of a standard uniformly distributed random variable
Note that the geometric and exponential distributions are the only memoryless waiting time distributions and that the latter is the limit of the former as the probability of success tends to zero and the number of trials tends to infinity, as shown by derivation 2.
We have previously seen that the minimum of \(n\) independent observations of an exponential random variable is also exponentially distributed^{[3]}, which is another property that it shares with the geometric distribution, as proven in derivation 3.
Here the
Note that, as we did for the Bernoulli distribution, we're constraining the probability of success to lie strictly between zero and one.
Program 1 plots the PMF for a range of success probabilities.
The CDF is even simpler, as shown by listing 3.
Program 2 demonstrates its correctness by comparing the PMF at \(k\) to the difference between the CDF at \(k\) and at \(k1\).
In order to cope with rounding errors, when calculating the inverse of the CDF we round down instead of up and then increment the result until the value of the CDF is no less than its argument, as done by the
Note that
Program 3 compares the inverse of the CDF's result to its argument.
To calculate the CF we shall exploit the facts that
and program 4 plots the real and imaginary parts of its results in black and red respectively.
Finally we have the random number generator for which we need constructors to allow the specification of standard uniform random number generators other than the default
Here we're simply calling the helper function
Program 5 illustrates the use of
You will find the implementations of these functions in GeometricDistribution.js and we shall take a look at another property of Bernoulli processes next time.
[2] How Long Must I Wait?, www.thusspakeak.com, 2015.
[3] The Longest Wait, www.thusspakeak.com, 2015.
We have already seen^{[2]} that if waiting times for memoryless events with fixed average arrival rates are continuous then they must be exponentially distributed and in this post we shall be looking at the discrete analogue.
A Series Of Unfortunate Events
Since the trials are independent, the probability of having \(k\) failures in a row is trivially
\[
(1p)^k
\]
If we represent the number of failures before a success with a random variable \(K\) then its probability mass function, or PMF, must be
\[
f_p(k) = \Pr(K = k) = (1p)^k \times p
\]
for nonnegative integer \(k\) and zero otherwise. To calculate its cumulative distribution function, or CDF, we simply add up values of the PMF with
\[
\begin{align*}
F_p(k) &= \Pr(K \leqslant k) = \Pr\left(K \leqslant \lfloor k \rfloor\right)\\
&= \sum_{j=0}^{\lfloor k \rfloor} f_p(j)
= \sum_{j=0}^{\lfloor k\rfloor} (1p)^j \times p
\end{align*}
\]
where the odd shaped brackets represent the floor of the value between them, being the largest integer that is less than or equal to it, and \(\sum\) is the summation sign. Now this is simply a geometric series and so we have
\[
F_p(k) =
\begin{cases}
\dfrac{p \times \left(1(1p)^{\lfloor k \rfloor +1}\right)}{1  (1p)} = 1  (1p)^{\lfloor k \rfloor + 1} & k \geqslant 0\\
0 & \text{otherwise}
\end{cases}\\
\]
Derivation 1: The Inverse CDF
Firstly, define
\[
\begin{align*}
c &= 1  (1p)^{k+1}\\
1c &= (1p)^{k+1}
\end{align*}
\]
Taking logarithms yields
\[
\begin{align*}
\ln (1c) &= (k+1) \times \ln (1p)\\
k &= \frac{\ln (1c)}{\ln (1p)}  1
\end{align*}
\]
The result follows from the fact that if \(c\) is greater than \(F_p\left(\lfloor k \rfloor\right)\) then \(k\) will be greater than \(\lfloor k \rfloor\) and therefore within the accumulated mass at \(\lceil k \rceil\).
\[
F_p^{1}(c) = \Bigg\lceil \frac{\ln (1c)}{\ln (1p)}  1 \Bigg\rceil
\]
where these odd shaped brackets represent the ceiling of the value between them, being the smallest integer that is greater than or equal to it, as shown by derivation 1.To make observations of \(K\) we can simply evaluate the inverse at the result of a standard uniformly distributed random variable
\[
\begin{align*}
U &\sim U(0, 1)\\
K &= F_p^{1}(U)
\end{align*}
\]
Finally, we have the characteristic function, or CF, which is defined by
\[
\hat{f_p}(t) = E\left[e^{itK}\right]
\]
where \(E\) represents the expected value of the term within the square brackets following it and \(i\) is the imaginary unit, and equals
Derivation 2: The Limit
Defining \(p = \frac{\lambda}{n}\) and \(X = \frac{K+1}{n}\) we have
\[
\begin{align*}
\Pr(X \leqslant x) &= \Pr(K \leqslant n \times x  1)\\
&= 1  \left(1\frac{\lambda}{n}\right)^{\lfloor n \times x  1\rfloor + 1}
\end{align*}
\]
Taking the limit as \(n\) tends to infinity yields
\[
\begin{align*}
\lim_{n \rightarrow \infty}\Pr(X \leqslant x)
&= \lim_{n \rightarrow \infty} \left(1  \left(1\frac{\lambda}{n}\right)^{n \times x}\right)\\
&= 1  \left(\lim_{n \rightarrow \infty} \left(1\frac{\lambda}{n}\right)^n\right)^x\\
&= 1  \left(e^{\lambda}\right)^x
= 1  e^{\lambda \times x}
\end{align*}
\]
which is the exponential CDF.
\[
\begin{align*}
\hat{f_p}(t) &= \sum_{k=0}^\infty (1p)^k \times p \times e^{itk}\\
&= \sum_{k=0}^\infty p \times \left((1p)\times e^{it}\right)^k
\end{align*}
\]
This is another geometric series and so we have
\[
\begin{align*}
\hat{f_p}(t) &= p \times \frac{1\left((1p) \times e^{it}\right)^\infty}{1(1p) \times e^{it}}\\
&= \frac{p}{1(1p) \times e^{it}}
\end{align*}
\]
provided that \(1p\) is less than one.Note that the geometric and exponential distributions are the only memoryless waiting time distributions and that the latter is the limit of the former as the probability of success tends to zero and the number of trials tends to infinity, as shown by derivation 2.
We have previously seen that the minimum of \(n\) independent observations of an exponential random variable is also exponentially distributed^{[3]}, which is another property that it shares with the geometric distribution, as proven in derivation 3.
Derivation 3: The Minimum
Firstly note that
\[
\Pr(K > k) = 1  \Pr(K \leqslant k) = (1p)^{\lfloor k \rfloor + 1}
\]
and secondly note that if the minimum of the \(K_i\) is greater than \(k\) then they must all be greater than \(k\) so that
\[
\begin{align*}
\Pr\left(\min_{1 \leqslant i \leqslant n} K_i > k\right)
&= \Pr\left(K_1 > k \wedge K_2 > k \wedge \dots \wedge K_n >k\right)\\
&= \Pr\left(K_1 > k\right) \times \Pr\left(K_2 > k\right) \times \dots \times \Pr\left(K_n >k\right)\\
&= (1p)^{\lfloor k \rfloor + 1} \times (1p)^{\lfloor k \rfloor + 1} \times \dots \times (1p)^{\lfloor k \rfloor + 1}\\
&= \left((1p)^{\lfloor k \rfloor + 1}\right)^n
= \left((1p)^{n}\right)^{\,\lfloor k \rfloor + 1}
\end{align*}
\]
where \(\wedge\) means and, and consequently
\[
\Pr\left(\min_{1 \leqslant i \leqslant n} K_i \leqslant k\right)
= 1  \left(1p^\prime\right)^{\,\lfloor k \rfloor + 1}
= F_{p^\prime}(k)
\]
where \(p^\prime = 1(1p)^n\).
The ak Geometric Distribution Functions
The implementation of the geometric distribution follows our usual conventions, as demonstrated byak.geometricPMF
in listing 1.Listing 1: ak.geometricPMF
function pmf(p, k) { if(k>=0 && k===ak.floor(k)) return Math.pow(1p, k)*p; return isNaN(k) ? ak.NaN : 0; } ak.geometricPMF = function() { var state = {p: 0.5}; var arg0 = arguments[0]; var f, ll; constructors[ak.nativeType(arg0)](state, arg0, arguments); f = function(k){return pmf(state.p, k);}; f.p = function(){return state.p;}; return Object.freeze(f); }; var constructors = {};
Here the
pmf
helper function performs the actual calculation and the state
object is initialised by dispatch to a constructors
object. Listing 2 gives those for default and specified success probability initialisation.Listing 2: The Constructors
constructors[ak.UNDEFINED_T] = function() { }; constructors[ak.NUMBER_T] = function(state, p, args) { var arg1 = args[1]; constructors[ak.NUMBER_T][ak.nativeType(arg1)](state, arg1, args); state.p = Number(p); if(!(state.p>0 && state.p<1)) { throw new Error('invalid p in ak.geometric distribution'); } }; constructors[ak.NUMBER_T][ak.UNDEFINED_T] = function(state) { };
Note that, as we did for the Bernoulli distribution, we're constraining the probability of success to lie strictly between zero and one.
Program 1 plots the PMF for a range of success probabilities.
Program 1: The Geometric PMF



The CDF is even simpler, as shown by listing 3.
Listing 3: ak.geometricCDF
ak.geometricCDF = function() { var state = {p: 0.5}; var arg0 = arguments[0]; var f; constructors[ak.nativeType(arg0)](state, arg0, arguments); f = function(k){return k<0 ? 0 : 1Math.pow(1state.p, ak.floor(k)+1);}; f.p = function(){return state.p;}; return Object.freeze(f); };
Program 2 demonstrates its correctness by comparing the PMF at \(k\) to the difference between the CDF at \(k\) and at \(k1\).
Program 2: The PMF Versus The CDF



In order to cope with rounding errors, when calculating the inverse of the CDF we round down instead of up and then increment the result until the value of the CDF is no less than its argument, as done by the
invCDF
helper function in listing 4.Listing 4: ak.geometricInvCDF
function invCDF(p, ln1p, c) { var k; if(isNaN(c)) return ak.NaN; if(c>=1) return ak.INFINITY; if(c<=0) return 0; k = ak.floor(Math.log(1c)/ln1p  1); while(1Math.pow(1p, k+1)<c) ++k; return k; } ak.geometricInvCDF = function() { var state = {p: 0.5}; var arg0 = arguments[0]; var ln1p, f; constructors[ak.nativeType(arg0)](state, arg0, arguments); ln1p = Math.log(1state.p); f = function(c){return invCDF(state.p, ln1p, c);}; f.p = function(){return state.p;}; return Object.freeze(f); };
Note that
ak.geometricInvCDF
is precalculating the logarithm of \(1p\) to avoid recalculating it every time the helper function is called.Program 3 compares the inverse of the CDF's result to its argument.
Program 3: The Inverse CDF



To calculate the CF we shall exploit the facts that
\[
e^{i\theta} = \cos \theta + i \times \sin \theta
\]
and
\[
\frac{1}{x  i \times y}
= \frac{x + i \times y}{(x + i \times y) \times (x  i \times y)}
= \frac{x + i \times y}{x^2 + y^2}
\]
so that we can directly initialise an ak.complex
object with the real and imaginary parts of its results. The cf
helper function does so for ak.geometricCF
in listing 5Listing 5: ak.geometricCF
function cf(p, t) { var ct = Math.cos(t); var st = Math.sin(t); var re = 1(1p)*ct; var im = (1p)*st; var d = re*re+im*im; return ak.complex(p*re/d, p*im/d); } ak.geometricCF = function() { var state = {p: 0.5}; var arg0 = arguments[0]; var f; constructors[ak.nativeType(arg0)](state, arg0, arguments); f = function(t){return cf(state.p, t);}; f.p = function(){return state.p;}; return Object.freeze(f); };
and program 4 plots the real and imaginary parts of its results in black and red respectively.
Program 4: The Geometric CF



Finally we have the random number generator for which we need constructors to allow the specification of standard uniform random number generators other than the default
Math.random
, as provided in listing 6.Listing 6: ak.geometricRnd
ak.geometricRnd = function() { var state = {p: 0.5, rnd: Math.random}; var arg0 = arguments[0]; var ln1p, f; constructors[ak.nativeType(arg0)](state, arg0, arguments); ln1p = Math.log(1state.p); f = function(){return invCDF(state.p, ln1p, state.rnd());}; f.p = function(){return state.p;}; 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; };
Here we're simply calling the helper function
invCDF
with the result of the state.rnd
generator.Program 5 illustrates the use of
ak.geometricRnd
by comparing a histogram of its results to a graph of the geometric PMF.
Program 5: ak.geometricRnd



You will find the implementations of these functions in GeometricDistribution.js and we shall take a look at another property of Bernoulli processes next time.
References
[1] One Thing Or Another, www.thusspakeak.com, 2020.[2] How Long Must I Wait?, www.thusspakeak.com, 2015.
[3] The Longest Wait, www.thusspakeak.com, 2015.
Leave a comment