Chapter 2: binomial and negative binomial plots

To generate the plots shown in Fig 2.4 we must first calculate binomial and negative binomial probabilities. Each of these requires us to compute factorials, which are defined as:

N! = factorial(N) = N x (N-1) x (N-2) x ... x 1.

These computations can easily lead to an overflow problem (see Appendix B), and we should make it a practice to compute factorials using the logarithmic gamma function. Without worrying about how the gamma function is computed, we can use the built-in matlab gamma function to compute factorials:

N! = factorial(N) = N x (N-1) x (N-2) x ... x 1 = gamma(N + 1),

where the log-factorial is therefore:

ln(N!) = gammaln(N + 1).

Further, we create a function called lnfactorial.m. This is done by creating a blank function under ‘New’, and filling the blank function with:

%lnfactorial.m

% USAGE: GF=lnfactorial(n)

% uses gamma functio to compute n!

%

% EXAMPLES:

% exp(lnfactorial(3)) %=6

% exp(lnfactorial(3.123)) %=7

function GF=lnfactorial(n)

GF=gammaln(n+1);

and then saving the function with the name lnfactorial

Furthermore, we see that the factorials appearing in (2.15) and (2.16) are used to compute the multiplicities, which enumerate the number of sequences of length N from which you can choose S successes (often abbreviated ‘N choose K’, where K stands in for the number of successes). Thus, we would also like to create a function lnchoosek.m that computes the log of the multiplicity:

% LNCHOOSEK.m

% USAGE: lnM=lnchoosek(N,K)

% computes lnM=ln[N!/K!(N-K)!], that is: ln[Binomial Coefficient], using the gammaln function

% [note that gammln supports fractional inputs N,K)

%

% OUTPUT: BC=ln[N!/K!(N-K)!]=ln(N!)-ln(K!)-ln([N-K]!)

% [the number of distinct configurations that have N successes in N observations]

% INPUTS: N=number of observations

% K=number of successes

function lnM=lnchoosek(n,k)

NUM=lnfactorial(n);

DEN=lnfactorial(k)+lnfactorial(n-k);

lnM=NUM-DEN; end

We can then use these functions to produce the binomial and negative binomial plots. We begin by computing the set of four binomial distributions shown in Fig. 2.4:

%%% Binomial distributions

Lbin=@(N,S,theta) lnchoosek(N,S)+S*log(theta)+(N-S)*log(1-theta);

Nnow=40;

thetas=[.1 .5 .7 .95];

Slist=cell(1,2); Lprob=cell(1,2);

Slist{1}=linspace(0,Nnow,Nnow+1)';

Slist{2}=linspace(0,Nnow,10*Nnow+1)';

Lprob{1}=Lbin(Nnow,Slist{1},thetas);

Lprob{2}=Lbin(Nnow,Slist{2},thetas);

We then plot the distributions on a single plot using ‘hold on’:

figure(1); clf; hold on

for n=1:4,

plot(Slist{2},exp(Lprob{2}(:,n)),'-','Color',.7*[1 1 1],'LineWidth',1); end

plot(Slist{1},exp(Lprob{1}(:,1)),'ko','MarkerFaceColor','k','MarkerSize',10)

plot(Slist{1},exp(Lprob{1}(:,2)),'kd','MarkerFaceColor',.65*[1 1 1],'MarkerSize',10)

plot(Slist{1},exp(Lprob{1}(:,3)),'ko','MarkerFaceColor',.35*[1 1 1],'MarkerSize',10)

plot(Slist{1},exp(Lprob{1}(:,4)),'ks','MarkerFaceColor','w','MarkerSize',10)

Similarly, we compute a set of negative binomial distributions, now holding the success rate constant at theta=0.5 and varying the number of total trials (N) on the abscissa for each of four target-numbers of successes, S:

%%% Negative binomial distributions

%%% selecting 4 S-values

Lnbin=@(N,S,theta) lnchoosek(N-1,S-1)+S*log(theta)+(N-S)*log(1-theta);

Nlist=cell(1,2); Lprob=cell(1,2);

Snow=[1 3 9 15];

thetas=0.5;

Nlist{1}=linspace(1,Nnow,Nnow)'; Nlist{2}=linspace(1,Nnow,10*Nnow)';

for n=1:length(Snow),

ntmp=Nlist{1}; ntmp(ntmp<Snow(n))=nan; Lprob{1}(:,n)=Lnbin(ntmp,Snow(n),thetas);

ntmp=Nlist{2}; ntmp(ntmp<Snow(n))=nan; Lprob{2}(:,n)=Lnbin(ntmp,Snow(n),thetas); end

Lprob{1}(isnan(Lprob{1}))=-inf; Lprob{2}(isnan(Lprob{2}))=-inf;

figure(2); clf; hold on

for n=1:4,

plot(Nlist{2},exp(Lprob{2}(:,n)),'-','Color',.7*[1 1 1],'LineWidth',1); end

plot(Nlist{1},exp(Lprob{1}(:,1)),'ko','MarkerFaceColor','k','MarkerSize',10)

plot(Nlist{1},exp(Lprob{1}(:,2)),'kd','MarkerFaceColor',.65*[1 1 1],'MarkerSize',10)

plot(Nlist{1},exp(Lprob{1}(:,3)),'ko','MarkerFaceColor',.35*[1 1 1],'MarkerSize',10)

plot(Nlist{1},exp(Lprob{1}(:,4)),'ks','MarkerFaceColor','w','MarkerSize',10)