One Thing Or Another

| 0 Comments

Several years ago we took a look at memoryless processes in which the probability that we should wait for any given length of time for an event to occur is independent of how long we have already been waiting. We found that this implied that the waiting time must be exponentially distributed[1], that the waiting time for several events must be gamma distributed[2] and that the number of events occuring in a unit of time must be Poisson distributed[3].
These govern continuous memoryless processes in which events can occur at any time but not those in which events can only occur at specified times, such as the roll of a die coming up six, known as Bernoulli processes. Observations of such processes are known as Bernoulli trials and their successes and failures are governed by the Bernoulli distribution, which we shall take a look at in this post.

Heads Or Tails?

The Bernoulli distribution represents successful trials with a one and unsuccessful trials with a zero. The probability mass function, or PMF, for a trial with an outcome of \(X\) having a probability of success \(p\) is therefore
\[ f_p(x) = \Pr(X = x) = \begin{cases} 1-p & x = 0\\ p & x = 1\\ 0 & \text{otherwise} \end{cases} \]
and its cumulative distribution function, or CDF, is
\[ F_p(x) = \Pr(X \leqslant x) = \begin{cases} 0 & x<0\\ 1-p & 0 \leqslant x < 1\\ 1 & x \geqslant 1 \end{cases} \]
The inverse of the CDF is
\[ F_p^{-1}(c) = \begin{cases} 0 & c \leqslant 1-p\\ 1 & \text{otherwise} \end{cases} \]
and the characteristic function, or CF, is
\[ \hat{f_p}(t) = E\left[e^{itX}\right] = (1-p) \times e^{it \times 0} + p \times e^{it \times 1} = 1-p + p \times e^{it} \]
where \(E\) represents the expected value of the term between the following square brackets and \(i\) is the square root of minus one.

Finally, the outcomes of Bernoulli trials can be modelled using observations of the random variable of the standard uniform distribution \(U(0, 1)\), having equal chances of yielding any particular value between zero and one, with
\[ \begin{align*} U &\sim U(0, 1)\\ \\ X &= \begin{cases} 0 & U \leqslant 1-p\\ 1 & \text{otherwise} \end{cases} \end{align*} \]
where \(\sim\) means drawn from.

The ak Bernoulli Distribution Functions

Given how simple the functions associated with the Bernoulli distribution are, it should come as no surprise that their implementations are trivial, as illustrated by ak.bernoulliPMF in listing 1.

Listing 1: ak.bernoulliPMF
ak.bernoulliPMF = function() {
 var state = {p: 0.5};
 var arg0  = arguments[0];
 var f;

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

 f = function(k) {
  if(k===0) return 1-state.p;
  if(k===1) return state.p;
  return isNaN(k) ? ak.NaN : 0;
 };
 f.p = function(){return state.p;};
 return Object.freeze(f);
};

var constructors = {};

We're using our usual scheme of dispatching to a constructors object to initialise the state and listing 2 gives the implementation of the default and specified probability of success initialisers.

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.bernoulli distribution');
 }
};

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

One thing that might be a little surprising is that success probabilities of zero and one are treated as errors, even though all of the distribution functions are well defined for them. The reason for this is that for many of the other distributions associated with Bernoulli processes they aren't and for the sake of consistency it makes sense to treat them as invalid for all of them.
Note that this is similar to treating a standard deviation of zero as invalid for the normal distribution since in both cases their random variables are constants.

The implementation of the CDF is given by listing 3

Listing 3: ak.bernoulliCDF
ak.bernoulliCDF = function() {
 var state = {p: 0.5};
 var arg0  = arguments[0];
 var f;

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

 f = function(k) {
  if(k<0) return 0;
  if(k>=1) return 1;
  return isNaN(k) ? ak.NaN : 1-state.p;
 };
 f.p = function(){return state.p;};
 return Object.freeze(f);
};

and its inverse by listing 4.

Listing 4: ak.bernoulliInvCDF
ak.bernoulliInvCDF = function() {
 var state = {p: 0.5};
 var arg0  = arguments[0];
 var f;

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

 f = function(c) {
  if(c<=1-state.p) return 0;
  return isNaN(c) ? ak.NaN : 1;
 };
 f.p = function(){return state.p;};
 return Object.freeze(f);
};

To implement the CF we note that
\[ e^{i\theta} = \cos \theta + i \times \sin \theta \]
so that
\[ \hat{f_p}(t) = 1-p + p \times e^{it} = 1 - p + p \times \cos t + i \times p \times \sin t \]
as exploited by ak.bernoulliCF to directly initialise an ak.complex object, as shown by listing 5.

Listing 5: ak.bernoulliCF
ak.bernoulliCF = function() {
 var state = {p: 0.5};
 var arg0  = arguments[0];
 var f;

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

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

Program 1 plots the Bernoulli CF for a range of success probabilities with the real part in black and the imaginary part in red.

Program 1: The Bernoulli CF

Finally, we have the random number generator ak.bernoulliRnd which is implemented in listing 6.

Listing 6: ak.bernoulliRnd
ak.bernoulliRnd = function() {
 var state = {p: 0.5, rnd: Math.random};
 var arg0  = arguments[0];
 var f;

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

 f = function() {return state.rnd()<state.p ? 1 : 0;}
 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 adding a couple of constructors to allow the use of uniform random number generators other than Math.random. Furthermore, since they never return one we're using a strictly less than comparison against the success probability rather than a less than or equal to comparison.
Program 2 demonstrates its use by plotting a bar representing the number of failures in red beneath the number of successes in green.

Program 2: Using ak.bernoulliRnd

Whilst the Bernoulli distribution itself isn't particularly interesting, it's fundamental to the construction of Bernoulli processes which have some rather more complicated properties, as we shall see over the next few months. In the meantime you can find all of these functions in BernoulliDistribution.js.

References

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

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

[3] ...And Then Three Came Along At Once!, www.thusspakeak.com, 2015.

Leave a comment