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

Outline
12.1 Probability Distributions
12.2 The three most important probability distributions
12.3 Example : Radioactive decay of radium

SetDirectory["/home/stump/Desktop/phy102/lesson12-probability_and_statistics"]
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.

/12.1/ Probability Distributions

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

Examples of random processes are flipping a coin, or throwing dice, or radioactive decays.

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

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

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

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

Associated with the space of values of x, there a probability distribution P(x).

▮ If X is a discrete variable, then P(x) is the probability that one "trial" ( ≡ one measurement of X) will give the value x; the sum of P(x) must be 1;      (Underscript[Σ, x]) P(x) = 1 . In the case of a discrete variable, P(x) is the probability to get the value x.

▮ If X is a continuous variable, then P(x) dx is the probability that the measured value of X will lie in the interval (x,x+dx); the integral of P(x) must be 1;
∫  P(x) dx = 1. In the case of a continuous variable, P(x) is the probability density.

Now we need to review some mathematical analysis of probability distributions P(x).

● Cumulative distributions

We might want to know the probability for finding X between certain limits,
e.g., x1 < X < x2. This probability is given by the cumulative,
or integrated distribution,

For a continuous variable, 
     Probability{x1 < X < x2} = (∫_x1^x2) P(x) dx   (1)

For a discrete variable, the integral is replaced by a sum,
     Probability{x1 < X < x2} = (Underscript[Σ, x1 < x < x2]) P(x) dx   (2)

The probability distribution must be normalized to 1,

      (∫_D) P(x) dx = 1   D = domain of X values   (3)

In words, Eq. (3) states that the probability is 1 (i.e., 100 percent) that an observation will yield one of the possible outcomes; or, in other words, the observation must a value of X in the domain D.

● Expectation Values

The expectation value of the random variable X, or a random variable function f(X), is defined by

(4) E[X] = (∫_D)  x P(x) dx ;    E[f(X)] = (∫_D)  f(x) P(x) dx

So E[X] is the expected value of one measurement of X = the average of many measurements.

As usual the discrete case is analogous, but replacing the integrals by summations.

● The Mean and Variance and Standard Deviation

These quantities are familiar in data analysis.

◼ Sample Mean and Theoretical Mean
   sample mean =(1/N) (Underoverscript[Σ, (i=1), N]) (Subscript[x, i])  =  ⟨x⟩
   theoretical mean = E[x]

◼ Sample Variance and Theoretical Variance
   sample variance = (1/N) (Underoverscript[Σ, (i=1), N]) ( (Subscript[x, i]) – ⟨x⟩ )^2
   theoretical variance =  (∫_D) ( x – ⟨x⟩ )^(2  )P(x) dx  = E[(X^2)] - ((E[X])^2)

◼ Standard Deviation
Standard deviation = square root of the variance
   = root mean square deviation

/12.2/ The three probability distributions that are most important in Physics

In Section 1 we general reviewed properties of probability distributions.

Now we will study three examples — the three most important probability distributions in Physics:
	◼ the binomial distribution ;
	◼ the Poisson distribution ;
	◼ the Gaussian distribution (or, normal distribution) .
Also, we we see how probability distributions are implemented in Mathematica.

Mathematica has many built-in probability distributions, along with several commands that can be used with a probability distribution. We will study these distributions:

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

PoissonDistribution[μ] (depends on a single parameter; μ = mean)

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

2.1 The Binomial Distribution

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

Now, the questions is what is the probability of  obtaining r successes in n trials?

This probability is given by the binomial distribution, B(n;p).

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, p])(n,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 unbiassed coin, p = q = (1/2) . The probability of getting “heads” r times after flipping the coin 20 times is (Subscript[P, 1/2])(20,r).

Some useful Mathematica commands

To specify the binomial distribution with n trials and p = probability of success in one trial, set
     dist  =  BinomialDistribution[n,p]
in the commands below.

?PDF

?Mean

?StandardDeviation

Exercises

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

     P[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]
  P[r_] = PDF [ dist, r] 
←←←
(A) flipping an unbiased coin
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. Also execute Mean[BinomialDistribution [n , p]  ].

(D) Calculate the standard deviation. Also execute StandardDeviation[BinomialDistribution [n , p] ]

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.

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]

2.2 The Poisson Distribution

There are several ways to derive the Poisson distribution.

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

Then what is the probability of observing r “successes” ?

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

The probability distribution for this limit is
   P(r) = ((μ^r)/(r!)) (e^(–μ))  .  the poisson distribution function

A process whose probability is (μ^r) (e^(–u)) /r! is called a Poisson process.

Properties of the Poisson distribution

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

◼ The figure below shows a bar chart of the Poisson distribution with mean μ = 7.

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 need to prove this theorem ...

Theorem
The limit of the binomial distribution B(n;p), as n ⟶∞  and p ⟶ 0 with n*p = μ fixed is the Poisson distribution Poisson(μ).
A graphical proof of the theorem.

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

First, plot B(n;p) (in red) for n=10 and p = 0.7, and Poisson[μ] 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 distribution. Now plot B(n;p) (in red) for n=20 and p = 0.35, and Poisson[μ] 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 see the two distributions are closer together. Now plot B(n;p) (in red) for n=40 and p = 0.175, and Poisson[μ] 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 get larger and larger, B(n;7/n) approaches Poisson[7], as claimed by the theorem.

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]
B=temp[[1,1,1]]
B=B /. {p ->μ/n}
B=(μ/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 P(μ,r)

EXERCISE 3
(A) Show that the normalization condition is true:
(Underoverscript[Σ, (r=0), ∞]) 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 P(5,r), P(10,r) and 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

     P(μ,r) = (((μ^r) (e^(–μ)))/(r!))

◼Now consider an 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 the 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, 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 δt.
So we effectively have a binary process. The number of trials is n. 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: binary process, n ->∞, p ->0,
and the expected number of successes is μ = Γ Δt = Γ δt ▫n = n▫p. 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.   Here Γ is the mean decay rate for one unstable atom.

Information copied from an unidentified web site (no copyright given)

Poisson processes

A Poisson process is one in which occurrences are randomly distributed in time, space or some other variable with the number of occurrences in any non-overlapping intervals statistically independent.

For example, naturally occurring gamma rays detected in a scintillation detector are randomly distributed in time, or chocolate chips in a cookie dough are randomly distributed in volume. For simplicity, we will limit our discussion to occurrences or “events" randomly distributed in time. The “events” will be decays of unstable particles in a radioactive sample.

The average decay rate will be denoted Γ;
then in any time interval Δt the expected number of decays is Γ Δt ≡ μ.

While the expected number of decays for a given experiment need not be an integer, the number of number n actually observed must be. Moreover, due to the randomness of the decays, n may be more or less than Γ Δt.

We will prove that the probability to observe n events in the time interval Δt is

	P( n ) = (μ^n) exp(-μ ) / n! 

where μ = Γ Δt. This is the Poisson distribution.

Another, more rigorous derivation of the Poisson probability distribution

Let P( r ; t ) = the probability that r decays will occur in the time interval t.

Now consider P ( r+1 ; t ).

In order to have r+1 decays in time t, there must be ...

(i) r decays in time t', for some time t’ with t' < t ;
this has probability P(r ;t' ).

(ii) one decay during the infinitesimal time interval from t' to t' + dt' ; this has probability Γ dt' .

(iii) no decays in the interval from t' +dt'  to t;
this has probability P(0 ; t – t' – dt') which we may approximate by
P(0; t – t') treating dt' as infinitesimal.

Now, P( r+1; t) is the product of the three probabilities, integrated over all possible values of t', i.e., t’ from from 0 to t :

	P( r+1 ; t ) = ((∫_0)^t)  P( r  ;t' ) Γ dt' P( 0 ; t – t' ).

Now, we know that

	 P( 0 ; t – t' ) = exp{ – Γ (t – t' ) }
	 
from the law of exponential decay, so

	P( r+1 ; t ) = ((∫_0)^t)  P( r  ;t' ) Γ exp{ – Γ (t – t' ) }  dt' .

Now it is straightforward to show that the solution of this recursion relation is 

	P( r ; t ) = (μ^(r )) exp{–μ} / r!	where μ = Γ t ;
	
and this is the Poisson distribution.

An example of a Poisson process is the sequence of radioactive decays in a sample of an unstable nuclear isotope. For example, consider a sample of Ra-226. This isotope is unstable—it undergoes alpha-decay—and the half-life is (Subscript[t, 1/2]) = 1600 years.

The probability per unit time for any single nucleus to decay, which is called the "decay rate", is   Γ = ln(2) / (Subscript[t, 1/2])  = xxx (sec^(-1)).
So if we use one second as a standard time interval, δt = 1 sec, then the probability of decay during one second is p = xxx << 1. Since this is very small, the process of radioactive decays is a Poisson process.

Note that only the mean μ appears in P(r;t) so that knowledge of n and p is not always necessary. The mean in refers to the mean number in units of the time.

I.e., µ = Γ Δt, and we can rewrite the Poisson distribution as

P( n ; Δt ) = (Γ ((Δt))^n) exp(–Γ Δt) /n!

EXERCISE 2
P(n) = (μ^(n )) exp(–μ) / n! 
(A) Calculate the theoretical mean of n.
(B) Also calculate the theoretical variance.

EXERCISE 3
Plot the Poisson distribution P(n) for μ = 5, 10, 40.
Put all three plots on one graph.
Remember, the domain of the Poisson distribution is Z ( = the non-negative integers.) So use DiscretePlot.

We see from Exercise 3 that as µ becomes large, the distribution becomes more and more symmetric and approaches a Gaussian form. We are replacing a discrete distribution by a continuous one. This brings us to the third important probability distribution—the Gaussian distribution.

2.3 The Gaussian or Normal Distribution

The Gaussian or normal distribution is the most familiar distribution in all the sciences. Measurement errors, 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 Gaussian is a continuous, symmetric distribution whose density is given by

EXERCISE 4

(A) Calculate the mean and variance for the Gaussian distribution with parameters μ and σ.

(B) Plot the Gaussian distribution for μ = 0, 1, 2 and σ = 1, 2, 3.

(C) The full width at half maximum (FWHM) for the Gaussian distribution is
FWHM = A σ, where A is a numerical constant. Determine A.

The areas under the Gaussian between integral intervals of sigma are important in the publication and interpretation of measurement errors. For example, if an experimental result is presented as x = a ± b, it means that the probability that the true results lies between a–b and a+b is 68%.

EXERCISE 5

(A) What is the area under the Gaussian distribution, for x ∈ ( μ–σ , μ+σ ).
This is the probability for x to lie within 1-sigma of the mean.

(B) Repeat the calculation for the 2-sigma range.

(C) Repeat the calculation for the 3-sigma range.

In high-energy physics, a fundamental discovery must pass the 5-sigma test.

Radium

All isotopes of radium have half-lives much shorter than the age of the Earth, so that any primordial radium would have decayed long ago. Radium nevertheless still occurs in the environment, because the isotopes 223 Ra, 224 Ra, 226 Ra, and 228 Ra are part of the decay chains of natural thorium and uranium isotopes; since thorium and uranium have very long half-lives, the radium “daughters” are continually being regenerated by their decay.[9] Of these four isotopes, the longest - lived is 226 Ra (half - life 1600 years), a decay product of natural uranium.Because of its relative longevity, 226 Ra is the most common isotope of the element, making up about one part per trillion of the Earth' s crust; essentially all natural radium is 226 Ra.[21] Thus, radium is found in tiny quantities in the uranium ore uraninite and various other uranium minerals, and in even tinier quantities in thorium minerals.One ton of pitchblende typically yields about one seventh of a gram of radium.[22] One kilogram of the Earth' s crust contains about 900 picograms of radium, and one liter of sea water contains about 89 femtograms of radium.[23]

The common historical unit for radioactivity, the curie, is based on the radioactivity of 226 Ra.[32]

Some of the biological effects of radium were apparent from the start. The first case of so-called “radium-dermatitis” was reported in 1900, only 2 years after the element’s discovery. The French physicist Antoine Becquerel carried a small ampoule of radium in his waistcoat pocket for 6 hours and reported that his skin became ulcerated. Pierre and Marie Curie were so intrigued by radiation that they sacrificed their own health to learn more about it. Pierre Curie attached a tube filled with radium to his arm for ten hours, which resulted in the appearance of a skin lesion, suggesting the use of radium to attack cancerous tissue as it had attacked healthy tissue.[71] Handling of radium has been blamed for Marie Curie’s death due to aplastic anemia. A significant amount of radium’s danger comes from its daughter radon: being a gas, it can enter the body far more readily than can its parent radium.[35]

Today, 226Ra is considered to be the most toxic of the quantity radioelements, and it must be handled in tight glove boxes with significant airstream circulation that is then treated to avoid escape of its daughter 222Rn to the environment. Old ampoules containing radium solutions must be opened with care because radiolytic decomposition of water can produce an overpressure of hydrogen and oxygen gas.[20]

●

The most common isotope of radium is 88(Ra^226) .
It is unstable and decays to a radon isotope by alpha-decay,
	(Ra^(226 ))⟶ (Rn^222) + α ; half life = 1600 years.
Radon is also unstable, so the decay of radium ultimately produces ionizing radiation,
which can excite fluorescent materials ( radio-luminescence ). In the early 20th century, watch dials were painted with radium to make them luminous. However, radium is toxic, so today radium has no practical applications outside of nuclear medicine.

In nature, radium is found in uranium ores in trace amounts—about 1/7 gram per ton.

EXERCISE 8
Consider a sample of pure Ra-226, amounting to 1.0(x10^(-6)) gram.
Let Γ = the decay rate for Ra ⟶ Rn + α.
The probability that a radium atom will decay during the infinitesimal time dt is Γ dt.
(A) From the half-life (1600 years) calculate Γ , and express Γ in (((sec))^(-1)) .
(B) Calculate the mean number of decays in 1 hour, and call this number A.
(C) Calculate the probability that the number of decays in 1 hour will be less than 0.25A.
(D) Calculate the probability that the number of decays in 1 hour will be more than 0.75A.
(E) Plot a graph that shows the probability distribution for the value of A.