12 - Probability and Statistics in Physics : Part I
▮
Week of April 1 - 5

Outline
12.1 Introduction and Terminology
12.2 The binomial distribution
12.3 The Poisson distribution
12.4 The normal distribution
12.5 An Exercise on Radioactive Decay

dir0="/home/stump/Desktop/phy102/"
SetDirectory[StringJoin[dir0,"/heap12"]]
FileNames[]

Probability and Statistics in Physics

Statistics is the branch of mathematics that is used in the analysis of experimental data.

For example, methods of “statistics” are necessary to account for the effects of experimental errors.

Also, “statistics” provides techniques for comparing data and theory.

Probability is a related branch of mathematics.

/12.1/ Introduction and Terminology

“Probability” involves the analysis of “random processes”.

Examples of random processes are flipping a coin, or throwing dice, or radioactive decays, or experimental errors in physical measurements.

The outcomes of a random process, repeated many times, gives a set of experimental data.

Let Y be a measurable quantity, i.e., a quantitative observable in an experiment.

Y is called a "random variable"  if the space of outcomes of the measurement is more than a single value.

Y may be continuous or discrete.
Discrete: flipping a coin; Y = heads or tails, two possible outcomes;
Discrete: rolling a die; Y = {1 2 3 4 5 6}, six possible outcomes;
Continuous: measuring a time interval; Y is a continuum(for an analog watch).

A “trial” consists of a single measurement of Y.

An “experiment” consists of many trials. Let n denote the number of trials in an experiment. The outcome of the experiment is a list of values {(Subscript[y, 1]), (Subscript[y, 2]), (Subscript[y, 3]), ..., (Subscript[y, n]) }.

The “mean” of an experiment is 1/n Σ (Subscript[y, i])   ≡ ⟨y⟩

The "variance”  is 1/n Σ ((Subscript[y, i]) – ⟨y⟩ )^2  = ⟨(y^2)⟩ – ((⟨y⟩)^2)

The “standard deviation” is SQRT [ variance ].

•••

Today we are going to study the three probability distributions that are the most important in physics. Also, we will see how probability distributions are implemented in Mathematica.

Mathematica has many built-in probability distributions, along with numerous commands that can be used with any specific probability distribution. We will study these three distributions:

BinomialDistribution[n,p]     (depends on two parameters;
n = the number of trials, p = the probability of “success” in one trial)

PoissonDistribution[μ] (depends on a single parameter;  μ  is the mean value of many measurements)

NormalDistribution[μ,σ] (depends on two parameters; μ = mean and σ = standard deviation)

12.2 The Binomial Distribution;   B( n,p )

The binomial distribution B( n,p ) is a probability distribution associated with a binary process.

Consider a random process, with observable X, in which each trial has only two possible outcomes. In other words the domain of X is D = {0, 1}. We’ll use the following language: the possible outcomes are  x = 1 (called “success”) or x = 0 (called “failure”). This process is called a binary process.

Now assume that for any single measurement, the probability of success = p and the
probability of failure = q = 1 – p.

What is the probability of  obtaining r successes in n trials?
The answer to this question is the binomial distribution function;
     (Subscript[P, B])( n,p ; r ) = the probability of getting r “successes” in n trials
Note that r is a variable with domain { 0 1 2 3 ... n } ; n and p are fixed parameters.

We are assuming that the probability of success (p) does not change from one trial to the next, and does not depend on previous or future trials. Then the probability of getting success (i.e., x = 1)  r times in n trials is

(Subscript[P, B])( n,p ; r) = (({{n}, {r}})) (p^r) (((1–p))^(n–r))    = the binomial probability distribution function
where
     (({{n}, {r}})) = ((n!)/(r! (n–r)!))      = the binomial coefficient

Again, p is the probability of success in a single trial.

Example.
A good example of the binomial distribution is flipping a coin. Say
     “heads” ⟺ success ⟺ x = 1;
     “tails” ⟺ failure ⟺ x = 0.
For an unbiased coin, p = q = (1/2) . The probability of getting “heads” r times after flipping the coin 20 times is (Subscript[P, B])( 20 , 0.5 ; r).

Some useful Mathematica commands

In Mathematica we specify the binomial distribution with n trials and p = probability of success in one trial, by
     dist  =  BinomialDistribution[n,p] .
That is, we specify the name of the probability distribution that we want to use, and  call it by some expression name; here I’ve called it ‘dist’. The values n and p will be chosen somewhere in the program. For example, for 20 coin flips with an unbiased coin I might use a replacement rule to specify the values,
    dist  =  BinomialDistribution[n,p]  /. { n -> 20, p -> 0.5 }

●● The probability distribution function (p.d.f.) with respect to the variable r is obtained with the command PDF.

?PDF

●● The mean and standard deviation. The probability distribution PDF[dist,r] depends on a variable r. We will be interested in the distribution mean and standard deviation of the variable r.

the “mean” of the p.d.f.  is defined by ...
    = the mean of r 
    = ⟨r⟩   ≡   (Underscript[Σ, r values]) r P( n,p ; r )

the “standard deviation” of the p.d.f.  is defined by ...
    = SQRT [ ⟨ r – ⟨r⟩ ⟩^2 ]
    = SQRT [  ⟨(r^2)⟩ – ((⟨r⟩)^2) ]
    =  the “width” of the distribution of r values.

Mathematica commands to obtain the mean and standard deviation are Mean and StandardDeviation.

?Mean

?StandardDeviation

Exercises

EXERCISE  1 : flipping an unbiased coin
In Mathematica the probability distribution function for the binomial distribution is defined by

     (Subscript[P, B])[r_]    =    PDF [ BinomialDistribution [n , p] , r]  ] 

where
n = the number of trials (Note: we can’t use N because it is reserved.)
p = the probability of success in one trial,
r = the number of successes.
→→→
     Or, we could program it this way:
     dist = BinomialDistribution [n , p]
      (Subscript[P, B])[r_]   = PDF [ dist, r] 
←←←
(A) Let n = 20 trials and p = 0.5. Plot P(r) versus r for all possible values of r.
Since r has only discrete values, r ∈ {0 1 2 3 ... 20}, plot the function with the command DiscretePlot (see the help file if necessary).

(B) Make a bar chart. (See ?BarChart for instructions.)

(C) Calculate the mean of the number of “heads” from the defining equation, mean = Σ r P(r); the sum is over all values in the domain.
Also execute Mean[BinomialDistribution [n , p]  ] for arbitrary n and p.

(D) Calculate the standard deviation of the number of “heads” from the defining equation, sd = Σ(((r-mean))^2) P(r).
Also execute StandardDeviation[BinomialDistribution [n , p] ] for arbitrary n and p.

s1

Remove["Global`*"]
bd=BinomialDistribution[n,p]
{Mean[bd],StandardDeviation[bd]}
% /. {n->20,p->0.5}

dist=bd /.  {n-> 20,p->1/2}
f[r_]=PDF[dist,r]
DiscretePlot[f[r],{r,0,20}]
list=Table[f[r],{r,0,20}]//N
BarChart[list]

EXERCISE  2 : flipping a biased coin

(A) Let n = 20 trials and p = 0.7. {I.e., x = 1 (“heads”) is more likely than x = 0 (“tails”).) Plot P(r) versus r for all possible values of r.
(Since r has only discrete values, r ∈ {0 1 2 3 ... 20}, plot the function with the command DiscretePlot (see the help file if necessary).

(B) Make a bar chart. (See ?BarChart for instructions.)

(C) Calculate the mean and standard deviation of the number of  “heads” for 20 trials with this biased coin.

s2

Remove["Global`*"]
bd=BinomialDistribution[n,p]
{Mean[bd],StandardDeviation[bd]}
% /. {n->20,p->0.7}

dist=bd /.  {n-> 20,p->7/10}
f[r_]=PDF[dist,r]
DiscretePlot[f[r],{r,0,20}]
list=Table[f[r],{r,0,20}]//N
BarChart[list]

12.3 The Poisson Distribution;  Poi(μ)

There are several ways to derive the Poisson distribution, Poi(μ).

First, the Poisson distribution is a limiting form of the binomial distribution.
Suppose in a binary process the probability p is very small ( i.e., p ⟶ 0 ) and the number of trials is very large ( i.e., n ⟶ ∞ ). We take these limits in such a way that the mean, µ = n▫p,  remains constant.

Then what is the probability of observing r “successes” ?  Call it (Subscript[P, P])(μ,r) .

Comments. We are taking the limit p ⟶ 0, so the probability of getting success in one trial is very small. However we are also taking the limit n ⟶ ∞, so we will get some successes during the whole experiment. In fact, we are keeping n▫p constant, = μ. So the mean number of successes will be μ.

The probability distribution for this limit is
   (Subscript[P, P])( μ , r ) = ((μ^r)/(r!)) (e^(–μ))  .  the Poisson distribution function

A process whose probability distribution is (μ^r) (e^(–u)) /r! is called a Poisson process. We’ll see later that radioactive decay is a Poisson process.

Properties of the Poisson distribution

◼ (Subscript[P, P])(μ,r) is the probability of getting r successes when the mean is μ and number of trials is very large.
◼ The domain of r  is   { 0,1,2,3,... } = Z = the set of nonnegative integers; the number of successes must be an integer; and there is no limit on the number of successes.
◼ The mean number of successes is μ ; and the standard deviation of r is (Sqrt[μ]).

◼ The exercise below will show a bar chart of the Poisson distribution with mean μ = 7.

EXERCISE  3 : The p.d.f. for “events” in a Poisson process with mean number = 7.
(“Event” means a successful observation.)

Consider a Poisson process with mean = 7. In any one experiment, the number of events could be any integer from 0 to ∞.   The probability distribution is obtained in Mathematica from the command

     (Subscript[P, p])[r_]  =  PDF[ PoissonDistribution[μ] , r ]

(A) Use DiscretePlot to plot the function from r = 0 to 20.

(B) Make a bar chart. (See ?BarChart for instructions.)

(C) Calculate the mean and standard deviation of the number of events for this probability distribution function.

s3

Remove["Global`*"]
cllist=Table[j,{j,0,20}]
list=Table[PDF[PoissonDistribution[7],r],{r,0,20}]//N
BarChart[list,PlotRange->All,
ChartLabels->Placed[cllist,Bottom]]
Mean[PoissonDistribution[7]]
StandardDeviation[PoissonDistribution[7]]

We will prove the following theorem ...

Theorem
The limit of the binomial distribution B(n,p), as n ⟶ ∞  and p ⟶ 0 with n▫p = μ fixed, is the Poisson distribution Poi(μ).

What we need to prove is that (Subscript[P, P])(μ,r) = limit (Subscript[P, B])( n,μ/n ; r ) as n ⟶ ∞ .

First, we can use Mathematica to produce a graphical proof.
That is not rigorous, but it is convincing.

● First, plot these two p.d.f.s  
(Subscript[P, B])( n,p ; r ) (in red) for n=10 and p = 0.7, and (Subscript[P, P])( μ ,r ) (in black) for μ = 7.

fB[r_]=PDF[BinomialDistribution[10,0.7],r];
fP[r_]=PDF[PoissonDistribution[7],r];
DiscretePlot[{fB[r],fP[r]},{r,0,20},
PlotStyle->{Red,Black},ExtentSize->{0.5,0.25},ImageSize->256]

● You see the difference between the two distributions.
Now plot (Subscript[P, B])( n,p ; r ) (in red) for n=20 and p = 0.35, and (Subscript[P, P])( μ ,r ) (in black) for μ = 7.

fB[r_]=PDF[BinomialDistribution[20,0.35],r];
fP[r_]=PDF[PoissonDistribution[7],r];
DiscretePlot[{fB[r],fP[r]},{r,0,20},
PlotStyle->{Red,Black},ExtentSize->{0.5,0.25},ImageSize->256]

● You should see that the two distributions are closer together.
Now plot (Subscript[P, B])( n,p ; r ) (in red) for n=40 and p = 0.175, and (Subscript[P, P])( μ ,r ) (in black) for μ = 7.

fB[r_]=PDF[BinomialDistribution[40,0.175],r];
fP[r_]=PDF[PoissonDistribution[7],r];
DiscretePlot[{fB[r],fP[r]},{r,0,20},
PlotStyle->{Red,Black},ExtentSize->{0.5,0.25},ImageSize->256]

● Now the two distributions are very close together. As n increases, B( n,7/n ) approaches Poi[7], as claimed by the theorem. This sort of proves the theorem, although it is not complete nor rigorous. Another indication that the theorem is true comes from comparing the means and standard deviations.

EXERCISE  4
(A) Calculate the mean and standard deviation of B( n,7/n ) (the binomial distribution) for n = 10, 20, 40, 100.
(B) Show that these numbers approach the mean and standard deviation of Poi(7).

Rigorous Proof of the Theorem

Start with the binomial distribution. The distribution function is
(Subscript[P, B])( n,p ; r ) = (({{n}, {r}})) (r^p) (((n-r))^(1-p)) = ((n!)/(r! (n–r)!))   (r^p) (((n-r))^(1-p)) 
We want mean n▫p = μ, so substitute p = μ/n.

Remove["Global`*"]
temp=PDF[BinomialDistribution[n,p],r]
PB=temp[[1,1,1]] (* picks out the function *)
PB=PB /. {p ->μ/n}
PB=(μ/n)^r*(1-μ/n)^(n-r)*Gamma[n+1]/Gamma[r+1]/Gamma[n-r+1]

Now note that I can write B in terms of 4 factors,

fac1=(μ^r)/Gamma[r+1]
fac2=(1-μ/n)^n
fac3=1/(n-μ)^r
fac4=n!/(n-r)!
fac4=Gamma[n+1]/Gamma[n-r+1]

Verify that fac1*fac2*fac3*fac4 = B.

test=PowerExpand[fac1*fac2*fac3*fac4]//FullSimplify
β=PowerExpand[B]//FullSimplify
Simplify[test/β,Assumptions -> {n>0,r>0}]

Now take the limit as n ⟶ ∞ .

Limit[fac1,n->∞]
Limit[fac2,n->∞]
fac3*fac4
Limit[fac3*fac4,n->∞]

The result is that the limit is (((μ^r) (e^(-μ)))/(r!)),
and that proves the theorem.

TESTS of the Poisson distribution (Subscript[P, P])(μ,r)

EXERCISE 5
(A) Show that the normalization condition is true:     (Underoverscript[Σ, (r=0), ∞]) (Subscript[P, P])(μ,r) = 1.
(B) Use Mathematica to calculate the expectation value of r.
(C) Use Mathematica to calculate the standard deviation of r.
(D) Plot (Subscript[P, P])(5,r), (Subscript[P, P])(10,r) and (Subscript[P, P])(20,r) superimposed on one graph. Use appropriate ranges for the axes.

EXAMPLE : Radioactive decays are a Poisson process

◼ A random process for which the probability distribution is the Poisson distribution is called a Poisson process.

◼ The Poisson distribution describes binary processes for which the single-trial probability of success is very small (p ⟶ 0  and the number of trials is very large (n ⟶ ∞) such that the expected number of successes is μ. The probability distribution function is

     (Subscript[P, P])(μ,r) = (((μ^r) (e^(–μ)))/(r!)) .

◼ Now consider this experiment. We have a sample that includes a radioactive material. Let N(t) = the number of unstable atoms in the sample at time t. The radioactive decays of the unstable atoms occur randomly with a decay rate Γ. This is the decay rate, so the probability for any one of these atoms to decay during a short time interval from time t to time t + δt is Γ δt.

◼ The expected number of decays in the sample during the time interval δt is (Subscript[δN, d]) = N(t) (Γ δt). Therefore N(t + δt) = N(t) – (Subscript[δN, d]) = (1 – Γ δt ) N(t). Or, we may write a differential equation,

(dN/dt) =  –ΓN

The solution is N(t) = N(0) (e^(–Γ t))     (the exponential decay law).

◼ Now we ask, what is the probability to observe r decays during a time interval Δt (which is not necessarily small)? Subdivide Δt into n segments of size Δ/n where n is large. We’ll define a “success” to be a decay that occurs during the time interval from t to t + δt. So we effectively have a binary process. The number of trials is n = Δt/δt. The probability of success for any one trial is p = Γ δt.  Now let n be large and p be small, by taking δt ≪Δt.

◼ We have the conditions stated in the Theorem:
a binary process, with n -> ∞ and p -> 0, and the expected number of successes
is μ = n ▫ p = (Δt/δt) ▫ (Γ δt)  = Γ Δt.    The process is a Poisson process.

◼ Now to answer the question in the paragraph above...
The probability to observe r decays during the time interval Δt is (μ^r) (e^(–μ)) /r! , where μ = Γ Δt.   Remember, Γ is the mean decay rate for one trial. μ is dimensionless, and Γ has units (sec^(–1)).

12.4 The Gaussian or Normal Distribution

The normal distribution is the most familiar distribution in all the sciences. Measurement errors in an experiment, such as instrumental errors, are generally described by the Gaussian distribution. And even if the Gaussian distribution is not strictly correct, it often provides a good approximation to the true governing distribution.

The normal distribution is a continuous, symmetric distribution whose probability density function is

(Subscript[P, n])( μ,σ ; x ) =  K  (e^(–(((x-μ))^2)/(2(σ^2))))    (7)

where K is a normalization constant.

The name of the distribution in Mathematica is NormalDistribition[μ,σ].

?NormalDistribution

It can be shown that the Poisson distribution can be approximated by a normal distribution if μ is large. If μ is greater than about 20, the Poisson distribution can be approximated by a normal distribution with mean μ and standard deviation (Sqrt[μ]). This, and other features of the normal distribution, will be studied in a later lesson.

EXERCISE 6.
Let x be the variable in a normal distribution.
(A) The domain of x is from –∞ to ∞ . Determine the normalization constant K in equation (7). (Use Mathematica to evaluate the integral for arbitrary μ and σ.)
(B) Calculate the mean and standard deviation for the normal distribution with parameters μ and σ.
(C) Plot the Gaussian distribution for {μ,σ} = {0,1} and {1,2} and {2,3}.
(D) The full width at half maximum (FWHM) for the Gaussian distribution is
FWHM = A σ, where A is a numerical constant. Determine A, both exactly and approximately to 3 decimal places.

FWHM = 2 σ (Sqrt[2 ln 2])

s6

Assuming [s>0,Integrate[Exp[-(x-m)^2/(2s^2)],{x,-Infinity,Infinity}]]

Remove[μ,σ]
Mean[NormalDistribution[μ,σ]]
StandardDeviation[NormalDistribution[μ,σ]]

pdf[μ_,σ_,x_]=PDF[NormalDistribution[μ,σ],x];
Plot[{pdf[0,1,x],pdf[1,2,x],pdf[2,3,x]},{x,-10,10}]

test[σ_,ξ_]=PDF[NormalDistribution[0,σ],ξ]
Solve[test[σ,x]==(1/2)*test[σ,0],x]
% /. C[1]->0

12.5 An Exercise on Radioactive Decay

EXERCISE 7.
Consider a radioactive source of (^137)Cs . The half life of (^137)Cs is 27 years.
(A) Calculate the decay rate Γ in units of (sec^(–1)).
(B) Consider a sample with mass = 1 μg = 1.0 ×(10^(–9)) kg. Calculate the number of Cs atoms. (Hint the atomic mass of Cs is 137.)
(C) What is the probability to observe r decays for this sample in τ = 3 seconds?
(D) Calculate the mean and standard deviation for the number of decays observed in 3 seconds.

s7

(A)
N(t) = N0 exp(-Γt) = 1/2 N0
Γ (Subscript[t, half]) = ln 2

thalf=27*365.25*24*3600;
Γ=Log[2]/thalf

(*B*)
u=1.66*^-27;
mass=1.0*^-9;
NCs = mass / (137 u)

(C)
τ = 3
p = Γ▫t = mean number of decays per particle
μ = mean number from sample = n▫p = NCs▫Γ▫t 
prob for r decays = PDF[ PoissonDistribution[μ], r]

(*D*)
Remove[τ,p,μ]
τ=3
p=NCs*Γ
μ=p*τ
Sqrt[μ]

appx=NormalDistribution[μ,Sqrt[μ]]
Mean[appx]
StandardDeviation[appx]