Binomial Distribution

from Statistical Physics using Mathematica
© James J. Kelly, 1996-2001

Initialization

Defaults and packages

ClearAll["Global`*"];
Off[General::spell, General::spell1]

$DefaultFont = {"Times", 12};
$TextStyle = {FontFamily -> "Times", 
          FontSize -> 12};

Needs["Graphics`Master`"]

Establish rules for manipulating logarithms and factorials

Stirling formulas

These Stirling formulas are derived in the notebook stirling.nb.

\(logStirling1 = {Log[\(n_!\)] -> 
          n\ Log[n] - n}; 
\(logStirling2 = {Log[\(n_!\)] -> 
          n\ Log[n] - n + 1\/2\ Log[2\ π\ n]}; \)\)
logStirling3 = {Log[n_!] -> 
          n Log[n] - n + (1/2) Log[2 π n] + 
            (1/(12 n))};

Rules for simplify Gamma functions

gammaRule = {Gamma[1 + n_] -> n!, 
        Gamma[n_] -> (n - 1)!};

Logarithm rules

These rules, which permit logarithms to be expanded and combined easily, are used to direct the simplification process toward the final form desired; they are provided for purely cosmetic effect.

\(logExpand = {Log[a_\ b_] -> Log[a] + Log[b], 
        Log[a_\/b_] -> Log[a] - Log[b], Log[a_\^n_] -> n\ Log[a], 
        Log[Rational[1, n_]] -> \(-Log[n]\)}; 
    logCombine = {c_. \ Log[a_] + c_. \ Log[b_] -> c\ Log[a\ b], 
        c_. \ Log[a_] - c_. \ Log[b_] -> c\ Log[a\/b], 
        Log[a_] + Log[b_] -> Log[a\ b], 
        Log[a_] - Log[b_] -> Log[a\/b]}; \)

Definition of the binomial coefficient and probability distribution

g[N_, n_] := ((N!)/(n! (N - n)!))

probability[n_, N_, p_] := 
      g[N, n] (p^n) (((1 - p))^(N - n))

Definition of the Gaussian probability distribution

gaussian[n_, m_, σ_] := 
      Exp[-((((n - m))^2)/(2  (σ^2)))]/((Sqrt[2  π ]) σ)

Evaluate accuracy of Stirling formulas

Accuracy for Log[n!]

The percent deviation between true and approximate evaluations of Log[n!] is evaluated for several Stirling formulas of successively higher order.  These deviations are plotted as functions of n below.

\(logN[n_] = Log[\(n!\)]; 
\(logN1[n_] = 
      Log[\(n!\)] /. logStirling1; \)\)
logN2[n_] = Log[n!] /. logStirling2; 
logN3[n_] = Log[n!] /. logStirling3;

errorLogN1[n_] := N[100 (logN1[n]/logN[n] - 1)]
errorLogN2[n_] := N[100 (logN2[n]/logN[n] - 1)]
errorLogN3[n_] := N[100 (logN3[n]/logN[n] - 1)]

LogPlot[
      Evaluate[Abs[{errorLogN1[n], errorLogN2[n], errorLogN3[n]}]], {n, 5, 
        100}, Frame -> True, 
      FrameLabel -> {"N", "percent error"}, 
      PlotLabel -> "log N!"];

The percent error becomes quite small for relatively small values of N.

Accuracy for log of binomial coefficients

The percent deviation between true and approximate evaluations of the binomial coefficient is evaluated for several Stirling formulas of successively higher order.  These deviations are plotted as functions of n below.

\(logG[n_, m_] = 
      Log[g[n, m]]; 
\(logG1[n_, m_] = \(Log[g[n, m]] //. logExpand\) /. 
        logStirling1; \)\)
logG2[n_, m_] = Log[g[n, m]] //. logExpand /. 
        logStirling2; 
logG3[n_, m_] = Log[g[n, m]] //. logExpand /. 
        logStirling3;

errorLogG1[m_] := 
        N[100 (logG1[2  m, m]/logG[2  m, m] - 1)];
errorLogG2[m_] := 
        N[100 (logG2[2  m, m]/logG[2  m, m] - 1)];
errorLogG3[m_] := 
      N[100 (logG3[2  m, m]/logG[2  m, m] - 1)]

LogPlot[
      Evaluate[Abs[{errorLogG1[n], errorLogG2[n], errorLogG3[n]}]], {n, 5, 
        100}, Frame -> True, 
      FrameLabel -> {"n/2", "percent error"}, 
      PlotLabel -> "log multiplicity"];

Analytic properties of binomial coefficients

In this section we examine some of the properties of the binomial coefficient using both our definition and the intrinsic Mathematica function.

(∑_(n = 0)^N) Binomial[N, n]

(∑_(n = 0)^N)(((-1))^n) Binomial[N, n]

(∑_(n = 0)^N) g[N, n]  /. gammaRule

(∑_(n = 0)^N)(((-1))^n) g[N, n]

Binomial distribution

Display binomial distribution

values = {N -> 50, p -> 0.7};

nmax = N /. values; 
    distribution = 
      Table[Evaluate[probability[n, N, p] /. values], {n, 1, 
          nmax}];

\(units[xmin_, xmax_] := Range[0, Floor[xmax], 10]; 

    BarChart[distribution, PlotRange -> {0, Max[distribution]}, 

      Ticks -> {units, Automatic}]; \)

Verify normalization and obtain moments using raising operator

Direct method

The leading moments of the binomial distribution are evaluated directly.  Mathematica must be instructed to perform the desired algebraic simplifications because it does not automatically assume that the variables should have values appropriate to probability distributions.   Nevertheless, it is not difficult to recover the familiar moments.

((∑_(n = 0)^N) probability[n, N, p]  // Simplify // 
          PowerExpand)/.gammaRule

nMean = ((∑_(n = 0)^N) n probability[n, N, p]  // 
              Simplify // PowerExpand) /. N! -> 
          N Gamma[N]

n2Mean  = ((∑_(n = 0)^N) (n^2) probability[n, N, 
                    p]  /. N! -> N Gamma[N])  // Simplify // 
        PowerExpand

variance = Simplify[n2Mean - (nMean^2)]

Raising operator method

The raising operator method is based upon the observation that the binomial expansion of (((p+q))^N) , where p+q=1, contains N terms which can each be interpreted as a probability.  If n is the exponent of p, then the term involving (p^n) is the probability for n successes in N trials.  The sum is normalized to unity and each successive application of the raising operator p (∂_p)x generates a higher moment of the distribution.

sumProb = (((p + q))^N)

raise[x_, p_] := p (∂_p) x

raise[sumProb, p] /. p + q -> 1

(raise[raise[sumProb, p], p]/.{p+q->1})//Collect[#, {p, N}]&

This method gives the same results for the moments of the distribution and is a useful method for analyzing more complicated distribution functions.

Gaussian distribution

Derive Gaussian approximation to binomial distribution

Next we develop a Gaussian approximation to the binomial distribution that is valid when N is large and p is not too close to either 0 or 1.  This is done by performing a Taylor series expansion of the Log of the binomial distribution to second order in n and applying the appropriate limiting conditions.  The lowest order Stirling approximation is used to simplify the resulting expressions.

logProb = Log[ probability[n,N,p] ] //.logExpand/.logStirling1

First we determine the peak of the distribution in the limit that both n and N are very large compared to unity.  This is done by examining the first derivative of logProb.

term1 = D[logProb,n]//.logCombine//Simplify

sol1 = Solve[term1==0, n][[1]]

Therefore, to this level of approximation, the most probable value of n is equal to the mean of the binomial distribution.  At this level of approximation, the constant term becomes:

term0=(logProb/.sol1//.logCombine//Simplify)/.logCombine//Simplify

Next, we evaluate the second-order contribution to the Taylor expansion.

term2  =  D[term1, n]  /. sol1  // Simplify

We now recognize this result as the reciprocal of the variance of the distribution.

term2 variance //Simplify

Combining these terms and exponentiating, we obtain an unnormalized Gaussian distribution.  The normalization constant is evaluated by integration.  The range of integration is extended to (-∞, ∞) because both N and N p are very large.  Note that in order to obtain the normalization without this step, we would have had to use the logStirling2 instead of logStirling1.

arg = -((((x - m))^2)/(2 (σ^2))); 
unnormalizedApprox = (E^arg)

normalizationFactor = 1/Integrate[unnormalizedApprox, {x,-∞,∞},Assumptions->{m∈Reals,σ>0}]

Therefore, we can approximate the binomial distribution by the familiar Gaussian distribution by identifying the most probable value and the variance with those of the binomial distribution.

gaussianApprox[n_,N_,p_] := gaussian[n, N p, Sqrt[N p(1-p)]]

Compare Gaussian and binomial distributions

We compare a Gaussian distribution to an asymmetric binomial distribution, where both distributions have the same mean and variance.  One finds that the Gaussian approximation is already pretty good even for rather asymmetric distributions with relatively small N.  The difference between the two distributions can be difficult to discern on a linear plot, but is more readily visible using a semilog plot.  By modifying the values below when can examine a wide variety of cases.

\(values = {N -> 50, p -> 0.7}; 

    nmax = N /. values; \)

points = 
      Table[{n, 
          Evaluate[probability[n, N, p] /. values]}, {n, 0, 
          nmax}];

plotPoints = 
      ListPlot[points, PlotStyle -> PointSize[0.01], 
        DisplayFunction -> Identity];

histogram = 
      BarChart[Table[
          Evaluate[probability[n, N, p] /. values], {n, 1, 
            nmax}], Ticks -> {units, Automatic}, 
        DisplayFunction -> Identity];

LogPlotPoints = 
      LogListPlot[points, PlotStyle -> PointSize[0.01], 
        DisplayFunction -> Identity];

PlotGaussianApprox = 
      Plot[Evaluate[gaussianApprox[n, N, p] /. values], {n, 
          0, nmax}, 
        PlotStyle -> {{}, Dashing[{0.02, 0.02}], 
            Dashing[{0.04, 0.01, 0.02, 0.01}]}, 
        DisplayFunction -> Identity];

LogPlotGaussianApprox = 
      LogPlot[Evaluate[
          gaussianApprox[n, N, p] /. values], {n, 0, nmax}, 
        PlotStyle -> {{}, Dashing[{0.02, 0.02}], 
            Dashing[{0.04, 0.01, 0.02, 0.01}]}, 
        DisplayFunction -> Identity];

Show[{plotPoints, PlotGaussianApprox, histogram}, 
      PlotRange -> {{0, nmax}, {0, 0.2}}, Frame -> True, 
      FrameLabel -> {"n", "probability"}, 
      PlotLabel -> "gaussian approximation", 
      DisplayFunction -> $DisplayFunction];

Show[{LogPlotPoints, LogPlotGaussianApprox}, Frame -> True, 
      FrameLabel -> {"n", "probability"}, 
      PlotLabel -> "gaussian approximation", 
      DisplayFunction -> $DisplayFunction];

Poisson distribution

There are many physical processes, such as the decay of radioactive nuclei, for which the number of trials N is very large but the probability for the event of interest is very small.  Under those circumstances the binomial distribution can become so asymmetric that the Gaussian distribution does not provide a suitable approximation.  Fortunately, one can easily derive a more useful approximation known as the Poisson distribution.

Derive as limit of binomial distribution

It is useful to introduce λ = N p as the mean value of n.

logPoisson  =  (logProb /. {p -> (λ/N)}) // 
        Limit[#, N -> ∞] &

Exp[logPoisson]//Simplify

It is customary to use Stirling's approximation in reverse to express the (n^(-n)) factor in terms of n!, although that is not really necessary.

poisson = 
      Exp[logPoisson] /. {n Log[n] -> Log[n!] + n} // 
        Simplify

Moments

(∑_(n = 0)^∞) Evaluate[poisson]

nMean = (∑_(n = 0)^∞) Evaluate[n poisson]

n2Mean = (∑_(n = 0)^∞) Evaluate[
          (n^2) poisson]

variance = Simplify[n2Mean - (nMean^2)]

Compare Poisson and Gaussian distributions

Since these distributions differ most for small numbers, we compare Gaussian and Poisson distributions with equal means and variances for a small value of the mean.

\(values2 = {λ -> 2}; 

    nmax = 5\ λ /. values2; \)

poissonDistribution = 
      Table[Evaluate[poisson /. values2], {n, 0, 
          nmax}];

poissonPlot = 
      BarChart[poissonDistribution, 
        PlotRange -> {0, 1.1 Max[poissonDistribution]}, 
        DisplayFunction -> Identity];

\(σ = \@λ /. values2; 
    gaussianPlot = 
      Plot[Evaluate[
          gaussian[
              n - 1, λ, σ] /. values2], {n, \(-nmax\), nmax}, Frame -> True, 
        FrameLabel -> {"n+1", "probability"}, 
        PlotRange -> {0, 1.1\ Max[poissonDistribution]}, 
        DisplayFunction -> Identity]; \)

units[xmin_, xmax_] := Range[Floor[xmin], Floor[xmax], 2]; 
    Show[{gaussianPlot, poissonPlot}, 
      PlotRange -> {{0, nmax}, {0, 1.1 Max[poissonDistribution]}}, 
      Ticks -> {units, Automatic}, 
      DisplayFunction -> $DisplayFunction];

Technical notes

1) There is an annoying complication when using BarChart to produce a histogram for a discrete probability distribution.  Although there is a finite probability for n = 0 successes, BarChart begins its numbering at n = 1 because it simply counts bins.  Therefore, to accurately compare the discrete and thecontinuous distributions, we must plot the continuous distribution as a function of n + 1.  There seems to be no easy way to adjust the binning employing by BarChart.  Although it should not be difficult to create a routine which produces the desired histogram plot, the expedient employed here suffices for our purposes.

Exercises

Compare Gaussian, Poisson, and binomial distributions

a) Write a module which graphically compares the Gaussian and binomial distributions for arbitrary choices of parameters.  In other words, the arguments should be n, N, and p and the output should be linear and semilog plots that compare the distributions.  Use your module to survey various interesting cases such as small N, symmetric versus asymmetric distributions, etc.

b) Repeat for the Poisson distribution.

Refined derivation of Gaussian distribution

Modify the derivation of the Gaussian distribution from the binomial distribution by using logStirling2.  Note that more care will be needed in enforcing the limiting conditions.  Show that this refined procedure gives a properly normalized distribution without applying the normalization factor "by hand".

Radioactive decay

Uranium-238 decays by alpha emission with a half-life of 4.47 × (10^9) years.  Assume that the probability for any nucleus to decay in a specified time interval Δt is independent of time and of the state of any other nucleus in the sample.  Imagine dividing the observation time into N equal time intervals which are so small that the probability of observing more than one decay is negligible  The mass of our sample is 1 mg and the atomic weight is 238.  
	a) Calculate the average number of alpha particles emitted per second by this sample.
	b) Plot the distribution P (n) where n is the number of disintegrations observed during one second.
	c) Compare Poisson, binomial, and Gaussian distributions under these conditions.