Generating the Figures in Primer3 of Otto and Day

Clear["Global`*"]

Figure P3.2 - Bernoulli histogram

BarChart[{0.5,0.5},ChartLabels->{0,1},PlotRange->{{0.25,2.5},{-0.05,1}},Ticks->{False,Automatic}]

BarChart[{0.2,0.8},ChartLabels->{0,1},PlotRange->{{0.25,2.5},{-0.05,1}},Ticks->{False,Automatic}]

Figure P3.3 - Binomial distribution

label=Table[x,{x,0,6}];

BarChart[Table[PDF[BinomialDistribution[6,0.5],x],{x,0,6}],PlotRange->{Automatic,{-0.02,0.55}},ChartLabels->label,Ticks->{False,Automatic}]

BarChart[Table[PDF[BinomialDistribution[6,0.1],x],{x,0,6}],PlotRange->{Automatic,{-0.02,0.55}},ChartLabels->label,Ticks->{False,Automatic}]

BarChart[Table[PDF[BinomialDistribution[6,0.9],x],{x,0,6}],PlotRange->{Automatic,{-0.02,0.55}},ChartLabels->label,Ticks->{False,Automatic}]

Figure P3.4 - Geometric distribution

Note that Mathematica represents the geometric distribution as the number of failures that occur before the xth success.  Because we are interested in the total number of events that occur until the first success, including that success, we add one (changing the plot range from 0-14 to 1-15):

label=Table[x,{x,1,15}];

BarChart[Table[PDF[GeometricDistribution[0.5],x],{x,0,14}],PlotRange->{Automatic,{-0.03,0.51}},AspectRatio->0.4,ChartLabels->label,Ticks->{False,Automatic}]

BarChart[Table[PDF[GeometricDistribution[0.1],x],{x,0,14}],PlotRange->{Automatic,{-0.03,0.51}},AspectRatio->0.4,ChartLabels->label,Ticks->{False,Automatic}]

Figure P3.5 - Negative binomial distribution

Note that Mathematica represents the negative binomial distribution as the number of failures that occur before the nth success.  Because we are interested in the total number of events that occur until the nth success, including these successes, we add n:

label=Table[x,{x,4,18}];

BarChart[Table[PDF[NegativeBinomialDistribution[4,0.5],x],{x,0,14}],PlotRange->{Automatic,{-0.03,0.51}},AspectRatio->0.4,ChartLabels->label,Ticks->{False,Automatic}]

label=Table[x,{x,8,22}];

BarChart[Table[PDF[NegativeBinomialDistribution[8,0.5],x],{x,0,14}],PlotRange->{Automatic,{-0.03,0.51}},AspectRatio->0.4,ChartLabels->label,Ticks->{False,Automatic}]

Figure P3.6 - Poisson distribution

label=Table[x,{x,0,14}];

BarChart[Table[PDF[PoissonDistribution[0.75],x],{x,0,14}],PlotRange->{Automatic,{-0.03,0.51}},AspectRatio->0.4,ChartLabels->label,Ticks->{False,Automatic}]

label=Table[x,{x,0,14}];

BarChart[Table[PDF[PoissonDistribution[7.5],x],{x,0,14}],PlotRange->{Automatic,{-0.03,0.51}},AspectRatio->0.4,ChartLabels->label,Ticks->{False,Automatic}]

Figure P3.7 - Poisson distribution

Show[
BarChart[Table[PDF[PoissonDistribution[115.2],x],{x,80,150}],PlotRange->{Automatic,{0,0.05}},AspectRatio->0.4,BarSpacing->None,Ticks->{{{0,"80"},{5,"85"},{10,"90"},{15,"95"},{20,"100"},{25,"105"},{30,"110"},{35,"115"},{40,"120"},{45,"125"},{50,"130"},{55,"135"},{60,"140"},{65,"145"},{70,"150"}},Automatic}],
Plot[0.04,{x,0,15}],
Plot[0.04,{x,57,70}]

]

Export["FigP3-7.eps",%]

Sum[PDF[PoissonDistribution[115.2],x],{x,0,95}]

Sum[PDF[PoissonDistribution[115.2],x],{x,137,Infinity}]

%+%%

Figure P3.9 - Uniform distribution

ListPlot[{{0.25,0},{0.25,1/1.5},{1.75,1/1.5},{1.75,0}},Joined->True,PlotRange->{{0,2},{0,2/1.5}},AspectRatio->0.4,Ticks->{{0,{0.25,"min"},{1.75,"max"}},{0,{1/1.5,"1/(max-min)"}}}]

Figure P3.10 - Expoential distribution

Plot[PDF[ExponentialDistribution[0.5],x],{x,0,14},PlotRange->{{0,14.5},{0,0.51}},AspectRatio->0.4]

Plot[PDF[ExponentialDistribution[0.1],x],{x,0,14},PlotRange->{{0,14.5},{0,0.51}},AspectRatio->0.4]

label=Table[x,{x,1,14}];

Show[
BarChart[Table[PDF[GeometricDistribution[0.5],x],{x,0,14}],PlotRange->{{0,14.5},{-0.03,0.51}},AxesOrigin->{0,0},AspectRatio->0.4,ChartLabels->label,Ticks->{False,Automatic}],
Plot[PDF[ExponentialDistribution[0.5],x],{x,0,14},PlotStyle->Thickness[0.01],
PlotRange->{{0,14.5},{0,0.51}},AxesOrigin->{0,0},AspectRatio->0.4]
]

Show[
BarChart[Table[PDF[GeometricDistribution[0.1],x],{x,0,14}],PlotRange->{{0,14.5},{-0.03,0.51}},AxesOrigin->{0,0},AspectRatio->0.4,ChartLabels->label,Ticks->{False,Automatic}],
Plot[PDF[ExponentialDistribution[0.1],x],{x,0,14},PlotStyle->Thickness[0.01],PlotRange->{{0,14.5},{0,0.51}},AxesOrigin->{0,0},AspectRatio->0.4]
]

Figure P3.11 - Gamma distribution

Note that Mathematica represents Gamma distributions differently than is standard in biology.

?GammaDistribution

Mean[GammaDistribution[b,1/a]]

Variance[GammaDistribution[b,1/a]]

Show[
Plot[PDF[GammaDistribution[0.5,4],x],{x,0,14},PlotRange->All],
Plot[PDF[GammaDistribution[1,2],x],{x,0,14},PlotRange->All],
Plot[PDF[GammaDistribution[2,1],x],{x,0,14},PlotRange->All],
PlotRange->{{0,5},{0,0.75}},AspectRatio->0.4]

Show[
Plot[PDF[GammaDistribution[0.5,8],x],{x,0,14},PlotRange->All],
Plot[PDF[GammaDistribution[1,4],x],{x,0,14},PlotRange->All],
Plot[PDF[GammaDistribution[2,2],x],{x,0,14},PlotRange->All],
PlotRange->{{0,10},{0,0.375}},AspectRatio->0.4]

Figure P3.12 - Normal distribution

?NormalDistribution

Show[
Plot[PDF[NormalDistribution[1,0.5],x],{x,-5,8},PlotRange->All],
Plot[PDF[NormalDistribution[1,1],x],{x,-5,8},PlotRange->All],
Plot[PDF[NormalDistribution[1,2],x],{x,-5,8},PlotRange->All],
PlotRange->{{-5,8},{0,1}}]

Export["FigP3-12a.eps",%];

?Plot3D

Plot3D[{z=(1/(2 π σx σy Sqrt[1-(ρ^2)])) Exp[-((((((x-μx)/(σx))))^2)+(((((y-μy)/(σy))))^2)-2 ρ (((x-μx)/(σx))) (((y-μy)/(σy))))/(2 (1-(ρ^2)))];z}/.{μx->1,μy->3,σx->1,σy->1,ρ->0.8},{x,-2,4},{y,0,6},PlotRange->{0,.3}]

Export["FigP3-12b.eps",%];

Figure P3.13 - Normal and binomial distributions

plotnorm=Plot[PDF[NormalDistribution[3,Sqrt[30 0.1 0.9]],x],{x,-1,14},PlotRange->{{-1,10},{0,0.3}},AspectRatio->0.4,PlotStyle->Thickness[0.01]]

The following is the same, but using (x-1) as the variable so that it lines up with the BarChart, which places the first bar at “1” by default:

plotnorm=Plot[PDF[NormalDistribution[3,Sqrt[30 0.1 0.9]],{x-1}],{x,0,15},PlotRange->{{0,11},{0,0.3}},AspectRatio->0.4,PlotStyle->Thickness[0.01],Ticks->{Table[{i,i-1},{i,0,10}]}]

Show[
BarChart[Table[PDF[BinomialDistribution[30,0.1],x],{x,0,10}],ChartLabels->Table[x,{x,0,10}],PlotRange->{{0,11},{-0.02,0.3}},AxesOrigin->{1,0},AspectRatio->0.4,BarSpacing->0.5,Ticks->{False,Automatic}],
plotnorm
]

Figure P3.14 - Log-normal distribution

Show[
Plot[Exp[-(Log[y]-m)^2/(2 s^2)]/Sqrt[2 Pi y^2 s^2]/.m->1/.s->0.5,{y,0,20},PlotRange->All],
Plot[Exp[-(Log[y]-m)^2/(2 s^2)]/Sqrt[2 Pi y^2 s^2]/.m->1/.s->1,{y,0,20},PlotRange->All],
Plot[Exp[-(Log[y]-m)^2/(2 s^2)]/Sqrt[2 Pi y^2 s^2]/.m->1/.s->2,{y,0,20},PlotRange->All],
PlotRange->{{-1,8},{0,0.65}}]

(E^(m+((s^2)/2)))/.m->1/.s->{0.5,1,2}//N

Simplify[PDF[LogNormalDistribution[m,s],x],{x>0}]

Simplify[Integrate[Exp[-(Log[y]-m)^2/(2 s^2)]/Sqrt[2 Pi y^2 s^2],{y,0,Infinity}],{s>0}]

The mean equals:

Simplify[Integrate[y Exp[-(Log[y]-m)^2/(2 s^2)]/Sqrt[2 Pi y^2 s^2],{y,0,Infinity}],Re[(s^2)]>0]

The variance equals:

Simplify[Integrate[y^2 Exp[-(Log[y]-m)^2/(2 s^2)]/Sqrt[2 Pi y^2 s^2],{y,0,Infinity}]-((E^(m+((s^2)/2))))^2,Re[(s^2)]>0]

?LogNormalDistribution

Mean[LogNormalDistribution[m,s]]

Note that the intercept with the vertical axis (y=0) always occurs at zero:

Simplify[Limit[Exp[-(Log[y]-m)^2/(2 s^2)]/Sqrt[2 Pi y^2 s^2],y->0],s>0]

Figure P3.15 - Beta distribution

Mean[BetaDistribution[a,b]]

Variance[BetaDistribution[a,b]]

Show[
Plot[PDF[BetaDistribution[1,1],x],{x,0,1},PlotRange->All],
Plot[PDF[BetaDistribution[0.2,0.2],x],{x,0,1},PlotRange->All],
Plot[PDF[BetaDistribution[5,5],x],{x,0,1},PlotRange->All],
PlotRange->{{0,1},{0,5}},AspectRatio->0.4]

Show[
Plot[PDF[BetaDistribution[1,5],x],{x,0,1},PlotRange->All],
Plot[PDF[BetaDistribution[0.2,1],x],{x,0,1},PlotRange->All],
Plot[PDF[BetaDistribution[5,25],x],{x,0,1},PlotRange->All],
PlotRange->{{0,1},{0,10}},AspectRatio->0.4]

Show[
Plot[PDF[BetaDistribution[5,1],x],{x,0,1},PlotRange->All],
Plot[PDF[BetaDistribution[1,0.2],x],{x,0,1},PlotRange->All],
Plot[PDF[BetaDistribution[25,5],x],{x,0,1},PlotRange->All],
PlotRange->{{0,1},{0,10}},AspectRatio->0.4]

Figure P3.16 - A new distribution

BarChart[{F,(1-F)F,(1-F)^2}/.F->0.8,ChartLabels->{1,2,3},PlotRange->{{0.25,3.5},{-0.05,1}},AspectRatio->0.5,Ticks->{False,Automatic}]

BarChart[{F,(1-F)F,(1-F)^2}/.F->0.2,ChartLabels->{1,2,3},PlotRange->{{0.25,3.5},{-0.05,1}},AspectRatio->0.5,Ticks->{False,Automatic}]

The mean equals:

1*F+2*(1-F)*F+3*(1-F)^2//Simplify

The variance equals:

(1-mean)^2*F+(2-mean)^2*(1-F)*F+(3-mean)^2*(1-F)^2/.mean->%//Factor

Supplementary Figure - Beta and Bayes

The prior probability (Beta) times the probability of observing the data (binomial), integrated over all possible values of p:

Simplify[PDF[BetaDistribution[a,b],p] PDF[BinomialDistribution[n,p],k],{0≤k≤n,0<p<1}]/Integrate[Simplify[PDF[BetaDistribution[a,b],p] PDF[BinomialDistribution[n,p],k],{0≤k≤n,0<p<1}],{p,0,1}]

Simplify[%,{Re[a+k]>0,Re[b-k+n]>0}]

Using the relation between gamma and sine functions, the result is a beta distribution with new parameters a+k and b+n-k:

Simplify[PDF[BetaDistribution[a+k,b+n-k],p] ,{0<p<1}]

FullSimplify[%/%%]

Show[
Plot[PDF[BetaDistribution[1,9],x],{x,0,1},PlotRange->All],
Plot[PDF[BetaDistribution[1+30,9+20],x],{x,0,1},PlotRange->All],
Plot[PDF[BetaDistribution[1+300,9+200],x],{x,0,1},PlotRange->All],
PlotRange->{{0,1},{0,20}}]