[latex]\newcommand{\pr}[1]{P(#1)} \newcommand{\var}[1]{\mbox{var}(#1)} \newcommand{\mean}[1]{\mbox{E}(#1)} \newcommand{\sd}[1]{\mbox{sd}(#1)} \newcommand{\Binomial}[3]{#1 \sim \mbox{Binomial}(#2,#3)} \newcommand{\Student}[2]{#1 \sim \mbox{Student}(#2)} \newcommand{\Normal}[3]{#1 \sim \mbox{Normal}(#2,#3)} \newcommand{\Poisson}[2]{#1 \sim \mbox{Poisson}(#2)}[/latex]

Binomial Distribution

Counting occurrences of some random event is a common way of arriving at a discrete random variable. One of the simplest counting distributions is the Binomial distribution. As an example, let [latex]X[/latex] be the number of females in a sample of 3 individuals drawn independently from a population with a 1:1 sex ratio. There are 8 possible outcomes for the sexes of the 3 individuals:

FFF, FFM, FMF, FMM, MFF, MFM, MMF, MMM.

Since the probability of a female or male is 0.5, these 8 possibilities are all equally likely, with
\[ \pr{\mbox{FFF}} = 0.5 \times 0.5 \times 0.5 = 0.125, \]
and so on. If we count the number of females then we find that [latex]X[/latex] has only 4 values, not all equally likely. We can define the probability function on [latex]X[/latex] by [latex]\pr{X = 0} = 0.125[/latex], [latex]\pr{X=1} = 0.375[/latex], [latex]\pr{X=2} = 0.375[/latex], and [latex]\pr{X=3} = 0.125[/latex].
Because the random variable is discrete, we can list the values like this, or give them in a table, as in the table below.

Probability distribution of female count in a sample of 3

The probability function is plotted in the following figure.

A Bernoulli trial is a random process which has two possible outcomes, usually labelled “success” and “failure”. The count of successes in a series of independent Bernoulli trials, with a constant probability of success, has what is called the Binomial distribution. If [latex]X[/latex] is the count from [latex]n[/latex] trials with success probability [latex]p[/latex] then we write [latex]\Binomial{X}{n}{p}[/latex].

The probability distribution in the previous table is Binomial(3,0.5) distribution. Let us suppose that the female probability was [latex]p[/latex] and see where the probabilities for Binomial(3,[latex]p[/latex]) come from. Assuming the trials are independent, the probability of getting three females, FFF, is
\[ \pr{X = 3} = p \times p \times p = p^3. \]
The probability of FFM, 2 females, is [latex]pp(1-p) = p^2(1-p)[/latex]. But we can also get 2 females through FMF or MFF, each also with probability [latex]p^2(1-p)[/latex], so
\[ \pr{X = 2} = 3p^2(1-p). \]
Similarly, there are 3 combinations which give 1 female, each with probability [latex]p(1-p)^2[/latex], so
\[ \pr{X = 1} = 3p(1-p)^2. \]
Finally,
\[ \pr{X = 0} = (1-p)^3. \]
This distribution is summarised in the following table. You can check that when [latex]p = 0.5[/latex] we get the values shown in our earlier table.

Probability function for the number of females in a sample of size 3

The general formula for [latex]\pr{X = x}[/latex] when [latex]\Binomial{X}{n}{p}[/latex] is
\[ \pr{X = x} = {n \choose x} p^x(1-p)^{n-x}, \]
where
\[ {n \choose x} = \;^nC_x = \frac{n(n-1)\cdots(n-x+1)}{x(x-1) \cdots 1} = \frac{n!}{x!(n-x)!} \]
is the Binomial coefficient, available on many calculators. The Binomial coefficient counts the number of ways in which [latex]x[/latex] objects can be chosen from [latex]n[/latex] possibilities when order doesn’t matter. For example, suppose we take a sample of 2 people from a group of 10. There are [latex]{10 \choose 2} = 45[/latex] different pairs of people we could choose.

The binomial distribution table gives Binomial probabilities for some simple values of [latex]p[/latex] and small values of [latex]n[/latex] for use in our examples. The cumulative binomial distribution table gives cumulative Binomial probabilities, which will also be useful (although they can be calculated by adding up the corresponding values in the original table. In practice software will calculate both of these for us.

Note that if [latex]p = 0.5[/latex] then the Binomial distribution is symmetric (see the previous figure and this figure in Chapter 17). If [latex]p \lt 0.5[/latex] then the distribution is skewed to the right (see this other figure in Chapter 17) while if [latex]p \gt 0.5[/latex] then it is skewed to the left.

Finite Populations

An obvious use of Binomial distributions is in modelling survey data. Suppose 50% of all university students are male and that we take a random sample of 100 students. The probability of an individual student being male is 0.5 and so our survey seems to be a series of Bernoulli trials and so the count of “yes” would be Binomial.

However, if the first person we survey is male then the probability in the remainder of the population will be slightly less than 0.5. Thus the trials will not be strictly independent. This effect will be small, especially if the population is much larger than the sample size. We will thus use the Binomial distribution to model such surveys.

If the population is finite and the size of the population is known then you can also use the hypergeometric distribution (in a later section) to calculate exact sample probabilities.

Sample Counts and Proportions

The rules in Chapter 10 can be used to derive formulas for the expected value and standard deviation of more complicated random variables, such as one with a Binomial distribution. This section shows the algebra needed for obtaining the formulas. You should read through this and try to follow the steps. There is no need to learn the proofs but it is useful to see that all of the formulas which we will encounter later in the book are not at all magical. They can all be shown using such fairly simple (though sometimes messy) algebra.

Sample Count

So let [latex]\Binomial{X}{n}{p}[/latex], the count of successes from [latex]n[/latex] independent Bernoulli trials each with success probability [latex]p[/latex]. To find the expected value and standard deviation of [latex]X[/latex], define new random variables [latex]B_1, B_2, \ldots, B_n[/latex] by
\[ B_j = \left\{ \begin{array}{ll}
1, & \mbox{ if }j\mbox{ th trial is success} \\
0, & \mbox{ if }j\mbox{ th trial is failure}. \\
\end{array} \right.
\]
Since each [latex]B_j[/latex] has only two possible outcomes it is easy to give its probability function, shown in the table below.

Probability function for [asciimath]B_j[/asciimath]

We can work out the expected value and standard deviation directly from the definitions. Firstly,
\[ \mean{B_j} = p \times 1 + (1-p) \times 0 = p. \]
For example, if we pick a person at random from a population with 60% females then this says on average we will get 0.6 of a female.

Note that [latex]p[/latex] has two roles here and this is what sometimes makes working with sample counts and proportions a bit confusing. In the previous table, [latex]p[/latex] was a probability, the chance of a success. In [latex]\mean{B_j} = p[/latex] it now has units because it tells us what to expect for [latex]B_j[/latex] and [latex]B_j[/latex] is, for example, the number of females from a single Bernoulli trial. We have a 0.6 chance of getting a female which results in an average of 0.6 females per trial. Always the same number but a different interpretation.

The variance of each [latex]B_j[/latex] is
\begin{eqnarray*}
\var{B_j} & = & p \times (1-p)^2 + (1-p) \times (0 – p)^2 \\
& = & p(1 – 2p + p^2) + (1-p)p^2 \\
& = & p – 2p^2 + p^3 + p^2 – p^3 \\
& = & p – p^2 \\
& = & p(1 – p)
\end{eqnarray*}
so
\[ \sd{B_j} = \sqrt{p(1 – p)}. \]
Again this is in units of whatever we are counting.

Now what is [latex]B_1 + B_2 + \cdots + B_n[/latex]? Each time there is a success the corresponding [latex]B_j[/latex] is 1. So, for example, if we get 7 successes out of 10 trials then we will be adding up 7 lots of 1. This, of course, is 7. This sum will always give the count of successes. That is,
\[ X = B_1 + B_2 + \cdots + B_n. \]
Since we know the expected value and standard deviation of the [latex]B_j[/latex], we can use this to find the expected value and standard deviation of [latex]X[/latex]. Firstly,
\begin{eqnarray*}
\mean{X} & = & \mean{B_1 + B_2 + \cdots + B_n} \\
& = & \mean{B_1} + \mean{B_2} + \cdots + \mean{B_n} \\
& = & p + p + \cdots + p \\
& = & np.
\end{eqnarray*}
For standard deviation, since the [latex]B_j[/latex] are independent we have
\begin{eqnarray*}
\var{X} & = & \var{B_1 + B_2 + \cdots + B_n} \\
& = & \var{B_1} + \var{B_2} + \cdots + \var{B_n} \\
& = & p(1-p) + p(1-p) + \cdots + p(1-p) \\
& = & np(1-p).
\end{eqnarray*}
so [latex]\sd{X} = \sqrt{np(1-p)}[/latex].

For example, if 60% of a population are female and we pick a random sample of 10 individuals then we would expect to get [latex]10 \times 0.6 = 6[/latex] females in our sample. The standard deviation of the sample count is
\[ \sd{X} = \sqrt{10 \times 0.6 \times 0.4} = 1.55. \]
Note that these numbers have units — we expect to find 6 females in the sample with a standard deviation of 1.55 females.

Sample Proportion

The sample proportion is
\[ \hat{P} = \frac{X}{n} \]
so we can use the formulas from Chapter 10 to find
\[ \mean{\hat{P}} = \mbox{E}\left(\frac{X}{n}\right) = \frac{1}{n} \mean{X} = \frac{1}{n} (np) = p. \]
This is not surprising. If you take a sample of any size from a population with 60% females then you would expect your sample to have 60% females. We thus call [latex]\hat{P}[/latex] an unbiased estimator of [latex]p[/latex]. The standard deviation is
\[ \sd{\hat{P}} = \mbox{sd}\left(\frac{X}{n}\right) = \frac{1}{n} \sqrt{np(1-p)} = \sqrt{\frac{p(1-p)}{n}}. \]
Note that [latex]\sd{\hat{P}}[/latex] is proportional to [latex]1/\sqrt{n}[/latex], just as for [latex]\sd{\overline{X}}[/latex].

Maximum Likelihood Estimation

The sample proportion seems a sensible estimate of the population proportion since we’ve seen that [latex]\mean{\hat{P}} = p[/latex], that on average the sample proportion gives the population proportion. This is the main justification for estimators that we will use in this book. However, there is another important method that you should be familiar with, known as maximum likelihood estimation (MLE).

To see how this works we start by asking a slightly different question to before. For example, suppose we pick a sample of [latex]n= 20[/latex] people and find that [latex]X = 16[/latex] drink coffee. How can we estimate the proportion, [latex]p[/latex], of all the people that our sample came from that drink coffee?

The maximum likelihood approach is to consider all possible values for [latex]p[/latex] and work out the probability of actually getting [latex]X = 16[/latex] people who drink coffee in our sample. From the Binomial formula this is
\[ {20 \choose 16} p^{16}(1-p)^{20-16} = 4845p^{16}(1-p)^4. \]
For example, if [latex]p = 0.6[/latex] then [latex]\pr{X=16} = 0.035[/latex], while if [latex]p = 0.75[/latex] then [latex]\pr{X=16} = 0.130[/latex]. Thus it is more likely to get 16 coffee drinkers if [latex]p = 0.75[/latex] than if [latex]p = 0.6[/latex]. The obvious question then is which value of [latex]p[/latex] makes it most likely to get 16 coffee drinkers in our sample?

You can solve this question using calculus but here we will just look at the figure below. This is a plot of the likelihood function, [latex]\pr{X=16 \; | \; p}[/latex], the probability of obtaining our data for different values of the unknown population proportion [latex]p[/latex]. We see the value that maximises the likelihood function is [latex]p = 0.8[/latex], so 0.8 is our maximum likelihood estimate of the proportion of coffee drinkers in our population.

This isn’t a terribly surprising example since the ML estimate is just [latex]\hat{p}[/latex], the sample proportion! However, MLE can be used in more complicated settings where there is no obvious estimate for the population parameter of interest.

Hypergeometric Distribution

We have noted that the Binomial distribution is not strictly correct as a model for sampling from a finite population. To be precise, once we have made our measurement of success or failure on the first individual, the distribution of the remainder will change slightly, resulting in observations that are not independent. If the population is sufficiently large, at least 10 times the sample size, then this effect will be negligible and so we have worked with this approximation quite happily.

However, it is possible to calculate exact probabilities for counts from finite populations. Suppose [latex]X[/latex] is the count of successes in a sample of size [latex]n[/latex] from a population of size [latex]N[/latex]. For the Binomial distribution we also needed the success probability [latex]p[/latex] to make calculations. However, here we can describe the whole distribution by saying that [latex]N = N_1 + N_2[/latex], where [latex]N_1[/latex] is the number of “successes” and [latex]N_2[/latex] is the number of failures. For the first sample we would have [latex]p = \frac{N_1}{N}[/latex] but this would change immediately after the first sample.

As a specific example, suppose we take a sample of [latex]n[/latex] = 10 people from a group with [latex]N_1[/latex] = 20 females and [latex]N_2[/latex] = 30 males. What is the probability of getting exactly [latex]X = 6[/latex] females in our sample?

Firstly, the number of different samples we could possibly get is the number of ways of choosing 10 things from a collection of [latex]20 + 30 = 50[/latex] things. This is [latex]{50 \choose 10} = 10272278170[/latex]. Each of these is equally likely and we want to know how many have 6 females.

Now there are 6 spots in our sample for females and there are 20 females to choose from. Thus the number of ways we can fill these 6 spots is [latex]{20 \choose 6} = 38760[/latex]. For each of these possibilities we need to see how many ways we can fill the remaining 4 spots with males. This is similarly [latex]{30 \choose 4} = 27405[/latex]. Together the total number of samples that would give 6 females is
\[ 38760 \times 27405 = 1062217800. \]
Thus the probability of 6 females in the sample of 10 is
\[ \pr{X = 6} = \frac{1062217800}{10272278170} = 0.1034. \]
A Binomial calculation with [latex]p = \frac{20}{50} = 0.4[/latex] fixed gives [latex]\pr{X = 6} = 0.1115[/latex], showing the error involved in the Binomial assumption of independence.

Following these steps in general gives the probability of getting [latex]x[/latex] successes out of [latex]n[/latex] trials as
\[ \pr{X = x} = \frac{{N_1 \choose x} {N_2 \choose n-x}}{{N \choose n}}. \]
If [latex]X[/latex] has these probabilities then we say that it follows a hypergeometric distribution.

Poisson Distribution

William Gosset, a chemist at the Guiness Brewery in Dublin, described an experiment in which he spread a thin layer of mixture containing yeast cells, gelatine and water on a plate of glass. A small area was divided into 400 squares of equal size and a hemocytometer was used to count the number of yeast cells in each square. Guinness did not allow employees to publish at the time, so his work was published using the pseudonym Student (1907). The experimental results are reproduced here in the Appendix.

The distribution of the 400 counts is shown in the figure above. How should we model this data? It looks like a Binomial distribution because it is a count of outcomes. But what are the values of [latex]n[/latex] and [latex]p[/latex]? We can imagine a very fine grid with [latex]n[/latex] holes in it, each the size of a yeast cell, and each hole has a small probability [latex]p[/latex] of having a yeast cell in it. Yet [latex]n[/latex] is rather indefinite, some very large number that we cannot really quantify, and [latex]p[/latex] is correspondingly tiny.

However, we can still talk about the average yeast cell count, [latex]np[/latex]. From the data in the Appendix we find a sample mean of [latex]\overline{x} = 4.68[/latex] yeast cells per unit of area. We will model this type of data by working directly with the mean.

Gosset derived a “law” related to the Binomial distribution for this case of events where [latex]p[/latex] is very small and [latex]n[/latex] is very large. Since [latex]p[/latex] is small it is often referred to as the “law of rare events”. It gives the probability of a count [latex]x[/latex] by
\[ \pr{X = x} = \frac{e^{-\lambda} \lambda^x}{x!}, \]
where [latex]\lambda[/latex] is the mean of the counts. The table below shows a comparison between the observed values and the expected values from this distribution when [latex]\lambda = 4.68[/latex]. There seems to be a reasonable match, which we will return to in Chapter 22.

Observed and expected counts of yeast cells

It turns out that Gosset had rediscovered a distribution that had earlier been described by the French mathematician Siméon Poisson in his 1837 book Recherchés sur la probabilité des jugements en matiére criminelle et matiére civile. This was a significant work but was not widely known in that century. In 1898 Ladislaus von Bortkiewicz published Das Gesetz der kleinen Zahlen (The Law of Small Numbers) which gave examples of a range of cheery settings in which the Poisson distribution arose, such as numbers of suicides and the numbers of deaths from horse kicks in the Prussian Army. Bulmer (1979) and Feller (1970) give further examples and applications.

As it was Bortkiewicz who first showed the wide applicability of the distribution, some authors have suggested it be called the “Bortkiewicz distribution”. However, if a random variable [latex]X[/latex] has the probability function given above then we say it has a Poisson distribution and write [latex]\Poisson{X}{\lambda}[/latex].

Aggregation and Over-dispersion

Since [latex]X[/latex] is discrete we can use the formulas from Chapter 10 to calculate the expected value and standard deviation of [latex]X[/latex]. However, unlike the Binomial distribution, here there is no upper limit for [latex]X[/latex] from a Poisson distribution. Even though it becomes very unlikely that large values could occur, [latex]\pr{X = x}[/latex] never equals 0. This just means that our sums are infinite sums. Such sums don’t always converge to a number but here they do and we find
\[ \mean{X} = \sum_{x=0}^{\infty} \frac{e^{-\lambda} \lambda^x}{x!} x = \lambda, \]
and then
\[ \var{X} = \sum_{x=0}^{\infty} \frac{e^{-\lambda} \lambda^x}{x!} (x – \lambda)^2 = \lambda. \]
We have implicitly used [latex]\mean{X} = \lambda[/latex] in the above examples. But the variance is also [latex]\lambda[/latex], and so looking at the ratio
\[ \mbox{CD } = \frac{s^2}{\overline{x}}, \]
the coefficient of dispersion, gives a simple way of seeing whether sample data comes from a Poisson distribution. This value, alternatively known as the index of dispersion or the variance-to-mean ratio, should be around 1 for Poisson data.

Often in biological studies there will be a larger number of extreme values than the Poisson distribution predicts. For example, if we were counting the spatial distribution of some vegetation then we would often find cells where a lot of vegetation is clumped, giving high counts. This aggregation gives a higher standard deviation, so we would find [latex]s^2 \gt \overline{x}[/latex], compared to the expected [latex]s^2 = \overline{x}[/latex] for a Poisson distribution. This would give a CD value greater than 1.

It is also possible to have a distribution where [latex]\sigma^2 \lt \mu[/latex]. Berkson et al. (1935) carried out a similar hemocytometer study to the one presented by Gosset but counting red blood cells instead of yeast cells. Over a number of experiments they found the general relation [latex]\sigma^2 = 0.85\mu[/latex]. This result is thought attributable to the forces acting between the red blood cells which encourages them to form a more uniform distribution across a surface, the opposite of the aggregation found in other biological settings. With a more uniform distribution, known as over-dispersion, the variability in the counts will be lower, leading to a lower [latex]s^2[/latex]. The result will be a value of CD less than 1.

Dispersion?

We have used terminology from ecology here where we take ‘dispersed’ to mean spread out in space. However, the term ‘dispersion’ is also used as a synonym for ‘variability’ in many other areas of statistics. In that case ‘over-dispersed’ would mean the variability is higher than it should be, corresponding to what we have called ‘aggregated’ above. Conversely, ‘under-dispersed’ would mean that the variability is less than it should be, precisely what we have called ‘over-dispersed’! This is something to be aware of but generally it will be clear from the context and the values of the statistics which terminology is being used.

Exercise 9

The figure below gives an example of spatial data, similar to the yeast cell data. As Gosset did, divide the area into a grid, such as an 8-by-8 grid of 64 squares. Count the number of individuals in each square and make an appropriate plot of the count data. Calculate the mean count and use it to estimate [latex]\lambda[/latex] in a Poisson distribution. Using the Poisson formula, calculate the estimated proportions of each count and compare these to your observations. Does a Poisson model of this data seem appropriate?

Binomial distribution

This table gives [asciimath]P(X=x)[/asciimath] for [latex]X\sim\mbox{Binomial}(n,p)[/latex]. For [asciimath]p \ge 0.50[/asciimath] reverse your question by counting failures instead of successes.

Cumulative Binomial distribution

This table gives [latex]P(X \ge x)[/latex] for [latex]X\sim\mbox{Binomial}(n,p)[/latex]. These probabilities will give quicker and more accurate results than can be obtained by summing the individual probabilities from the previous table. For [latex]p \ge 0.50[/latex] reverse your question by counting failures instead of successes.