Lab 12: Plotting Probability Distributions

In this lab, you will use Mathematica to calculate the expected
values of various distributions, to draw random numbers from 
a distribution, and to plot probability distributions.

We will focus on the binomial distribution and two distributions
that can be used to approximate the binomial distribution
(the Poisson and the normal distributions).

The distributions

Mathematica knows several distributions.  First we have to 
load special packages that contain these distributions.
Shift-Enter the following:

<<Statistics`DiscreteDistributions`
<<Statistics`ContinuousDistributions`
<<Graphics`Graphics`
Off[General::"spell1"]
Off[General::"spell"]

Once these distributions are loaded into memory,
Mathematica represents an entire distribution in a command,
such as BinomialDistribution[n, p]:

?BinomialDistribution

(The Binomial distribution requires that you give it information
about the number of trials, n, and the probability of a success
in each  trials, p.)

To find the probability of observing a number of  successes,
say k, we use the command PDF (which stands
for probability density function):

PDF[BinomialDistribution[n,p],k]

The if statements ensure that k is a positive integer and
Binomial[n,k] is "n choose k" or just n!/(k!(n-k)!).

Here is the information that Mathematica has on the
Poisson Distribution:

?PoissonDistribution

(The Poisson distribution only requires information about the mean, mu.)

PDF[PoissonDistribution[mu],k]

Here is the information that Mathematica has on the
Normal Distribution:

?NormalDistribution

(The Normal distribution requires that you give it information
about the mean, mu, and standard deviation, sigma, of the distribution.
sigma = Sqrt[variance])

PDF[NormalDistribution[mu,Sqrt[var]],k]

Means and Variances

Mathematica knows the mean, variance, and 
skewness (= a measure of asymmetry) for all the common distributions.

For the Binomial:

Mean[BinomialDistribution[n, p]]

Variance[BinomialDistribution[n, p]]

Skewness[BinomialDistribution[n, p]]

For what value of p is the binomial symmetrix (zero skew)?

For the Poisson:

Mean[PoissonDistribution[mu]]

Variance[PoissonDistribution[mu]]

Skewness[PoissonDistribution[mu]]

For the Normal Distribution:

Mean[NormalDistribution[mu,Sqrt[var]]]

Variance[NormalDistribution[mu,Sqrt[var]]]

Skewness[NormalDistribution[mu,Sqrt[var]]]

Take-home points

Look at the above means, variances, and skewnesses
to verify the following points.

Setting mu=n*p in the Poisson distribution,
the Poisson distribution can be used to approximate the
binomial distribution when p is small.
(In this case, the Poisson and the binomial have nearly
identical mean, variance, and skew.)

Setting mu=n*p, and var=n*p*(1-p) in the normal distribution,
the normal distribution can be used to approximate the
binomial distribution when n is large and/or when p is near 1/2.
(In this case, the normal and the binomial have nearly
identical mean, variance, and skew.)

We'll look at these points graphically in the last sections.

Biological Model

In this section, we use Mathematica to draw random numbers from a binomial
distribution to simulate the fate of a rare beneficial mutation.

SCENARIO:  In a population of n haploid adults, a new mutation appears in a single 
copy (p=1/n).  It improves fitness by an amount s, which we'll assume means that 
individuals carrying the mutation produce more offspring.

Using the haploid model of selection from class, the frequency of the new
mutant among offspring would be:

newfreq[p_, s_] = 
      p (1 + s)/(p (1 + s)  +  (1 - p))

Although many offspring are produced, only n survive to adulthood.  We can
thus draw a random number from the Binomial Distribution to determine how many
copies of the allele are found among the surviving offspring.  Here we'll assume that the allele
increases fecundity by 10% in a population of size n=1000:

s = 0.1;
n = 1000;

Random[BinomialDistribution[n, newfreq[1/n, s]]]

If you enter the above command over and over again, you'll see that the allele is often
lost, even though it confers a fitness advantage.

We can use a "While" loop to determine when the allele is lost, if it is lost, and 
when the allele is fixed, if it is fixed.

?While

(Make sure that the following coding makes sense. MutantNumber stands for the number
of mutants among the adults and is assumed to equal one right after the mutation appears. )

MutantNumber = 1;
time = 0;
While[
      MutantNumber > 0  &&  
        MutantNumber < n, 
MutantNumber = 
        Random[BinomialDistribution[n, 
            newfreq[MutantNumber/n, s]]]; 
time = 
        time + 1
]
If[MutantNumber == 0,  
      Print["Allele lost at time ", time]]
If[MutantNumber == n,  
      Print["Allele fixed at time ", time]]

We can use a "For" loop to repeat this process, e.g., 100 times, to assess how often the 
mutation is lost and how often the mutation is fixed:

numlost = 0;
numfixed = 0;
For[trial = 1, 
      trial < 1001, trial++, 
MutantNumber = 
        1; 
time = 
        0; 
While[
          MutantNumber > 0  &&  
            MutantNumber < n, 
MutantNumber = 
            Random[BinomialDistribution[n, 
                newfreq[MutantNumber/n, s]]]; 
time = 
            time + 1
]
If[
          MutantNumber == 0, 
          numlost = numlost + 1]
If[
          MutantNumber == n,  
          numfixed = 
            numfixed + 1]
]
numlost
numfixed

CLASS:  Write the % of times that you observed fixation on the board.

How do these simulations compare to Haldane's solution for the probability of fixation?  

Haldane used the Poisson distribution to approximate the Binomial, which is 
a reasonable assumption given that the mutation is initially rare.

Haldane found that the fixation probability, P, solves:

s = 0.1;

Solve[1 - P  ==  Exp[-(1 + s) P], P]

Although there are two numerical solutions, one is near zero and is not the one of
interest.  The second is gives the probability of fixation (P=0.176 for s=0.1).
This is a more accurate solution than 2s (P=0.2 for s=0.1), which assumes that s is small.

Graphs

In this section, we'll use bar charts to compare the Poisson and Normal distributions
to the Binomial distribution.

Before you begin, in which of the following cases do you expect 

(a) the Poisson to be a very good approximation to the binomial
(b) the normal distribution to be a very good approximation to the binomial

Let n=50 and p=0.02

nvalue=50;
pvalue=0.02;
mean=nvalue*pvalue;
var=nvalue*pvalue*(1-pvalue);

plot1=BarChart[Table[{PDF[BinomialDistribution[nvalue,pvalue],k],k},
	{k,0,20}],PlotRange->{Automatic,{0,0.5}}]

plot2=BarChart[Table[{PDF[PoissonDistribution[mean],k],k},
	{k,0,20}],PlotRange->{Automatic,{0,0.5}},BarStyle->{RGBColor[0,1,0]}]

plot3=Plot[PDF[NormalDistribution[mean,Sqrt[var]],k],
	{k,0,20},PlotRange->{Automatic,{0,0.5}}]

Here is the normal distribution shown as a bar chart
(the value at 5 technically shows the proportion of the
normal distribution that lies between 4.5 and 5.5):

plot3=BarChart[Table[{
	Integrate[PDF[NormalDistribution[mean,Sqrt[var]],z],{z,k-0.5,k+0.5}],k},
	{k,0,20}],PlotRange->{Automatic,{0,0.5}},BarStyle->{RGBColor[0,0,1]}]

Show[plot1,plot2]

(You may want to reorder these plots to see the other plots in front.)

Show[plot1,plot3]

Let n=50 and p=0.1

nvalue=50;
pvalue=0.1;
mean=nvalue*pvalue;
var=nvalue*pvalue*(1-pvalue);

plot1=BarChart[Table[{PDF[BinomialDistribution[nvalue,pvalue],k],k},
	{k,0,20}],PlotRange->{Automatic,{0,0.5}}]

plot2=BarChart[Table[{PDF[PoissonDistribution[mean],k],k},
	{k,0,20}],PlotRange->{Automatic,{0,0.5}},BarStyle->{RGBColor[0,1,0]}]

plot3=Plot[PDF[NormalDistribution[mean,Sqrt[var]],k],
	{k,0,20},PlotRange->{Automatic,{0,0.5}}]

Here is the normal distribution shown as a bar chart
(the value at 5 technically shows the proportion of the
normal distribution that lies between 4.5 and 5.5):

plot3=BarChart[Table[{
	Integrate[PDF[NormalDistribution[mean,Sqrt[var]],z],{z,k-0.5,k+0.5}],k},
	{k,0,20}],PlotRange->{Automatic,{0,0.5}},BarStyle->{RGBColor[0,0,1]}]

Show[plot1,plot2]

(You may want to reorder these plots to see the other plots in front.)

Show[plot1,plot3]

Although the Normal distribution appears to provide a reasonable match to the Binomial in this case, a major difference is that the Normal can give negative values.  Using PDF[NormalDistribution[mean,Sqrt[var]],z], how could you determine the fraction of times that the Normal Distribution would give you a negative value?

Let n=50 and p=0.5

nvalue=50;
pvalue=0.5;
mean=nvalue*pvalue;
var=nvalue*pvalue*(1-pvalue);

plot1=BarChart[Table[{PDF[BinomialDistribution[nvalue,pvalue],k],k},
	{k,15,35}],PlotRange->{Automatic,{0,0.2}}]

plot2=BarChart[Table[{PDF[PoissonDistribution[mean],k],k},
	{k,15,35}],PlotRange->{Automatic,{0,0.2}},BarStyle->{RGBColor[0,1,0]}]

plot3=Plot[PDF[NormalDistribution[mean,Sqrt[var]],k],
	{k,15,35},PlotRange->{Automatic,{0,0.2}}]

Here is the normal distribution shown as a bar chart
(the value at 5 technically shows the proportion of the
normal distribution that lies between 4.5 and 5.5):

plot3=BarChart[Table[{
	Integrate[PDF[NormalDistribution[mean,Sqrt[var]],z],{z,k-0.5,k+0.5}],k},
	{k,15,35}],PlotRange->{Automatic,{0,0.2}},BarStyle->{RGBColor[0,0,1]}]

Show[plot1,plot2]

(You may want to reorder these plots to see the other plots in front.)

Show[plot1,plot3]

Questions

When n=50 and p=0.02, the normal approximation to the binomial was poor.
Based on the Means and Variances section and the graphs,
which of the following applies in this case:
   (a) normal distribution has a mean that is too high
   (b) normal distribution has a variance that is too high
   (c) normal distribution has a variance that is too low
   (d) normal distribution has a skew that is too low (too symmetric).

When n=50 and p=0.5, the Poisson approximation to the binomial was poor.
Based on the Means and Variances section and the graphs,
which of the following applies in this case:
   (a) Poisson distribution has a mean that is too high
   (b) Poisson distribution has a variance that is too high
   (c) Poisson distribution has a variance that is too low
   (d) Poisson distribution has a skew that is too low (too symmetric).

A rule of thumb that is often used is that the normal distribution
will approximate the binomial distribution well as long as n>32
and p is not very near 0 or 1.