If At First You Don't Succeed

One Thing Or Another

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.

A Series Of Unfortunate Events

Since the trials are independent, the probability of having \(k\) failures in a row is trivially

\[ (1-p)^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) = (1-p)^k \times p \]

for non-negative 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} (1-p)^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-(1-p)^{\lfloor k \rfloor +1}\right)}{1 - (1-p)} = 1 - (1-p)^{\lfloor k \rfloor + 1} & k \geqslant 0\\ 0 & \text{otherwise} \end{cases}\\ \]

Derivation 1: The Inverse CDF

Firstly, define

\[ \begin{align*} c &= 1 - (1-p)^{k+1}\\ 1-c &= (1-p)^{k+1} \end{align*} \]

Taking logarithms yields

\[ \begin{align*} \ln (1-c) &= (k+1) \times \ln (1-p)\\ k &= \frac{\ln (1-c)}{\ln (1-p)} - 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\).

because of which the distribution of \(K\) is called the geometric distribution. The inverse of the CDF is given by

\[ F_p^{-1}(c) = \Bigg\lceil \frac{\ln (1-c)}{\ln (1-p)} - 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 (1-p)^k \times p \times e^{itk}\\ &= \sum_{k=0}^\infty p \times \left((1-p)\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((1-p) \times e^{it}\right)^\infty}{1-(1-p) \times e^{it}}\\ &= \frac{p}{1-(1-p) \times e^{it}} \end{align*} \]

provided that \(1-p\) 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) = (1-p)^{\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)\\ &= (1-p)^{\lfloor k \rfloor + 1} \times (1-p)^{\lfloor k \rfloor + 1} \times \dots \times (1-p)^{\lfloor k \rfloor + 1}\\ &= \left((1-p)^{\lfloor k \rfloor + 1}\right)^n = \left((1-p)^{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(1-p^\prime\right)^{\,\lfloor k \rfloor + 1} = F_{p^\prime}(k) \]

where \(p^\prime = 1-(1-p)^n\).

The ak Geometric Distribution Functions

The implementation of the geometric distribution follows our usual conventions, as demonstrated by ak.geometricPMF in listing 1.

Listing 1: ak.geometricPMF

function pmf(p, k) {
 if(k>=0 && k===ak.floor(k)) return Math.pow(1-p, 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
var n = 100; var p0 = 0.9; var p1 = 0.1; var dp = (p1-p0)/n; var h = new Array(20); var ub = h.length; var x_ticks = [0,ub/2,ub]; var y_ticks = [0,p0/2,p0]; var k; bar.bounds([[0,0],[ub,p0]]); ak.asynchLoop ( function body() { var pmf = ak.geometricPMF(p0); for(k=0;k<h.length;++k) h[k] = pmf(k); bar.clear(); bar.plot(h, {fillStyle:'green'}); bar.axes([0,0], [x_ticks, y_ticks]); p0 += dp; return --n>0; }, undefined, undefined, 50 );

Program 1: The Geometric PMF

var n = 100;
var p0 = 0.9;
var p1 = 0.1;
var dp = (p1-p0)/n;
var h = new Array(20);
var ub = h.length;
var x_ticks = [0,ub/2,ub];
var y_ticks = [0,p0/2,p0];
var k;

bar.bounds([[0,0],[ub,p0]]);

ak.asynchLoop (
 function body() {
  var pmf = ak.geometricPMF(p0);
  for(k=0;k<h.length;++k) h[k] = pmf(k);
  bar.clear();
  bar.plot(h, {fillStyle:'green'});
  bar.axes([0,0], [x_ticks, y_ticks]);

  p0 += dp;
  return --n>0;
 },
 undefined,
 undefined,
 50
);

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 : 1-Math.pow(1-state.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 \(k-1\).

Program 2: The PMF Versus The CDF
var p = 0.15; var f = ak.geometricPMF(p); var F = ak.geometricCDF(p); var n = 16; var k, fk, Fk0, Fk1, dF, d; console.writeln('k\tf(k)\tdF(k)\tf(k)-dF(k)'); for(k=0;k<n;++k) { fk = f(k); Fk0 = F(k); Fk1 = F(k-1); dF = (Fk0-Fk1); d = fk-dF; console.write(k, '\t'); console.write(fk.toFixed(4), '\t'); console.write(dF.toFixed(4), '\t'); console.write(d.toExponential(4), '\n'); }	k f(k) dF(k) f(k)-dF(k) 0 0.1500 0.1500 -2.7756e-17 1 0.1275 0.1275 -5.5511e-17 2 0.1084 0.1084 -1.3878e-17 3 0.0921 0.0921 -4.1633e-17 4 0.0783 0.0783 -1.3878e-17 5 0.0666 0.0666 -4.1633e-17 6 0.0566 0.0566 -2.7756e-17 7 0.0481 0.0481 -4.1633e-17 8 0.0409 0.0409 7.6328e-17 9 0.0347 0.0347 -4.1633e-17 10 0.0295 0.0295 -6.9389e-18 11 0.0251 0.0251 -6.9389e-18 12 0.0213 0.0213 -3.4694e-17 13 0.0181 0.0181 2.0817e-17 14 0.0154 0.0154 -1.7347e-17 15 0.0131 0.0131 2.7756e-17

Program 2: The PMF Versus The CDF

var p = 0.15;

var f = ak.geometricPMF(p);
var F = ak.geometricCDF(p);

var n = 16;
var k, fk, Fk0, Fk1, dF, d;

console.writeln('k\tf(k)\tdF(k)\tf(k)-dF(k)');

for(k=0;k<n;++k) {
 fk = f(k);
 Fk0 = F(k);
 Fk1 = F(k-1);
 dF = (Fk0-Fk1);
 d = fk-dF;
 console.write(k, '\t');
 console.write(fk.toFixed(4), '\t');
 console.write(dF.toFixed(4), '\t');
 console.write(d.toExponential(4), '\n');
}

k	f(k)	dF(k)	f(k)-dF(k)
0	0.1500	0.1500	-2.7756e-17
1	0.1275	0.1275	-5.5511e-17
2	0.1084	0.1084	-1.3878e-17
3	0.0921	0.0921	-4.1633e-17
4	0.0783	0.0783	-1.3878e-17
5	0.0666	0.0666	-4.1633e-17
6	0.0566	0.0566	-2.7756e-17
7	0.0481	0.0481	-4.1633e-17
8	0.0409	0.0409	7.6328e-17
9	0.0347	0.0347	-4.1633e-17
10	0.0295	0.0295	-6.9389e-18
11	0.0251	0.0251	-6.9389e-18
12	0.0213	0.0213	-3.4694e-17
13	0.0181	0.0181	2.0817e-17
14	0.0154	0.0154	-1.7347e-17
15	0.0131	0.0131	2.7756e-17

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(1-c)/ln1p - 1);
 while(1-Math.pow(1-p, 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(1-state.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 pre-calculating the logarithm of \(1-p\) 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
var p = 0.15; var F = ak.geometricCDF(p); var I = ak.geometricInvCDF(p); var n = 16; var k, Fk, IFk; console.writeln('k\tF(k)\tInv(F(k))'); for(k=0;k<n;++k) { Fk = F(k); IFk = I(Fk); console.write(k, '\t'); console.write(Fk.toFixed(4), '\t'); console.write(IFk, '\n'); }	k F(k) Inv(F(k)) 0 0.1500 0 1 0.2775 1 2 0.3859 2 3 0.4780 3 4 0.5563 4 5 0.6229 5 6 0.6794 6 7 0.7275 7 8 0.7684 8 9 0.8031 9 10 0.8327 10 11 0.8578 11 12 0.8791 12 13 0.8972 13 14 0.9126 14 15 0.9257 15

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 5

Listing 5: ak.geometricCF

function cf(p, t) {
 var ct = Math.cos(t);
 var st = Math.sin(t);
 var re = 1-(1-p)*ct;
 var im = (1-p)*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
var p = 0.25; var tb = 20; var cf = ak.geometricCF(p); var n = 100; var re = new Array(n2+1); var im = new Array(n2+1); var i, t, z; graph.bounds([[-tb,-1],[tb,1]]); graph.axes([0,0], [[-tb,-tb/2,tb/2/tb], [-1,-0.5,0.5,1]]); for(i=-n;i<=n;++i) { t = tb * i/n; z = cf(t); re[i+n] = [t, z.re()]; im[i+n] = [t, z.im()]; } graph.plot(re); graph.plot(im, {strokeStyle:'red'});

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(1-state.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
var m = 10; var draws = 10; var times = 1000; var p = 0.5; var pmf = ak.geometricPMF(p); var rnd = ak.geometricRnd(p); var sum = new Array(m); var geo = new Array(m); var max = 0; var c = 0; var i; for(i=0;i<m;++i) sum[i] = 0; ak.asynchLoop ( function body() { var k; for(i=0;i<draws;++i) { k = rnd(); if(k<m) { sum[k] += 1; max = Math.max(sum[k], max); } ++c; } for(k=0;k<m;++k) { geo[k] = [k+0.5, pmf(k)c]; } bar.bounds([[0,0],[m,max1.1]]); graph.bounds([[0,0],[m,max*1.1]]); bar.plot(sum, {fillStyle:'green'}); graph.plot(geo); return --times>0; } );

Program 5: ak.geometricRnd

var m = 10;
var draws = 10;
var times = 1000;

var p = 0.5;
var pmf = ak.geometricPMF(p);
var rnd = ak.geometricRnd(p);

var sum = new Array(m);
var geo = new Array(m);

var max = 0;
var c = 0;
var i;

for(i=0;i<m;++i) sum[i] = 0;

ak.asynchLoop (
 function body() {
  var k;

  for(i=0;i<draws;++i) {
   k = rnd();
   if(k<m) {
    sum[k] += 1;
    max = Math.max(sum[k], max);
   }
   ++c;
  }

  for(k=0;k<m;++k) {
   geo[k] = [k+0.5, pmf(k)*c];
  }

  bar.bounds([[0,0],[m,max*1.1]]);
  graph.bounds([[0,0],[m,max*1.1]]);
  bar.plot(sum, {fillStyle:'green'});
  graph.plot(geo);
  return --times>0;
 }
);

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.


AKCalc	ECMA


Endarkenment	Turning Sixteen

Search

If At First You Don't Succeed

Previous:

A Series Of Unfortunate Events

The ak Geometric Distribution Functions

References

Previous:

Categories:

Tags:

Leave a comment

Categories

Gallimaufry

Links

	Chrome 26+
	Firefox 20+
	Internet Explorer 9+

Sign In

Search

If At First You Don't Succeed

Previous:

A Series Of Unfortunate Events

The ak Geometric Distribution Functions

References

Previous:

Categories:

Tags:

Leave a comment

Tag Cloud

Categories

Gallimaufry

Links