CPS713 APPLICATIONS OF COMPUTATIONAL SCIENCE: STATISTICS FOR PHYSICS

Geoffrey Fox

Syracuse University

Northeast Parallel Architectures Center

111 College Place Syracuse, New York 13244-4100

(315) 443-2163 gcf@npac.syr.edu

http://www.npac.syr.edu

February, 1996

Statistics for Physics

Abstract of Statistics for Physics

• These notes were developed in a set of courses taught at Caltech, 1980–1989.

• They describe statistics with a particular emphasis on parameter estimation, as occurs in many physics experimental analysis problems.

• Although the notes were developed for physics majors, they can be read by any computational scientist.

• We cover the following topics:

• basic definitions and results,

• important discrete distributions,

• estimation of parameters,

• minimization,

• goodness of fit,

• and generation of random numbers.

Useful Books for Statistics Part of Course

• G. P. Yost, “Lectures on Probability and Statistics”

— copies available from G. Fox.

• W. T. Eadie, et al. (North Holland), “Statistical Methods in Experimental Physics” — a little sophisticated for this course, but technically excellent (called EDJRS).

• W. H. Press, et al., “Numerical Recipes — The Art of Scientific Computing”, Chapter 7, 13, 14 — Computer Related Issues. This has C, FORTRAN, and PASCAL programs available.

• J. Mathews and R. L. Walker (Benjamin), “Mathematical Methods of Physics”, Chapter 14.

• Cramer: Mathematical Methods of Statistics, published by Princeton.

Statistics for Physics

Table of Contents

1 Basic Definitions and Results  5
1.1 Basic Definitions and Results  6
1.2 Bayes Theorem  11
1.3 Continuous Distributions  20
1.4 Functions of a Random Variable  23
1.5 Properties of Joint Distribution Functions  27
1.6 Means, Moments, Correlations  29
1.7 Central Limit Theorem  43
2 Feefoefum in Las Vegas:  62

Combinatoric and Discrete Probability

2.1 Binomial and Poisson Distributions  63
2.2 Birth and Death Processes  78
3 Estimation of Parameters  89
3.1 Maximum Likelihood Principle  90
3.2 χ^{2} Method  121
3.3 Counter Physics  132
3.4 Bubble Chamber  139

4 Minimization  157
4.1 Statement of Problem  158
4.2 Derivative Formalism for χ^{2}  159
4.3 The Derivative Method for χ^{2}  165
4.4 The Real World and χ^{2} Minimization 168
5 Goodness of Fit  180
5.1 The χ^{2} Distribution  181
5.2 Confidence Level  185
5.3 Bartlett’s S Function  189
5.4 Hypothesis Testing  191
5.5 Student’s t Distribution  194
5.6 Robust Estimation  197
6 Generation of Random Numbers 211

Statistics for Physics  1: Basic Definitions and Results

1: BASIC DEFINITIONS AND RESULTS 1.1 Basic Definitions and Results

Is Probability Intuitive

or Mathematical?

• Formal Basis

• A major difficulty in formulating the theory of probability is that one must relate it to the real world.

• Cramer amusingly notes that some of the axiomatic approaches to probability are comparable to definitions of geometry in terms of the fundamental idea of a point as a limit of a spot of chalk dust of smaller and smaller size.

• I will first dispose of the 24 lectures that I suffered on formal probability by noting that the theory can be based on the following axioms:

Basic Definitions—Borel Sets

• Consider a variable point X in space Y where Y is some subset of ℜ_{k} (ℜ_{k} is the k dimensional real Euclidean space).

• Further, consider the family of all Borel sets S in Y.

• You will, of course, recall that a Borel set is gotten by performing (finite or denumerably many times) on intervals [a, b], (a, b), (a, b], [a, b) … the operations of addition, subtraction or multiplication (and continuing process denumerably often).

Axioms of Probability

• Axiom 1: To every set S corresponds a non-negative number P (S) which is called the probability of the event X belonging to S.

• Axiom 2: P (Y) = 1

• Axiom 3: P (S) is a completely additive set function: i.e., P (S_{1} + S_{2} + …) = P (S_{1}) + P (S_{2}) … + … if the S_{i} are disjoint. S_{i}S_{j} = 0, i ≠ j

Immediate Deductions from

Axioms of Probability

• Note that it follows from the axioms that

• (a) 0 ≤ P (S) ≤ 1: any S

Proof: P (S) ≥ 0 was postulated. Take S_{1} = S, S_{2} = Y − S = S^{∗} (the complement of S) in Axiom 3. P (S_{1}) + P (S_{2}) = P (Y) = 1, but P (S_{2}) ≥ 0 ⇒ P (S_{1}) = P (S) ≤ 1.

• (b) P (S + 1 + S_{2}) ≤ P (S_{1}) + P (S_{2}): any S_{1}, S_{2} with equality if S_{1} and S_{2} are disjoint. Proof:

S_{1} + S_{2} = (S_{1} − S_{1} S_{2}) + S_{2}
P (S_{1} + S_{2}) = P (S_{1} − S_{1} S_{2}) + P (S_{2})

are sets disjoint. But P (S_{1} − S_{1} S_{2}) ≤ P (S_{1}) — see (a). This proves desired result.

From Mathematical to Intuitive

Definition of Probability

• Anyhow, the variable point X is now called a random variable and P (S) is the probability function of X.

• To relate this definition to the outside world, we must resort to the vanishing chalk dust approach as espoused by Mathews and Walker.

• If we flip an (unbiased) coin N times, we expect

lim_{N→∞} \frac{N_{h}}{N} = \frac{1}{2}

where N_{h} is the number of times it turns up heads.

• We define

P (heads) = lim_{N→∞} \frac{N_{h}}{N} = \frac{1}{2}

• This is a random variable with only two possible values — heads or tails. Alternately, we could say it takes values x = 0 and 1.

• Then, in our previous language any Borel set containing 0 and not 1 has P (S) = 1/2, any containing 0 and 1, P (S) = 1, etc.

• This illustrates the relation of intuitive frequency definition of probability and that based on mathematical measure.

1.2 Bayes Theorem: Conditional

Probabilities—1) Frequency Approach

• Take two possible (Borel) subsets of Y ∈ ℜ_{k}. Call them A and B. A could be point x = 0. B could be point x = 1 only.

• Generally, Y is divided into four subspaces by A and B, and following the chalk dust method—imagine N experiments, which give the following number of occurrences N_{i} for events to fall in the four subspaces.

N_{1} = n(A − AB)
N_{2} = n(B − AB)
N_{3} = n(AB)
N_{4} = n ((A + B)^{∗} complement of (A + B))

where n(S) is number of events falling in subspace, S.

CPS713 STAT G.C.Fox, Feb 96  11  gcf@npac.syr.edu

Statistics for Physics  1: Basic Definitions and Results

Bayes Theorem: Conditional Probabilities—2) Frequency and

\underline{Probability}

• Then N_{1} + N_{2} + N_{3} + N_{4} = N.

• We previously defined P (S) and now we expect

P (A) = lim_{N→∞} \frac{N_{1} + N_{3}}{N}
P (B) = lim_{N→∞} \frac{N_{2} + N_{3}}{N}
P (A + B) = lim_{N→∞} \frac{N_{1} + N_{2} + N_{3}}{N}  (1)

• This is the usual relationship of probabilities to frequencies.

CPS713 STAT G.C.Fox, Feb 96  12  gcf@npac.syr.edu

Bayes Theorem: Conditional

Probabilities—3) Bayes Law

• Now we define conditional probability P (S/A) (read as probability of S given A has occurred) as P (SA)/P (A).

• Note this is a probability according to axioms, if we restrict “bounding space” Y to be A.

• Then, in terms of frequencies,

P (A/B) = \frac{N_{3}}{N_{2} + N_{3}} as N → ∞
P (B/A) = \frac{N_{3}}{N_{1} + N_{3}} as N → ∞  (2)

whence

P (B)P (A/B) = P (A)P (B/A)  (3)

which is Bayes Law.

• Of course Equation (3) is obvious from definition: we didn’t need frequencies to prove it.

Estimation of Image

Scanning Efficiency—I

• Consider (example taken from p. 11 EDJRS) the scanning of a set of pictures from either a bubble chamber or the swimsuit issue of Sports Illustrated. Let

• n_{tot} — pictures with some property of interest • n(A) — number found by an undergraduate

during finals week

• n(B) — number found by a graduate student working on a theory in 26 dimensions

• Here, A: Scan Number One B: Scan Number Two

• We can tag members of C = AB (A and B).

Example: Estimation of Image

Scanning Efficiency—II

• Then, P (C) = P (A)P (B) if we assume scans A, B are independent.

• (in practice not true as say short tracks or larger swimsuits would be systematically missed by either graduate or undergraduates!).

• Therefore, making independence assumption

\frac{n(C)}{n_{tot}} = \frac{n(A)}{n_{tot}} · \frac{n(B)}{n_{tot}}

or

n_{tot} = \frac{n(A)n(B)}{n(C)}

is desired estimate of total number of events.

• In practice, second scan will usually only be done on a subset of data sample.

• Note if A, B independent

P (A/B) = \frac{P (AB)}{P (B)} = P (A)

independent of B.

CPS713 STAT G.C.Fox, Feb 96  15  gcf@npac.syr.edu

Statistics for Physics  1: Basic Definitions and Results

Example: Analysis of Unreliable Test

• Consider (page 18, Yost) set B to categorize health of people and set A to be results of unreliable test.

• B has two possible results.

B_{1} = 0 ill
B_{2} = 1 well

• Consider the situation with A = “Test Positive”

• Bayes gives for B = B_{1} or B_{2}

P (B/A) = \frac{P (B)P (A/B)}{P (A)}

• Apply separately B = B_{1}, or B_{2}, i.e., B = B_{1} + B_{2} and P (B/A) = P (B_{1}/A) + P (B_{2}/A).

• Then given data proposed by Yost

P (B_{2}) P (A/B_{2}) P (B_{1}) P (A/B_{1})

P (B/A) = \frac{[.999 × 0.02δ(1 − x) + .001 × .99δ(x)]}{[P (A) = .001 × .99 + .02 × .999]}
P (B_{1}/A) = “ill” = .047
P (B_{2}/A) = “well” = .953
CPS713 STAT G.C.Fox, Feb 96  16  gcf@npac.syr.edu

The Bayes Controversy

Naive Application for Estimation—I

• Now we have derived Bayes theorem where A and B are (Borel) sets in a space populated by a random variable.

• Now a useful “application” of theorem is take the following choice of sets:

A: to be Borel set given a set of experimental

results—this is correct. B: to be a set of hypotheses, e.g., a set in space of

values of a theoretical parameter—this will turn out to be flaky.

Consider an example:

• We measure scattering angles θ and wish to find α

• when probability of observing θ is P (θ) = \frac{1}{2(1+α cos θ)}.

• Here, A is a set in θ space.

• B is a set in α space.

The Bayes Controversy

Naive Application for Estimation—II

• Then experimentalists measure P (A/B) i.e., values of θ subject to known (by the Almighty and those on fourth floor of Lauritsen—the theoretical physics group at Caltech) single value of α. Bayes says remarkably, that:

P (B/A) = P (A/B)\frac{P (B)}{P (A)}

• This formula for P (B/A) gives you what you really want—the “probability distribution of α” consequent from measurements of θ.

The Bayes Controversy

Naive Application for Estimation—III

• Problems with this application of Bayes Theorem are that:

• a) α is not a random variable—it only takes one value. Answer: Our knowledge of α can be considered as a random variable.

• b) What is P (B)? We will discuss this so called a priori knowledge of B later on—it essentially summarizes all previous (other) knowledge of α gained from previous (other) experiments.

• c) What is P (A)? The probability of measurement of θ given any possible α.

• Nobody knows what this is!

• You essentially determine it by normalizing P (B/A) or equivalently note that relative probabilities

\frac{P (B_{1}/A)}{P (B_{2}/A)} = \frac{P (A/B_{1}) P (B_{1})}{P (A/B_{2}) P (B_{2})}

are indeed independent of the noxious P (A).

CPS713 STAT G.C.Fox, Feb 96  19  gcf@npac.syr.edu

1.3 Continuous Distributions

General (Borel) Formalism

• To be formal and reveal relics of a primeval mathematical training: take random variables in ℜ_{1} and consider special Borel set ξ ≤ x: It has probability F (x) with F (∞) = 1.

• Now, in some cases, we can differentiate F (x) and so write:

F (x) = x ; ∫ ; −∞ p(x)dx

where p(x) is called the probability density.

• One can clearly interpret p(x)dx as the probability that random variable x takes values in range x to x + dx.

• Consider examples:

• (a) In the case of tossed coin, let heads be x = 1 and tails x = 0. Then for (an unbiased) coin

p(x) = \frac{1}{2} [δ(x) + δ(x − 1)]

Gaussian Distribution—I

• (b) Another important example is

p(x) = \frac{1}{\sqrt{2π}σ} exp [−\frac{1}{2}\frac{(x − μ)^{2}}{σ^{2}}]  (4)

• This is the probability density of the normal or Gaussian distribution—so called because it was first discussed by de Moivre in 1733.

Statistics for Physics  1: Basic Definitions and Results

Properties of Gaussian Distribution

• Note that we have ensured in Equation (4) that

+∞ ; ∫ ; −∞p(x)dx = 1

• Also, the associated distribution function is:

F (x) = Φ(x) = \frac{1}{\sqrt{2π}σ} x ; ∫ ; −∞dy exp [−\frac{(y − μ)^{2}}{2σ^{2}}]

• which cannot be done analytically except by relating it to the error integral:

erf(x) = \frac{2}{\sqrt{π}} x ; ∫ ; 0exp (−y^{2}) dy

Statistics for Physics  1: Basic Definitions and Results

1.4 Functions of a Random Variable

• If x is a random variable, then so is y = f(x).

• Let x have probability distribution p(x) and y probability distribution q(y).

• Then the probability of x lying in [x, x + dx] is p(x)dx:

• Let [y, y + dy] be associated interval in y.

• Clearly, p(x)dx = q(y)dy and

• \frac{dy}{dx} = |f′(x)|,

• so

p(x) = q (f(x)) |f′(x)|  (5)

• There are some obvious modifications of this if the x ↔ y relation is not single-valued.

Joint Probability Distributions

• We can now consider the important extension to more than one variable.

• Suppose we have two random variables x and y each taking values in ℜ_{1}.

• Then we can conceive of a random variable z which is the ordered pair (x, y).

• Even if we know the probability distributions of x and y

x → p(x)dx
y → q(y)dy

• The distribution of z is of course not defined.

• Suppose that z has density r(x, y) so that probability that z lies in volume element dx dy is r(x, y)dxdy.

• A priori, all we know is that r satisfies

∫ r(x, y)dy = p(x)
∫ r(x, y)dx = q(y)

• These are called marginal distributions (EDJRS p. 18).

Statistics for Physics  1: Basic Definitions and Results

Example of Joint Probability

Distributions: Bayes Law

for Densities—I

• Let event X be x lies in [x, x+dx]

• Let event Y be y lies in [y, y+dy]

• Let T (X/Y) = t(x|y)dx be probability of X given Y.

• By definition

T (X/Y) = \frac{P r(XY)}{P r(Y)} = \frac{r(x, y)}{q(y)} \frac{dx dy}{dy}

• or t(x|y) = r(x, y)/q(y) a conditional distribution (EDJRS p. 18).

• Similarly, let S(Y/X) = s(y|x)dy be probability of Y given X.

• We find

s(y|x) = \frac{r(x, y)}{p(x)}

CPS713 STAT G.C.Fox, Feb 96  25  gcf@npac.syr.edu

Example of Joint Probability

Distributions: Bayes Law

for Densities—II

• or eliminating r(x, y) from equations for t(x|y) and s(y|x), we find

t(x|y)q(y) = s(y|x)p(x)  (6)

• This is a density version of Bayes law.

• Example of a function of a random variable:

• Let y = f(x). Now y is a known function of x, and clearly s(y|x) = δ(y − f(x)).

• Then using previous equation for s(y|x), we can in this case find

r(x, y) = δ [y − f(x)] p(x).

1.5 Properties of Joint

Distribution Functions

• (i) Independence:

• If the distribution of z satisfies

r(x, y) = p(x)q(y)  (7)

• then we say that x and y are statistically independent.

• (ii) Function of Two Variables

• Suppose, rather than z—an ordered pair in ℜ_{2}, we consider h(x, y) in ℜ_{1}: a conventional function of x and y.

• Then w = h(x, y) is also a random variable with distribution d(w) where

d(w) = ∫∫ dx dy δ [w − h(x, y)] r(x, y) (8)

• This may perhaps be proved very prettily: it is obvious if you just do one—say x-integral in Equation (8), and use the same argument as already given for functions of a random variable.

Statistics for Physics  1: Basic Definitions and Results

Means—Moments, etc.

for One-dimensional Distributions

• Given a one-dimensional distribution, we can define:

• mean: ⟨x⟩ = +∞ ; ∫ ; −∞xp(x)dx

• n’th moment: ⟨x^{n}⟩ = ∞ ; ∫ ; −∞x^{n}p(x) dx

• expectation: ⟨f(x)⟩ = ∞ ; ∫ ; −∞f(x)p(x)dx which “follows” from n’th moment by Taylor expansion of f(x). This can be considered as the mean of the random variable z = f(x).

• n’th central moment:

⟨(x − ⟨x⟩)^{n}⟩ = +∞ ; ∫ ; −∞(x − ⟨x⟩)^{n} p(x)dx

• Variance: V_{x} = ⟨(x − ⟨x⟩)^{2}⟩ = 2nd central moment • Standard deviation σ_{x}: V_{x} = σ_{x}^{2} where σ_{x} is used much more than V_{x}.

• Note:

V_{x} = ⟨x^{2}⟩ − 2⟨x⟩⟨x⟩ + ⟨x⟩^{2} = ⟨x^{2}⟩ − ⟨x⟩^{2} (9)

a result which is used continuously.

CPS713 STAT G.C.Fox, Feb 96  28  gcf@npac.syr.edu

1.6 Means–Moments–Correlations

for Multidimensional Distributions I

• Given n-dimensional distribution p(x_{1}, … x_{n}):

• The mean now becomes an n-dimensional vector

⟨x_{k}⟩ = ∫ ⊗ dx_{i} x_{k} p(x_{1} … x_{n})

• The variance becomes an n × n moment matrix M

M_{kl} = ∫ ⊗ dx_{i} [x_{k} − ⟨x_{k}⟩] [x_{l} − ⟨x_{l}⟩] p (x_{1} … x_{n})

• Note for we can show—analogously to Equation (9)—that

M_{kl} = ⟨x_{k}x_{l}⟩ − ⟨x_{k}⟩ ⟨x_{l}⟩  (10)

• Also note that M is symmetric and positive semi-definite. [This follows because y_{k} M_{kl} y_{l}—any y—clearly has a ≥ 0 integrand].

• Now the diagonal terms of M are called variances which are ≥ 0:

• Their square roots are standard deviations again.

• Further, when k ≠ l, we put ρ_{kl} = M_{kl}/\sqrt{M_{kk} M_{ll}}, which is ≤ 1 and is called the correlation coefficient.

CPS713 STAT G.C.Fox, Feb 96  29  gcf@npac.syr.edu

Means—Moments—Correlations

for Multidimensional Distributions II

• Note that if ρ_{kl} = 0—all k ≠ l—then we say that x_{1} … x_{n} are uncorrelated.

• Further, if x_{1} … x_{n} are independent, then Equation (10) shows at once that they are uncorrelated.

• The converse does not follow of course (except under special circumstances, e.g., if x_{i} Gaussianly distributed).

• One can define higher order tensors M_{klmn...} but these are not too important as in usual applications x_{1} … x_{n} are n functions of N >> n observations.

• Then central limit theorem—to which we will come—implies that only the mean and M_{kl} are important asymptotically (in N).

Errors and Moment Matrices

• The Moment Matrix is used to quantify errors in experimental measurements.

• Often, the results of an experiments on x_{1} … x_{n} are quoted in terms of means ± errors:

• These errors are standard deviations = \sqrt{diagonal elements of moment matrix}.

• Now, this is throwing away information—the off diagonal terms in M—which can be important.

• Thus, suppose some years after the experiment measures and reports \underline{x}, a theorist decides that particular linear combinations \underline{y} = T \underline{x} (T an n × n matrix) are useful.

• What are errors on y?

⟨y_{k}y_{l}⟩ = T_{kk′} ⟨x_{k′}x_{l′}⟩ T_{l′l}^{T}  or
M_{y} = T M_{x}T^{T}  (11)

• To find diagonal terms in M_{y}—i.e., the errors^{2} on our new theorist’s y—requires off diagonal terms in M_{x}.

Density Matrix Element Example—I

• An example of importance of off-diagonal error matrix elements can be found from high energy physics in the decay of a spin 1 particle.

• One can define the angular distribution by density matrix elements

x_{1} = ρ_{00}
x_{2} = ρ_{11}
x_{3} = ρ_{10}

in some kinematic coordinate system (Lorentz frame)!

• Some years ago, these were always given in one frame—the so-called t channel frame.

• “Recently” (1975—these really are old notes!), it was found that the values in another frame (the “s channel”) were more meaningful physically.

• Unfortunately, latter values, y, were related by precisely such a linear transformation y = T x and one can no longer estimate errors in y from recorded experimental data.

Density Matrix Element Example—II

• The x_{i} and y_{i} in this example correspond to expressing angular decay

W (θ, φ) = 1 + 3 ; ∑ ; i=1 x_{i}X_{i}(θ, φ),  or
W (θ, φ) = 1 + 3 ; ∑ ; i=1 y_{i}Y_{i}(θ, φ)

• where y_{i} are linear combinations of the x_{i} corresponding to the X_{i}(θ, φ) being linear combinations of the Y_{i}(θ, φ).

• The details are different, but you can see the essential points from working through the choice

X_{1} = cos θ Y_{1} = P_{1}(cos θ)
X_{2} = cos^{2} θ Y_{2} = P_{2}(cos θ)
X_{3} = cos^{3} θ Y_{3} = P_{3}(cos θ)

• where P_{l} are Legendre Polynomials.

• This calculation is left as an exercise for the reader.

Statistics for Physics  1: Basic Definitions and Results

Examples of Correlations 2D Track Finding—I

1. Measurement of Tracks in Chambers

• Track is x = a + bz (a =intercept, b = slope), where

b = \frac{(x_{2} − x_{1})}{(z_{2} − z_{1})}
a = \frac{z_{2}x_{1} − z_{1}x_{2}}{z_{2} − z_{1}}

• x_{i} are measured as independent random variables, but the combinations a and b are clearly correlated.

⟨(b − ⟨b⟩)(a − ⟨a⟩)⟩ = −\frac{z_{2}σ_{1}^{2} − z_{1}σ_{2}^{2}}{(z_{1} − z_{2})^{2}}

• where x_{i} has standard deviation σ_{i}.

CPS713 STAT G.C.Fox, Feb 96  34  gcf@npac.syr.edu

Examples of Correlations

2D Track Finding—II

• (Note: If σ_{1} = σ_{2} and z_{1} = −z_{2}, i.e., origin z midpoint between chambers, then a and b are uncorrelated.)

• However, the origin, z = 0 (in our experiment E110 at Fermilab), was a convenient surveyor’s mark and so was not chosen to minimize correlation.

• In the case of many measurement chambers, the slope and intercept will still be correlated but there will be no simple way to get rid of it.

• In practice, it is important to keep full correlation matrix.

• For instance, estimation of best vertex from intersection of several tracks needs the full matrix.

Examples of Correlations

Exponential Fits—I

• Consider, as another example, our Measurement of Regge Pole Intercept

• In a Fermilab experiment, E350, of times gone by, we measured

\frac{dσ}{dt dx} = g(t)(1 − x)^{1−2α(t)}

where x ≤ 1(.5 ≤ x ≤ 1 in the experiment).

• g(t) is an unknown but fixed (at fixed t) constant of little interest.

• α(t) is extremely interesting.

Examples of Correlations

Exponential Fits—II

• In a (non-linear) fit, g(t) and α(t) are strongly correlated

• so that 1 ; ∫ ; .5g(t)(1 − x)^{1−2α(t)} is approximately

preserved.

• This integral sums all measurements, and is very well determined experimentally, and so defines best determined mix of theoretical parameters.

• Unfortunately, we only want α(t) (predicted by theory) and can not use fact that a particular linear combination of α(t) and g(t) is predicted much better than either parameter on its own.

Errors in Linear Combinations

of Uncorrelated Variables

• In practice, when off diagonal terms are not given, one assumes they are zero, i.e., variables are uncorrelated.

• This assumption/assertion gives:

σ_{y_{i}}^{2} = n ; ∑ ; j=1 (T_{ij})^{2} σ_{x_{j}}^{2}  (12)

which is the usual rule of addition of errors in quadrature.

• One has to emphasize that this “usual” rule is only valid where x_{j} uncorrelated, and usually x_{j} are correlated!

Statistics for Physics  1: Basic Definitions and Results

Error in Lack of

Correlation Assumption

• If y = a_{1}x_{1} + a_{2}x_{2},

• then σ_{y}^{2} = a_{1}^{2}σ_{1}^{2} + a_{2}^{2}σ_{2}^{2} +2a_{1}a_{2}σ_{1}σ_{2}ρ_{12},

• and all one can say in general is that σ_{y}^{2} ≤ (|a_{1}σ_{1}| + |a_{2}σ_{2}|)^{2}

• The uncorrelated estimate is σ̃_{y} = \sqrt{a_{1}^{2}σ_{1}^{2} + a_{2}^{2}σ_{2}^{2}}

• and the true σ_{y} satisfies 0 ≤ σ_{y} ≤ \sqrt{2}σ̃_{y}

• where extremes are only attained when |a_{1}σ_{1}| = |a_{2}σ_{2}| which is a most dangerous case.

• Other ratios, |a_{1}σ_{1}|/|a_{2}σ_{2}| ≠ 1 are better, e.g., a_{1}σ_{1} = 0 ⇒ σ_{y} = σ̃_{y}!

• For n variables, all one can say rigorously is 0 ≤ σ_{y} ≤ \sqrt{n}σ̃_{y}

which is clearly pretty bad for large n.

CPS713 STAT G.C.Fox, Feb 96  39  gcf@npac.syr.edu

Statistics for Physics  1: Basic Definitions and Results

Parameters of a Gaussian Distribution—I

• Suppose x is Gaussian, i.e.,

p(x) = \frac{1}{\sqrt{2π}σ} exp \{−\frac{[x − μ]^{2}}{2σ^{2}}\}

• Then the mean is given by

⟨x⟩ = +∞ ; ∫ ; −∞x p(x)dx
= μ + +∞ ; ∫ ; −∞(x − μ)p(x)dx
= μ

at last integral clearly vanishes from antisymmetry.

• So, fundamental parameter, μ, in Gaussian is just the mean.

Parameters of a

Gaussian Distribution—II

• Variance

= \frac{1}{\sqrt{2π}σ} +∞ ; ∫ ; −∞(x − μ)^{2}p(x)dx
= \frac{σ^{3}}{\sqrt{2π}σ} +∞ ; ∫ ; −∞dy \frac{y^{2}}{σ^{3}} exp [\frac{−y^{2}}{2σ^{2}}]

• Putting y = x − μ and multiplying, and dividing by σ^{3}.

V ariance = \frac{σ^{3}}{\sqrt{2π}σ} \frac{∂}{∂σ} \{ +∞ ; ∫ ; −∞e^{−y^{2}/2σ^{2}}dy\}
= \frac{σ^{2}}{\sqrt{2π}} \frac{∂}{∂σ} [\sqrt{2π}σ]
= σ^{2}

• so fundamental Gaussian parameter σ is standard deviation—as its nomenclature suggested.

Addition of Independent

Random Variables

• Suppose y = \frac{1}{n} n ; ∑ ; i=1x_{i}

• where the x_{i} are independent and they all have the same distribution with mean μ and standard deviation σ.

• Then clearly ⟨y⟩ is also μ.

• Further, to find standard deviation, calculate

⟨(y − μ)^{2}⟩ = \frac{1}{n^{2}} _{i,j}^{∑} ⟨(x_{i} − μ) (x_{j} − μ)⟩
= \frac{σ^{2}}{n}

as all off diagonal terms vanish.

• So it is precisely true that y has mean μ and standard deviation \frac{σ}{\sqrt{n}}.

• Later, we will show that as n → ∞ (i.e., not exactly), all other moments of y vanish faster than 1/n

• and so y becomes Gaussian whatever the original distributions of x.

1.7 Leading Up to the Central

Limit Theorem Moment

Generating Function—I

• The Moment Generating Function which is also called characteristic function, is defined in terms of probability distribution p(x) by

φ(k) = +∞ ; ∫ ; −∞dx e^{ikx}p(x)
= 1 + ik⟨x⟩ − \frac{k^{2}}{2!} ⟨x^{2}⟩ + …

• Note the Fourier inversion theorem gives

p(x) = \frac{1}{2π} +∞ ; ∫ ; −∞dx e^{−ikx}φ(k)

Moment Generating Function of a

Sum of Two Random Variables

• In an important application of moment generating functions, put h = x + y where x and y are independent.

• From joint probability distribution, it follows that the distribution density of h is

w(h) = +∞ ; ∫ ; −∞ +∞ ; ∫ ; −∞δ(h − x − y) p(x) q(y) dx dy

φ_{h}(k) = +∞ ; ∫ ; −∞dh w(h)e^{ikh}
= +∞ ; ∫ ; −∞ +∞ ; ∫ ; −∞e^{ik(x+y)} p(x) q(y) dx dy
= φ_{x}(k)φ_{y}(k)

• So characteristic functions multiply for addition of two independent random variables.

• This result can obviously be extended to sum of n independent random variables.

Moment Generating Function III

• Expanding the exponential in definition of φ(k) give

φ(k) = _{n}^{∑} k^{n}E\frac{(ix)^{n}}{n!}

• where E is the expectation or mean value.

φ(k) = 1 + i⟨x⟩k − \frac{⟨x^{2}⟩k^{2}}{2}, etc.

• For a Gaussian distribution

φ(k) = \frac{1}{\sqrt{2π}σ} +∞ ; ∫ ; −∞exp(ixk) exp [−\frac{1}{2}\frac{(x − μ)^{2}}{σ^{2}}]
= exp \{iμk − \frac{1}{2}k^{2}σ^{2}\}

• EJDRS define the Cumulant as

K(k) = ln φ(k)

• I have never used this! Note that

K(k) = ⟨x⟩ik − \underbrace{(⟨x^{2}⟩ − ⟨x⟩^{2})}_{σ^{2}} \frac{k^{2}}{2}

• i.e., as x → x + μ only shifts K by additive constant ikμ, we only find “central” moments for coefficients of k^{2} and above.

CPS713 STAT G.C.Fox, Feb 96  45  gcf@npac.syr.edu

Central Limit Theorem Statement

\underline{Laplace}

• Let y = \frac{1}{n} n ; ∑ ; i=1x_{i}: x_{i} independent

• (a) If x_{i} are Gaussian, then so is y: in fact, if y Gaussian exactly, then so must all the x_{i} be Gaussian.

• (b) De Moivre in 1733 showed that if x_{i} binomially distributed, then y was always Gaussian for large x (even though Gauss didn’t exist in 1733).

• (c) Laplace, in 1812, showed that de Moivre’s result was generally true for all distributions of x_{i}.

• Central Limit Theorem of Laplace—Cramer p. 215 If x_{1} … x_{n} are independent random variables with the same distribution, then y is asymptotically Gaussian with mean μ and standard deviation σ/\sqrt{n} (where x_{i} has mean μ and standard deviation σ, and both μ and σ exist, i.e., are finite).

Statistics for Physics  1: Basic Definitions and Results

Central Limit Theorem

Statement—Liapounoff

• This extended theorem reads:

• If x_{1} … x_{n} are independent random variables with possibly different distributions, and:

μ_{i} = ⟨x_{i}⟩
σ_{i}^{2} = ⟨(x_{i} − μ_{i})^{2}⟩
ρ_{i}^{3} = ⟨|x_{i} − μ_{i}|^{3}⟩

• where ρ_{i} must exist (not necessary in Laplace’s version where all x_{i} had same distribution).

• Let

ρ^{3} = n ; ∑ ; i=1 ρ_{i}^{3}  and
Σ^{2} = n ; ∑ ; i=1 σ_{i}^{2}

• Then if lim_{n→∞} ρ/Σ → 0,

• then y is Gaussianly distributed with • mean μ = ∑_{i=1}^{n} μ_{i}/n and

• standard deviation Σ/n.

• Note the lim_{n→∞} ρ/Σ does → 0 when ρ exists in Laplace’s case as ρ ∝ n^{1/3} and ∑ ∝ n^{1/2}.

CPS713 STAT G.C.Fox, Feb 96  47  gcf@npac.syr.edu

Proof of Central Limit Theorem

\underline{Laplace—I}

• Now we will prove the first (simple) version of the central limit theorem.

• From our discussion of characteristic functions

φ_{y}(k) = [φ_{x} (\frac{k}{n})]^{n}

where φ_{y} and φ_{x} are characteristic functions of y and x respectively.

• Now,

φ_{x}(k) = exp(ikμ) ∫ dxe^{ik(x−μ)}p(x)
= exp(ikμ) [1 − \frac{k^{2}σ^{2}}{2} + O (k^{3})]

Proof of Central Limit Theorem

\underline{Laplace—II}

• So, the characteristic function of y is:

φ_{y}(k) = exp(ikμ) [1 − \frac{k^{2}σ^{2}}{2n^{2}} + O (\frac{k}{n})^{3}]^{n}
= e^{ikμ} \{exp [−\frac{k^{2}σ^{2}}{2n^{2}} + O (\frac{k}{n})^{3}]\}^{n}
= exp(ikμ) exp \{−\frac{k^{2}σ^{2}}{2n} + O (\frac{k^{3}}{n^{2}})\}

• Now, \frac{k^{3}}{n^{2}} → 0 faster than \frac{k^{2}}{n} as n → ∞ and the Fourier transform of a Gaussian, is well known to also be a Gaussian,

• so we find probability distribution:

p(y) = \frac{1}{\sqrt{2π}} \frac{\sqrt{n}}{σ} exp \{−\frac{n(x − μ)^{2}}{2σ^{2}}\}

• and y is Gaussianly distributed as claimed.

CPS713 STAT G.C.Fox, Feb 96  49  gcf@npac.syr.edu

Monte-Carlo Integration and

the Central Limit Theorem—I

• Let’s consider Monte-Carlo Integration—an important numerical technique—from a formal point of view.

• We want to evaluate

I = b ; ∫ ; a f(x)dx

• We can surely write

I = (b − a) +∞ ; ∫ ; −∞f(x)p(x)dx

• where

p(x) = \frac{1}{(b − a)}  a ≤ x ≤ b
= 0  outside [a, b].

• and

∫ p(x) dx = 1

Statistics for Physics  1: Basic Definitions and Results

Monte-Carlo Integration and

the Central Limit Theorem—II

• Let \overline{I} = \frac{1}{n} n ; ∑ ; i=1f (x_{i}) where x_{i} are distributed according to p(x).

• Now apply central limit theorem to z_{i} = f (x_{i}).

• z_{i} are independent identically distributed random variables with

• mean μ = b ; ∫ ; af(x)p(x)dx and

• variance σ^{2} = b ; ∫ ; af^{2}(x)p(x)dx − μ^{2}.

• Then according to the Central Limit Theorem • \overline{I} → \frac{I}{(b−a)} ± \frac{σ}{\sqrt{n}}

• or I can be estimated as (b − a)\overline{I} with error (b − a)\frac{σ}{\sqrt{n}}.

Monte-Carlo Integration and

the Central Limit Theorem—III

• Clearly, a good Monte Carlo evaluation requires small σ (i.e., small “variance” of f).

• We discuss this in numerical methods notes.

• Note that in any Monte Carlo calculation, one should calculate the error as well as value of integral.

• Estimation of σ involves calculating

b ; ∫ ; a f^{2}(x)p(x)dx = ⟨f^{2}(x)⟩

• which is estimated using Central Limit Theorem as

\frac{1}{n} n ; ∑ ; i=1 f^{2} (x_{i})

• so one should always accumulate ∑f^{2} as well as ∑f when doing a Monte Carlo integral.

Extension of Central Limit Theorem

to Correlated Random Variables—I

• An important extension of the central limit theorem (necessary when discussing moments) is

• Let x_{i} be independent random variables with the same distribution.

• Let f(x), g(x) be any function of these random variables.

• Suppose we define means, standard deviation and correlations.

μ_{f} = ⟨f⟩ ,  μ_{g} = ⟨g⟩
σ_{f}^{2} = ⟨(f − μ_{f})^{2}⟩
σ_{fg} = ⟨(f − μ_{f}) (g − μ_{g})⟩
σ_{g}^{2} = ⟨(g − μ_{g})^{2}⟩
M = ⟨ \begin{matrix}σ_{f}^{2} & σ_{fg} \\ σ_{fg} & σ_{g}^{2} \end{matrix} ⟩  moment matrix

• Finally, set

F = \frac{1}{n} n ; ∑ ; i=1 f (x_{i})
G = \frac{1}{n} n ; ∑ ; i=1 g (x_{i})

Extension of Central Limit Theorem

to Correlated Random Variables—II

• Then, the new Central Limit Theorem states that the pair (F, G) is asymptotically a (two-dimensional) Gaussian distribution with probability distribution

p(F, G) ∝ exp \{−\frac{n}{2}(y − μ)^{T}B(y − μ)\}

• where we have lapsed into a vector notation:

y = [ \begin{matrix}F \\ G \end{matrix} ] , μ = [ \begin{matrix}μ_{f} \\ μ_{g} \end{matrix} ] and matrix B^{−1} = M.

• Proof: Just as in ordinary central limit theorem, except you consider two-dimensional characteristic function φ (k_{1}, k_{2}).

• This extension is important as it allows you to use correlations between functions calculated from the same data (observations).

Central Limit Theorem for Quotients

• Corollary of correlated central limit theorem: The quotient F/G is asymptotically Gaussian with mean μ_{f}/μ_{g} and standard deviation

\frac{⟨f⟩}{⟨g⟩}\sqrt{\frac{1}{n} ⟨(\frac{f}{⟨f⟩} − \frac{g}{⟨g⟩})^{2}⟩}

• Examples of above result are:

• (a) G = 1, then F/G has mean= μ_{f}, σ = \sqrt{\frac{1}{n}σ_{f}^{2}} , i.e., we recover the ordinary central limit theorem

• (b) F = G, then F/G has mean= 1, σ = 0. This is what one would hope.

Statistics for Physics  1: Basic Definitions and Results

(i) Example from Histogramming

(Physics) Events—I

• Suppose we select a sample of events from some set,

• e.g., the sample may be those falling into some bin of a histogram.

• We wish to estimate fraction of events with given property.

• Define random variable f by

f = 1 event passes selection criteria
= 0 event fails criterion.

• Write f = θg • g = 1 Event in Given Histogram

• g = 0 Event Not in Histogram

• θ is only defined for events that fall in histogram, and

• θ = 1 event in particular bin

• θ = 0 event does not fall in bin

• Let ⟨θ⟩ = p, and ⟨f⟩ = p⟨g⟩

CPS713 STAT G.C.Fox, Feb 96  56  gcf@npac.syr.edu

(i) Example from Histogramming

(Physics) Events—II

• Here, we have a hierarchical situation, where we have a total N events of which N_{hist} fall in the particular histogram.

\frac{F}{G} = \frac{N ; ∑ ; i=1 θ_{i} g_{i}}{N ; ∑ ; i=1 g_{i}}

= \frac{Number of Events in Bin}{Number of Events in Histogram}

• This estimate of fraction of events in histogram has standard deviation

\frac{⟨f⟩}{⟨g⟩} \sqrt{\frac{1}{N} ⟨(\frac{f}{⟨f⟩} − \frac{g}{⟨g⟩})^{2}⟩}

= p \sqrt{\frac{1}{N ⟨f⟩^{2}} ⟨(gθ − g p)^{2}⟩}
= p \sqrt{\frac{⟨g^{2}⟩}{N⟨f⟩^{2}} ⟨(θ − p)^{2}⟩}

• θ is binomial and ⟨θ^{2}⟩ = ⟨θ⟩ = p

• ⟨(θ − p)^{2}⟩ = p(1 − p)

Statistics for Physics  1: Basic Definitions and Results

(i) Example from Histogramming

(Physics) Events—III

• Also, ⟨g^{2}⟩ = ⟨g⟩

• And so error of estimate is \sqrt{p(1 − p)}/\sqrt{N⟨g⟩}

• where we can estimate

• \sqrt{N⟨g⟩} = N_{hist}− the number of events in histogram.

• Actually, we can “derive” this result more simply using just the single indicator variable θ on the events that fall in the histogram.

• This is a little flaky because the value N_{hist} is a random variable (with an error ±\sqrt{N_{hist}}) as it comes from selecting events from the larger sample N.

• The tortuous argument above can convince you that it doesn’t matter that N_{hist} is a random variable!

• The naive argument which gives correct answer is

⟨θ⟩ = \frac{1}{N_{hist}} N_{hist} ; ∑ ; i=1 θ (x_{i})
⟨θ^{2}⟩ = \frac{1}{N_{hist}} N_{hist} ; ∑ ; i=1 θ^{2} (x_{i}) = ⟨θ⟩
σ = ⟨θ^{2}⟩ − ⟨θ⟩^{2} = p − p^{2}

• and so the error is σ/\sqrt{N} in agreement with above (σ is standard deviation of θ).

CPS713 STAT G.C.Fox, Feb 96  58  gcf@npac.syr.edu

(ii) Analysis of Scattering

Experiments—I

• Suppose a scattering experiment observes θ with probability p(θ).

• Then typically, one would expand p(θ) = ∞ ; ∑ ; n=0a_{n}P_{n}(cos θ) as a series of Legendre Polynomials P_{n}(cos θ).

• Then,

⟨P_{n_{1}}(cos θ)⟩ = ∫ P_{n_{1}}(cos θ)p(θ)dθ
= a_{n_{1}} \frac{2}{2n + 1}

• Therefore, a_{n} is proportional to ⟨P_{n_{1}}(cos θ)⟩, which needs to be calculated as averages over observed values, θ_{k}.

⟨P_{n}(cos θ)⟩ = \frac{1}{N} N ; ∑ ; k=1 P_{n} (cos θ_{k})

(ii) Analysis of Scattering

Experiments—II

• Even though P_{n_{1}}(cos θ) and P_{n_{2}}(cos θ) are orthogonal polynomials, the means ⟨P_{n_{1}}(cos θ)⟩, ⟨P_{n_{2}}(cos θ)⟩ are not uncorrelated.

• Both are calculated from same set {θ_{k}} and fluctuate in a correlated fashion as observations θ_{k} fluctuate.

• We can calculate correlations

f = P_{n_{1}}(cos θ)
g = P_{n_{2}}(cos θ)
F = \frac{1}{N} _{k}^{∑} P_{n_{1}} (cos θ_{k})
G = \frac{1}{N} _{k}^{∑} P_{n_{2}} (cos θ_{k})

• F and G will be independent if ∫ p(θ)P_{n_{1}}(cos θ)P_{n_{2}}(cos θ) = 0.

• This is not in general true.

• Thus, our formalism gives correlation between a_{n_{1}} and a_{n_{2}}.

Statistics for Physics  1: Basic Definitions and Results

(iii) Example of Weighted Histograms

• Consider a histogram in which N events with weight w_{i}(i = 1 … N) and x-value x_{i} fall.

• This occurs in Monte Carlo calculation of ^{\frac{∫ x w(x)dx}{∫ w(x)dx}}.

⟨x⟩ = \frac{N ; ∑ ; i=1w_{i} x_{i} = S_{3}}{N ; ∑ ; i=1w_{i} = S_{1}}

• with f = wx, g = w.

• The error in this calculation is

σ = \frac{S_{3}}{S_{1}}\sqrt{\frac{1}{N} ⟨(\frac{Nwx}{S_{3}} − \frac{Nw}{S_{1}})^{2}⟩}
σ^{2} = \frac{S_{3}^{2}}{S_{1}^{2}} N ⟨[\frac{w^{2}x^{2}}{S_{3}^{2}} − \frac{2w^{2}x}{S_{3}S_{1}} + \frac{w^{2}}{S_{1}^{2}}]⟩
= \frac{1}{S_{1}^{2}} [S_{6} − 2⟨x⟩S_{4} + ⟨x⟩^{2}S_{2}]
S_{2} = N ; ∑ ; i=1 w_{i}^{2}
S_{4} = N ; ∑ ; i=1 w_{i}^{2} x_{i}
S_{6} = N ; ∑ ; i=1 w_{i}^{2} x_{i}^{2}

CPS713 STAT G.C.Fox, Feb 96  61  gcf@npac.syr.edu

2: FEEFOEFUM IN LAS VEGAS: COMBINATORIC AND DISCRETE

Statistics for Physics  2: Feefoefum in Las Vegas

PROBABILITY OVERVIEW Combinatoric and Discrete

Probability Overview

2.1 Binomial & Poisson Distributions

• We now come to a change of pace and consider more combinatorial aspects of probability.

• One of the nicest references for this is a two-volume book by Feller, who has a host of beautiful examples.

Bernoulli Trials—I

• A series of Bernoulli Trials is the generalized coin tossing problem.

• We have a sequence of independent “trials”—each of which has only two outcomes:

• Success or “heads”: probability p

• Failure or “tails”: probability q

where, of course p + q = 1.

• Afficiando’s of affirmative action and women’s lib will note that generations of math books have downgraded tails compared to heads.

• Anyhow, we follow in their footsteps and set b(k, n, p) as the probability of k successes in n trials.

Bernoulli Trials II

• b(k, n, p) is easily calculated by writing out all possible outcomes as, for instance,

\underbrace{ht h tt …}_{n}

• There are ( \begin{matrix}n \\ k \end{matrix} ) configurations with k “h” ’s and

each has probability p^{k}q^{n−k}. So:

b(k, n, p) = ( \begin{matrix}n \\ k \end{matrix} ) p^{k}q^{n−k}  (13)

• This is the binomial distribution because b(k, n, p) are just terms in binomial expansion

1 = (p + q)^{n} = n ; ∑ ; k=0 b(k, n, p).

Example of Bernoulli Trials

• A course consists of 50 lectures given to 20 people, all of whom get up before 9:00 a.m. and attend all classes.

• A grand total of 1,000 smiles is observed.

• What is the chance that Frank Nagy (a physics graduate student in times gone by) smiled exactly once during the whole course?

• Solution: Use previous formalism, setting each trial as a smile: and success as Frank Nagy smiled.

• A priori, each smile has equal probability of belonging to each person so p = 1/20, q = 19/20, and n = 1000.

• Thus: P r{Frank Nagy smiles once and once only}

= ( \begin{matrix}1000 \\ 1 \end{matrix} ) · \frac{1}{20} · (\frac{19}{20})^{999}
= 50 e^{999 log_{e} 19/20} which isn’t very big
≈ 50 e^{−51}
≈ 10^{−19}

Poisson Limit of Binomial Distribution

• The Poisson Limit is a very important special case of the binomial distribution, and it describes anything from telephone calls to events observed in particle physics.

• Take limit in b(k, n, p) of n → ∞ and p → 0 such that np = λ remains finite.

• k is also fixed.

b (k, n, \frac{λ}{n}) = \frac{n!}{k!(n − k)!} (\frac{λ}{n})^{k} (1 − \frac{λ}{n})^{n}

• Now \frac{n!}{(n−k)!} = n(n − 1) … (n − k + 1) → n^{k} as n → ∞ for fixed k. (You can also get this using Stirling’s formula for the gamma function).

• Also,

(1 − \frac{λ}{n})^{n−k} ≈ (1 − \frac{λ}{n})^{n} → e^{−λ}. So,
b (k, n, \frac{λ}{n}) → \frac{n^{k}}{k!} \frac{λ^{k}}{n^{k}} e^{−λ}
= P (k, λ) = e^{−λ} \frac{λ^{k}}{k!}

• which is the Poisson Distribution P (k, λ) of k events.

• Note: ∞ ; ∑ ; k=0P (k, λ) = 1 as it should.

Example of Poisson Distribution

• Given circumstances of previous Binomial example, what is chance that Frank Nagy smiles exactly once in a given lecture:

• Solution: Again, we set each trial = a smile and we have 1,000 of them. Now Success = “Frank Nagy smiled this smile in our given lecture” and it is not very probable p = \frac{1}{20.50} = \frac{1}{1000}. So, n is large, p is small, but np = λ = 1 is reasonable.

• So, the Poisson approximation to the exact (in this example) Binomial distribution is valid and chance is e^{−1}, i.e., it is quite probable that Frank smiles once and once only in a given lecture.

• The binomial estimate is

( \begin{matrix}1000 \\ 1 \end{matrix} ) \frac{1}{1000} (\frac{999}{1000})^{999} = [\frac{999}{1000}]^{999}

which is indeed ≈ e^{−1}.

Statistics for Physics  2: Feefoefum in Las Vegas

Compound Poisson Distribution

• Let y = N ; ∑ ; i=1y_{i}

• y_{i} is Poisson with mean μ

• N is Poisson with mean λ

• P r(y = integer r) is ∞ ; ∑ ; N=0P r(Ntakes given value) × P r ( N ; ∑ ; i=1y_{i}(N fixed has r events)).

• As for fixed N, y is Poisson with mean Nμ, we find

P r ( y = integer r)
= ∞ ; ∑ ; N=0\frac{λ^{N}e^{−λ}}{λ!} · \frac{(Nμ)^{r}e^{−Nμ}}{r!}

• More generally, we could give y_{i} a non-Poisson distribution. One example of non-Poisson y_{i}, is:

• The physics of charged particle ionization—p.53 EDJRS

Particle Cluster Production—I

• Particularly interesting is sensitivity to assumption of which y_{i} are chosen to be Poisson.

• Note the Poisson distribution with mean (number of objects) = λ, and where each object counts four units (telephone calls, etc.), is much broader than ordinary Poisson of mean 4λ, where each object gives a single unit.

• This is of importance in particle physics where

• We believe that q\underline{q} “clusters” are produced independently in a high energy collision.

• These clusters are (assumed to be) Poissonly distributed and decay into 4 → 6 final particles.

• One can show that if particles were produced in a Poisson fashion singly, then

• f_{2} = σ^{2}− mean for final particles is zero.

• However, in cluster model where one can assume y_{i} are distributed Gaussianly with the 4 → 6 mean. f_{2} is always > 0.

• Experimentally, one finds f_{2} > 0 except when perturbed by kinematic constraints.

• i.e., cluster model and experiment is broader than naive Poisson.

Particle Cluster Production II

• To take a particularly simple case, consider a Poisson distribution of clusters—each of which always contains four particles.

• For a random variable x with a Poisson distribution

Mean = λ σ^{2} = λ
f_{2}(x) = σ_{x}^{2} − ⟨x⟩ = 0

• Now, let random variable y = 4x

⟨y⟩ = 4λ
⟨y^{2}⟩ = ⟨(4x)^{2}⟩ = 16 ⟨x^{2}⟩
= 16 (⟨x^{2}⟩ − ⟨x⟩^{2} + ⟨x⟩)
= 16 (σ_{x}^{2} + ⟨x⟩^{2})
= 16 (λ + λ^{2})
f_{2}(y) = ⟨y^{2}⟩ − ⟨y⟩^{2} − ⟨y⟩
= 16 (λ + λ^{2}) − 16λ^{2} − 4λ
= 12λ > 0

• i.e., y is always broader than a Poisson of mean 4λ, which has f_{2} = 0.

Example of Cherenkov Light

• Cherenkov Light is given off by a charged particle in certain media.

• The amount of radiation is used to distinguish different particles in physics experiments.

• Total pulse height observed y = N_{γ} ; ∑ ; i=1y_{i}.

• y_{i} = response in phototube of a single photon—probably Gaussian (i.e., limit of large Poisson distribution).

• N_{γ} = number of photons given off in Cherenkov medium which is certainly Poisson.

Additivity Property of Binomial

or Poisson Distributions

• These two distributions (Binomial and Poisson) share an important additivity property, which is

• α), obvious and β), leads to confusion when used in conjunction with central limit theorem (see later).

• (a) Suppose the independent random variables x_{1} and x_{2} have binomial distributions with parameters (n_{1}, p) and (n_{2}, p) respectively.

• Here, x_{i} = k if there are k successes.

• Then it is obvious from source of binomial distribution that x_{1} + x_{2} has a binomial distribution with parameter (n_{1} + n_{2}, p).

• (b) Suppose independent random variables y_{1} and y_{2} have Poisson distributions λ_{1} and λ_{2}.

• Then y_{1} + y_{2} has Poisson distribution with parameter λ_{1} + λ_{2}.

• This is obvious either algebraically, or by using • the result (a) above with n_{1} = λ_{2}/p and

λ_{2} = λ_{2}/p and let n_{i} → ∞, p → 0, fixed λ_{1}, λ_{2}.

Gaussian Limit of Binomial Distribution

• Matthews and Walker confused me by discussing two limits of binomial distribution.

• We have already discussed one of these—the Poisson limit—which is really a distinctive limit. However, we can also see a Gaussian limit which is just a special case of central limit theorem.

• Although this limit has a particular interpretation because of the additivity property we just discussed.

• Note also that central limit theorem was discovered in general in 1812 while de Moivre discovered its binomial special case in 1733.

• Let x_{i} be independent random variables—each of which has the same b(k, 1, p) distribution,

• i.e., for each i : x_{i} = 0 is failure and x_{i} = 1 is success with relative probabilities 1 − p and p respectively.

• Then y = N ; ∑ ; i=1x_{i} is in fact, k, the number of successes in n Bernoulli trials.

Gaussian Limit of Binomial

and Poisson Distributions—II

• Then (α) from additivity of binomial distributions; y is binomial with distribution b(k, n, p),

• (β) from standard central limit theorem, y is Gaussian with mean nμ and standard deviation \sqrt{n}σ.

• Here μ and σ are mean and standard deviation of b(k, 1, p) distribution,

• i.e.,

μ = 0 · (1 − p) + 1 · p = p
σ^{2} = ⟨x^{2}⟩ − ⟨x⟩^{2}
= 0 · (1 − p) + 1 · p − p^{2} = p(1 − p)

• We conclude that for fixed p, and n → ∞, b(k, n, p) tends to Gaussian distribution with mean np and standard deviation \sqrt{p(1 − p)n}.

CPS713 STAT G.C.Fox, Feb 96  75  gcf@npac.syr.edu

Gaussian Limit of Binomial

and Poisson Distributions—III

• We can do the same trick for the Poisson distribution.

• Fix λ finite and let x_{i} be independent with each x_{i} as Poisson random variable with parameter λ.

• Then from additivity y = x_{1} … + x_{n} has Poisson distribution with parameter nλ = Λ.

• Meanwhile, from central limit theorem, y is Gaussian with mean nμ and standard deviation \sqrt{n}σ.

• Here again μ and σ are mean and standard deviation of individual distribution.

• Here mean is given by

μ = ∑ exp(−λ)k\frac{λ^{k}}{k!}
= λe^{−λ}\frac{∂}{∂λ} (\frac{∑ λ^{k}}{k!})
= λe^{−λ}\frac{∂}{∂λ} (e^{λ})
= λ

Gaussian Limit of Binomial

and Poisson Distributions—IV

• Standard deviation of Poisson is given by

σ^{2} = ⟨x^{2}⟩ − λ^{2}
= e^{−λ} (\frac{λ∂}{∂λ})^{2} (e^{λ}) − λ^{2}
= e^{−λ} (\frac{λ∂}{∂λ}) (λe^{λ}) − λ^{2}
= λ + λ^{2} − λ^{2} = λ = μ

• So, for very large Λ, any Poissonly distributed random variable y has mean Λ and standard deviation \sqrt{Λ}.

• We will soon show that counts in any counting experiment (e.g., number of raindrops in a bowl or particles in a high energy physical experiment) is described by a Poisson random variable with (large) parameter Λ = number of expected counts.

• The above result then translates into well known result that counting experiment with N events has error \sqrt{N}.

2.2 Birth and Death

Processes Overview—I

• We will briefly sketch some techniques which at their most elementary, justify Poisson for a counting experiment.

• The more advanced applications are important in many practical problems (telephones (random calls) and biology (random birth/deaths)), involving happenings distributed more or less randomly in times.

• We start with a discrete (as opposed to continuous) example: Telephone calls.

• Suppose your favorite secretary receives N telephone calls in an hour.

• We wish to find the probability that she received k calls in a given time t

Statistics for Physics  2: Feefoefum in Las Vegas

Telephone Calls as a Birth Process—I

• Consider example given on previous page.

• Consider every second.

• ∃ probability p = N/3600 of receiving one call:

• ∃ essentially no probability of receiving > 1 call.

• Thus, each second can be described by variable y_{i}(= 0 no, = 1 one call) which has Poisson distribution with parameter λ = N/3600.

• The number of calls in time t seconds is t ; ∑ ; i=1y_{i} which from additivity is also Poisson with parameter

t ; ∑ ; i=1λ_{i} = pt.

• So

P_{k}(t) = e^{−pt}\frac{(pt)^{k}}{k!}

is probability of receiving k calls in t seconds.

Telephone Calls as a Birth Process—II

• Note mean number received is pt or for t = 3600, just N as input data claimed.

• One could also derive this by taking each second as a Bernoulli trial with probability p of success.

• Then we get distribution b(k, t, p) which has t → ∞, p → 0: pt finite, i.e., we can employ usual Poisson limit.

Basic Birth Process Formalism—I

• It is natural to take limit of time interval (second in above example) as zero, and let t be a continuous and not a discrete valued variable.

• Suppose we are studying events (e.g., telephone calls as above) that

• α) occur independently of previous history of world,

• β) in time [t, t + dt].

• The probability of an event happening is λdt + o(dt),

• and the probability of > 1 event happening is o(dt).

• Define P_{k}(t) as the probability of k events happening in time 0 to t.

Basic Birth Process Formalism—II

• Then, the probability of k events at time t + dt is

P_{k}(t + dt) = P_{k}(t)[1 − λdt] + P_{k−1}(t)λdt + o(dt)

• or

\frac{dP_{k}}{dt}(t) = −λP_{k}(t) + λP_{k−1}(t) : k ≠ 0

• For k = 0, the above is modified as P_{0−1} term doesn’t exist, i.e.,

\frac{dP_{0}(t)}{dt} = −λP_{0}(t) : k = 0

• Solving the last gives P_{0}(t) = e^{−λt} and then induction shows

P_{k}(t) = e^{−λt}\frac{(λt)^{k}}{k!}

and we have achieved Poisson again.

• It is now obvious why any counting experiments for random events should have a Poisson Distribution.

General Birth and Death Process

• Consider (Feller, p. 404) a system with a denumerable number of states E_{n} whose probability of occupation varies with t.

• For Poisson E_{n} is “we have observed n events.”

• Note only one E_{n} occupied at a time

• and ∃ chance λdt of a transition E_{n} → E_{n+1}.

• This is a pure birth process.

• We can generalize λ to λ_{n}(t): a function of n and t.

• Also, we can have death processes which also allow relapses E_{n} to E_{n−1} with chance μ_{n}(t)dt. So:

• change E_{n} → E_{n+1} chance λ_{n}dt + o(dt)
•  E_{n} → E_{n−1} chance μ_{n}dt + o(dt)
• Any other change chance o(dt)

• We can, at once, derive a similar differential equation for P_{n}(t): the probability of system being in state E_{n}.

\frac{dP_{n}(t)}{dt} = − (λ_{n} + μ_{n}) P_{n}(t)
+ λ_{n−1}P_{n−1}(t) + μ_{n+1}P_{n+1}(t)

• with μ_{0} = 0 and λ_{−1} = 0 so that above is still correct with n = 0.

CPS713 STAT G.C.Fox, Feb 96  83  gcf@npac.syr.edu

Example of Birth/Death Process—I

• Consider a telephone exchange with an infinite number of lines and let E_{n} be state “n lines are busy.”

• Then we apply above formalism with

• μ_{n} = nμ where μ is fixed (1/μ is average length of call),

• but λ_{n} = λ independent of n.

• One can find the steady state solution—independent of both t and initial state:

P_{n}(t) = e^{−λ/μ}\frac{(\frac{λ}{μ})n}{n!}

as probability of n lines being busy.

Example of Birth/Death Process—II

• It is also easy to do—the seemingly harder problem—of a finite number of telephone lines in exchange.

• There is an example of this sort in problem set about parking lots.

• Similar techniques can be used to decide on a number of servicemen needed to repair a set of machines of prescribed reliability. Here,

• Death = machine starts working

• Birth = machine breaks down, and

• E_{n} = n machine’s broken.

• There is a beautiful discussion of this on pages 416–420 of Feller.

Markov Chains

• One can generalize either the discrete or continuous time cases by considering a system with again states E_{n} and now not

• only transitions E_{n} → E_{n±1}

• but a general probability p_{jk} (discrete) or p_{jk}dt (continuous) of transitions from E_{j} to E_{k} at each trial.

• This theory includes general random walk problem and is discussed in Volume 1, Feller.

• The mathematics is not very deep—solving recurrence relations—and the techniques are not useful for the basic problem in physics (statistics of data reduction).

• So I will not consider further. People with biological bent should, however, read further in Feller.

• This Markov theory is relevant in discussions of the so called Metropolis Method in Monte Carlo integration.

3: ESTIMATION OF PARAMETERS

Statistics for Physics  3: Estimation of Parameters

Estimation of Parameters—References

• These notes (Estimation of Parameters) refer to some additional references:

• F. T. Solmitz, 1964, Ann. Rev. Nucl. Sci, 14, 375.

• and some old notes I have copies of:

• J. Orear, “Notes on Statistics for Physicists,” UCRL-8417 (1958).

• S. Yellin, DESY 72/12: Preprint (1972).

• K. McDonald; CTSL Internal Report #57.

Statistics for Physics  3: Estimation of Parameters

3 Overview of Estimation

of Parameters

• So far, we have been concerned with rather straight-forward (academic) problems involving random variables with known probability distributions.

• Now we turn to a fundamental scientific problem:

• Given some observations x_{1} … x_{n} (= values of a random variable) of some physical quantity,

• how do we extract theoretical parameters describing systems we have observed, and how do we estimate the error in our determination.

• This will require rather more careful thought as to the relation of probability theory and the real world.

• We will first describe the principle of Maximum Likelihood—everything else can be derived from this.

3.1 Maximum Likelihood Principle

Illustrated by Measurement Errors—I

• Take the example of a rod quivering in an earthquake

• x_{1}, x_{2} … x_{n} are n measurements of the length of a rod.

• Owing to a Gaussian earthquake (or a Gaussian, i.e., normal drunkenness) occurring during the measurement, the x_{i} are scattered about the true value.

• You visit your favorite theorist (or do it yourself after sobering up) who believes that the probability distribution of x_{i} is Gaussian with mean α_{0} = true length and standard deviation σ_{0}.

• Our superstar knows the value of σ_{0} but not of α_{0}.

• We may write down the (relative) probability of each measurement x_{i} as:

p_{i} (x_{i}) = c exp \{−\frac{1}{2σ_{0}^{2}} (x_{i} − α_{0})^{2}\}

Maximum Likelihood Principle

Illustrated by Measurement Errors—II

• The likelihood function for whole experiment of n events is defined as:

L (x_{1} … x_{n}) = n ; ∏ ; i=1 c exp \{−\frac{1}{2σ_{0}^{2}} (x_{i} − α_{0})^{2}\}
= c̃exp \{−\frac{1}{2σ_{0}^{2}} _{i}^{∑} (x_{i} − α_{0})^{2}\}

• c and c̃ are calculable constants to normalize probability. This is—up to a constant—the probability of the x_{i} given,

• the true length α_{0},

• and certain theoretical assumptions or prejudices that we represent as T.

L (x_{i} … x_{n}) = P r {x_{1} … x_{n}|α_{0}, T }

using conditional probability notation.

• T is, for instance, the notion that distribution is Gaussian and the value of σ_{0} taken (the latter being calculable from the strength of earthquake or amount of beer drunk).

Enter Mr. Bayes

for Conditional Probablities—I

• Now what we really want is the probability distribution of an estimate α of α_{0}.

• We can formally get this by use of Bayes law:

• First extend p_{i} (x_{i}) = exp \{−\frac{1}{2σ_{0}^{2}} (x_{i} − α)^{2}\} to be defined for any value of α; and not just for α = true value α_{0}.

• What does this mean, as only one true α_{0}?

• Introduce: P (α|x_{i}, T) = Probability that true length is α given measurements x_{i} and prejudice T.

• Then by Bayes—or rather a slight extension of it,

P (α|x_{i}, T) = P (x_{i}|α, T) \frac{P (α|T)}{P (x_{i}|T)}  (14)

• which only differs from the ordinary Bayes Law by extra conditional |T in each probability.

Enter Mr. Bayes

for Conditional Probablities—II

• Equation (14) is the essence of the maximum likelihood principle.

P (α|x_{i}, T) = P (x_{i}|α, T) \frac{P (α|T)}{P (x_{i}|T)}

• Equation (14) converts the likelihood L (x_{1} … x_{n}) = P (x_{i}|α, T)

• which is a function or probability of x_{i} for fixed α, T

• into P (α|x_{i}, T) which is the desired probability distribution of α for given measurements x_{i} and prejudices T.

Interpretation of Bayes Formulation

of Likelihood—I

• Consider Equation (14)

P (α|x_{i}, T) = P (x_{i}|α, T) \frac{P (α|T)}{P (x_{i}|T)}

• (α ) We can ignore P (x_{i}|T) because it is independent of α and fixed for our given x_{i} and T.

• This has the practical consequence that the likelihood is never normalized sensibly.

• (β ) P (α|T) philosophically represents a priori knowledge of α before experiment (measurement of x_{i}) performed.

• This is shaky—not just because we don’t really know how to define what |T (“given T ”) means

• but, for instance, suppose we wish to parameterize the simplest idea—“complete ignorance” of α before experiment performed.

Interpretation of Bayes Formulation

of Likelihood—II

• Natural parameterization of ignorance is P (α|T) = constant.

• But suppose we put β = α^{2} and then β has a priori probability ϕ(β/T)

ϕ(β|T)dβ = P (α|T)dα, i.e.,
ϕ(β|T) = \frac{1}{2\sqrt{β}} ≠ constant

• So the concept of “complete ignorance” is not well defined, as change of variables changes concept.

• Jeffreys (a big bad English-Oxford University Press book cited in Mathews and Walker) spends 100’s of pages trying to avoid this difficulty.

• His efforts are, in my opinion, so much garbage.

• Fortunately, there is no trouble for large n because of both an obvious numerical reason (see next slide (The Likelihood Maximizes) and a not so obvious mathematical reason (see slide (General Maximum Likelihood Methods)).

The Likelihood Maximizes the

Quivering Rod—I

• Consider, again, the explicit form of the likelihood of a rod quivering in the earthquake.

L (x_{1} … x_{n}) = c exp \{−\frac{1}{2σ^{2}} n ; ∑ ; i=1 (x_{i} − α)\}^{2}

• According to Equation (14), we now regard L as a function of α and so rewrite it:

L(α) = c exp \{−\frac{n}{2σ^{2}} (α − \frac{1}{n} n ; ∑ ; i=1x_{i})^{2}
− \frac{1}{2σ^{2}} [_{i}^{∑} x_{i}^{2} − \frac{1}{n} (Σ_{i}x_{i})^{2}]\}

• On completing the square to write exponent as quadratic in α.

Statistics for Physics  3: Estimation of Parameters

The Likelihood Maximizes the Quivering Rod—II

• The last term is a constant independent of α, and • so viewed as a function of α, L is just a Gaussian which maximizes at α = \frac{1}{n} n ; ∑ ; i=1x_{i} and has width \frac{σ}{\sqrt{n}}.

Statistics for Physics  3: Estimation of Parameters

The Likelihood Maximizes the

Quivering Rod—III

• So as L = P r (α|x_{i}, T)

• we say that maximum likelihood estimate of α is

α^{∗} = \frac{1}{n} n ; ∑ ; i=1x_{i}

• and its standard deviation is σ/\sqrt{n}.

• Note that whatever a priori P (α|T) you put in,

• the shape of L · P (α|T) is essentially independent of P (α|T) as n → ∞

• and we get a Gaussian peaked at α^{∗} riding on top of the constant value P (α^{∗}|T).

General Maximum Likelihood

\underline{Methods—I}

• The above is glib in many ways—especially in its treatment of errors.

• We will repair this later with real live theorems.

• First, we state the obvious generalization.

• The general maximum likelihood method is to put P (α|T) = 1 in Equation (14)

• and form L(α) = Π_{i=1}^{n}p (x_{i}, α)

• for any measurement x_{i} whose probability p (x_{i}, α) depends on value x_{i} and theoretical parameter(s) α.

• The maximum likelihood estimate of α is the value α^{∗} where L(α^{∗}) is a maximum.

• The error in estimate is the standard deviation of the distribution of L as a function of α, i.e.,

L(α) = exp \{const. − \frac{1}{2Δα^{2}} (α − α^{∗})^{2}\}

Statistics for Physics  3: Estimation of Parameters

General Maximum Likelihood

Methods—II

• One must issue some obvious warnings.

• For instance, L could have several maxima

• For instance, this would happen if in earthquake problem you were still drunk in morning and decided to find \sqrt{α}.

• Then \sqrt{α} would peak at ±\sqrt{α^{∗}}.

• Less obviously one should remember that theorems (to come) only say that maximum likelihood method works for large n.

• Then it is a universally good method.

• However, for small n it is much more likely to give nonsense than many other methods. Thus, …

CPS713 STAT G.C.Fox, Feb 96  100  gcf@npac.syr.edu

Statistics for Physics  3: Estimation of Parameters

Warning Example of Misuse of Likelihood Method

• Consider probability distribution appropriate for polarization experiment:

p(φ, α) = \frac{(1 + α cos φ)}{2π} : −1 ≤ α ≤ 1.

• If we had but one event (a very small experiment!)

L(α) = \frac{(1 + α cos φ_{1})}{2π}

• if cos φ_{1} > 0 : α = 1 (i.e., end of range)

maximizes L

• or cos φ_{1} < 0 : α = −1 (other end of range)

maximizes L

• If you have two events:

L(α) = \frac{(1 + α(cos φ_{1} + cos φ_{2}) + α^{2} cos φ_{1} cos φ_{2})}{2π}

• This has a maximum at end of range, again, unless cos φ_{1} cos φ_{2} < 0.

Practical Advice on Use

of Maximum Likelihood—I

• So usually people warn you by saying:

• “If distribution of α not Gaussian, then do not simply quote maximum and standard deviation, rather plot the likelihood.”

• Unfortunately, this is dubious, because for small n, nobody has proved any theorems and P (α|T) philosophy important

• while for large n, theorems are proven on the relevance of the likelihood method

• but they also prove L is Gaussian!

Practical Advice on Use

of Maximum Likelihood—II

• So correct warning is: “If distribution not Gaussian—everything is OK if either:

• (α) Just many sharp humps as in \sqrt{α} example, data can be fitted with > 1 value of parameter.

• (β) If several parameters, can be Gaussian in some and flat in others if depends weakly (or not at all!) on them.

• (γ) Gaussians maybe distorted by effects due to end of range.

• If (α), (β), or (γ) don’t apply and distribution of α isn’t Gaussian, then one can only be unhappy!

• You should apply for a new grant and do a better experiment!

Determination of Lifetime Example

of Maximum Likelihood—I

• All books give the example of determining the lifetime τ_{0} of a particle from the observation of decay times t_{1} … t_{n} of n decays.

• The probability distribution for decay problem is Poisson:

p(t)dt = \frac{1}{τ} exp (−\frac{t}{τ} ) dt

• Now, dt is meaningless—it is independent of τ and only affects normalization of likelihood.

• The unknown P (α|T), P (x_{i}|T) terms in Equation (14) made this arbitrary anyway.

• We must, however, keep terms in normalization (\frac{1}{τ} here) that depend on parameters to be determined.

• The likelihood of n decays is:

L(τ) ∝ \frac{1}{τ} e^{−\frac{t_{1}}{τ}} … \frac{1}{τ} e^{−\frac{tn}{τ}}, or
L(τ) = exp [−\frac{1}{τ} n ; ∑ ; i=1t_{i} − n ln τ]

where we lazily convert ∝ into an = sign.

Determination of Lifetime Example

of Maximum Likelihood—II

• The maximum of L is given by:

\frac{dL}{dτ} = 0 at τ = τ^{∗} or
\frac{1}{τ^{∗^{2}}} n ; ∑ ; i=1t_{i} − \frac{n}{τ^{∗}} = 0 or
τ^{∗} = \frac{1}{n} n ; ∑ ; i=1t_{i}

• or the maximum likelihood estimate is just the mean:

• such a simple result is only true because of particular Poisson distribution used in the example.

Central Limit Estimate of Error

• We can use this example to illustrate an interesting theoretical point about the error, and now compare two methods of finding an error.

• (a) According to the central limit theorem, the variable y = \frac{1}{n} n ; ∑ ; i=1t_{i} is a sum over independent random variables and has mean ⟨y⟩ = τ_{0} and standard deviation σ.

σ^{2} = \frac{1}{n} \{⟨t_{i}^{2}⟩ − ⟨t_{i}⟩^{2}\}
= \frac{1}{n} \{\frac{1}{τ_{0}} ∫ t^{2}e^{−\frac{t}{τ_{0}}} dt − τ_{0}^{2}\} or
σ = \frac{τ_{0}}{\sqrt{n}}  (15)

• The mean and standard deviation of y refer to distribution of values of y in many experiments—each of n observations.

Maximum Likelihood Estimate of

Error in Lifetime—I

• (b) However, maximum likelihood method of finding error is, off hand, quite different.

• Rather than looking at distribution of y over infinitely many experiments each of n observations,

• It considers τ dependence of L for one given experiment of n observations.

Maximum Likelihood Estimate of

Error in Lifetime—II

• Then we found L became Gaussian in τ, and calling maximum likelihood standard deviation σ̃

L(τ) ∝ exp \{−\frac{1}{2σ̃^{2}}(τ − τ^{∗})^{2}\}

• where we can find width σ̃—width of Gaussian by:

\frac{1}{σ̃^{2}} = −\frac{d}{dτ^{2}} ln L(τ)|_{τ=τ^{∗}}
= \frac{2}{τ^{3}} n ; ∑ ; i=1 t_{i} − \frac{n}{τ^{2}}|_{τ=τ^{∗}} or
σ̃ = \frac{τ^{∗}}{\sqrt{n}}  (16)

• as τ^{∗} = τ_{0} + O (\frac{1}{n}) we see that Equation (15) and (16) are in fact the same asymptotically (large n, that is).

Comparison of Two Error

Estimates in Maximum Likelihood—I

• Note that physics books (Solmitz is an exception) don’t mention method (a) of finding error in maximum likelihood method.

• On the other hand, mathematics books never mention (b) as a serious method.

• They do, however, prove (a) and (b) are the same for large n (as we indicated in example above) for general distribution.

• As far as I can see, method (a) is only correct definition of error.

• It follows that method (b) is wrong when it differs from method (a).

• (b) does differ from (a) for small n when, in particular, likelihood isn’t Gaussian.

• So again warnings for small n (“plot likelihood”) are very dubious.

Comparison of Two Error

Estimates in Maximum Likelihood—II

• The correct statement is that your results are simply unreliable for small n, and error is not calculable from likelihood function alone.

• The error is defined for small n, but must be calculated as in (a)—this is general as not possible to do with just the single experiment with only a few (n) observations.

• Note that when math books discuss maximum likelihood estimate for small n—it is only the value α^{∗} they discuss and not the shape of L and the method (b) of finding error.

Statistics for Physics  3: Estimation of Parameters

Theorem on Asymptotic Validity

of Maximum Likelihood

Method—Conditions and

Statement of Theorem

• This is described in Cramer Chapters 22 and 23 • Theorem: Let P (x_{i}, α) be a normalized probability for a result x_{i} given theoretical parameter α. We have n such observations x_{i}:

• L = Π_{i=1}^{n}P (x_{i}, α) is usual likelihood function.

• Let α^{∗} be maximum of L with respect to α.

• There are three basic results that can be proven.

Theorem on Asymptotic Validity

of Maximum Likelihood Method

Conditions—First Assertion

• (a) If we fix n and conduct many such series of n observations, the mean of α^{∗} satisfies:

⟨α^{∗}⟩ = α_{0} + O (\frac{1}{n})

• Note ⟨α^{∗}⟩ is not = α_{0} exactly (α_{0} is true theoretical parameter).

• If for any estimate ⟨α_{estimate}⟩ = α_{0}, then we say estimate is unbiased:

• if ⟨α_{estimate}⟩ = α_{0} + b(α_{0}) then we have a biased estimate with bias b(α_{0}).

• So maximum likelihood estimate is biased for finite n and is only asymptotically unbiased.

• However, bias is smaller than Standard Deviation and so it is of no real importance.

• Bias is discussed in detail by Cramer.

Theorem on Asymptotic Validity

of Maximum Likelihood Method

Conditions—Second and Third

\underline{Assertions}

• (b) We can regard α^{∗} as a random variable: then its standard deviation is O (\frac{1}{\sqrt{n}}) and to this order we can correctly calculate standard deviation in the “physics way” already discussed (i.e., from α dependence of L(α)).

• (c) There is no other estimate of α_{0} which has a smaller standard deviation asymptotically in n, i.e., the maximum likelihood method is the best as long as n is large.

Statistics for Physics  3: Estimation of Parameters

Proof of Maximum Likelihood Theorem Part (a)—I

• Can be found in Cramer, p. 479 or Hoel p. 379.

• For (a) and (b), write

ln L(α) = n ; ∑ ; i=1 ln P (x_{i}, α)

• Expand L about true theoretical value.

ln L(α) = L_{0}+(α−α_{0})L_{1}+\frac{(α − α_{0})^{2}}{2}L_{2}+\frac{(α − α_{0})^{3}}{6}L_{3}+…

L_{0} = n ; ∑ ; i=1 ln P (x_{i}, α_{0})
L_{1} = n ; ∑ ; i=1 \frac{∂}{∂α} ln P (x_{i}, α_{0})
L_{2} = n ; ∑ ; i=1 \frac{∂^{2}}{∂α^{2}} ln P (x_{i}, α_{0}), etc.

• We get maximum likelihood estimate, by calculating ∂ ln L/∂α^{∗} = 0 or

L_{1} + (α^{∗} − α_{0})L_{2} = 0 i.e.,
α^{∗} = α_{0} − \frac{L_{1}}{L_{2}}  (17)
CPS713 STAT G.C.Fox, Feb 96  114  gcf@npac.syr.edu

Statistics for Physics  3: Estimation of Parameters

Proof of Maximum Likelihood

Theorem Part (a)—II

• The error in this is O \{(α^{∗} − α_{0})^{2}\frac{L_{3}}{L_{2}} \} which is O (\frac{1}{n}) as we will soon show that (α^{∗} − α_{0})^{2} is O (\frac{1}{n}).

• Now L_{1} and L_{2} are random variables and Equation (17) defines α^{∗} as a random variable.

• We can calculate from (17):

⟨α^{∗}⟩ = α_{0} − ⟨L_{1}/L_{2}⟩
= α_{0} − \frac{⟨\frac{∂}{∂α} ln P (x, α_{0})⟩}{⟨\frac{∂^{2}}{∂α^{2}} ln P (x, α_{0})⟩}

where we use theorem on F/G, proved at end of Section 2.

CPS713 STAT G.C.Fox, Feb 96  115  gcf@npac.syr.edu

Proof of Maximum Likelihood

Theorem Part (a)—III

• Note, this form for ⟨α^{∗}⟩ is only valid to O(1/n).

⟨α^{∗}⟩ = α_{0} − \frac{⟨\frac{∂}{∂α} ln P (x, α_{0})⟩}{⟨\frac{∂^{2}}{∂α^{2}} ln P (x, α_{0})⟩}

• Now, the numerator of formula for ⟨α^{∗}⟩ is ⟨\frac{∂}{∂α} ln P (x, α_{0})⟩ which is by definition

∫ dxP (x, α_{0})\frac{∂}{∂α} ln P (x, α_{0}) = ∫ dx\frac{∂}{∂α}P (x, α)|_{α=α_{0}}
= \frac{∂}{∂α} [∫ dxP (x, α)] |_{α=α_{0}}
= \frac{∂}{∂α} · 1
= 0

as P is truely normalized.

• This proves part (a) that ⟨α^{∗}⟩ = α_{0} with corrections from both of order L_{3}/L_{2} and L_{1}/L_{2} ≠ quite 0 terms which are both of order 1/n.

• These corrections have been calculated by Yellin.

Statistics for Physics  3: Estimation of Parameters

Proof of Maximum Likelihood

Part (b): Formula for Error—I

• Now to prove part (b).

• First consider “physics” method.

• This takes standard deviation of α^{∗} as \sqrt{−\frac{1}{L_{2}}} or rather this evaluated at α = α^{∗} not α = α_{0} but these are same asymptotically.

σ_{physics}^{2} = \frac{−1}{n ; ∑ ; i=1 [\frac{∂^{2}}{∂α^{2}} ln P (x_{i}, α_{0})]}
⟨σ_{physics}^{2}⟩ = \frac{−1}{n ⟨\frac{∂^{2}}{∂α^{2}} ln P (x_{i}, α_{0})⟩}  where
⟨\frac{∂^{2}}{∂α^{2}} ln P (x, α)⟩ = ∫ dxP \{\frac{−(∂P/∂α)^{2}}{P^{2}} + \frac{∂^{2}P/∂α^{2}}{P}\}

• Now second term is just \frac{∂^{2}}{∂α^{2}} ∫ P and is zero again because of normalization condition. So,

⟨σ_{physics}^{2}⟩ = \frac{1}{\{n∫ dx (\frac{∂P}{∂α} )^{2} /P \}}  (18)

CPS713 STAT G.C.Fox, Feb 96  117  gcf@npac.syr.edu

Proof of Maximum Likelihood

Part (b): Formula for Error—II

• We can calculate the correct statistical standard deviation in α^{∗} by taking Equation (17) and using our F/G theorem with f_{i} = \frac{∂}{∂α} ln P (x_{i}, α), g_{i} = \frac{∂^{2}}{∂α^{2}} ln P (x_{i}, α). Generally

σ^{2} = \frac{1}{n} ⟨(\frac{f}{⟨g⟩} − \frac{g⟨f⟩}{⟨g⟩^{2}})^{2}⟩

• but ⟨f⟩ = 0—proved above, so: σ^{2} = \frac{1}{n} \frac{⟨f^{2}⟩}{⟨g⟩^{2}}.

• It is easy to see that

⟨f^{2}⟩ = ∫ \frac{(\frac{∂P}{∂α} )2}{P^{2}} · P dx

and this is ⟨g⟩.

Statistics for Physics  3: Estimation of Parameters

Proof of Maximum Likelihood

Part (b): Formula for Error—III

• So σ_{statistical}^{2} is just \frac{1}{n⟨g⟩} which is exactly equal to ⟨σ_{physics}^{2}⟩ from Equation (18).

• The difference between σ_{physis}^{2} and ⟨σ_{physics}^{2}⟩ is no sweat:

• Thus, from central limit theorem ⟨σ_{physics}^{2}⟩ and σ_{physics}^{2} differ only by relative error O (\frac{1}{\sqrt{n}}) which tends to 0 as n → ∞.

• One can normally ignore such “errors” in the error estimate.

Multiplication of Experiments

• If we have two experiments that are independent but depend on the same theoretical parameter α. Then

L_{1} = n_{1} ; ∏ ; i=1 p(x_{i}, α) is likelihood of experiment # 1
L_{2} = n_{2} ; ∏ ; i=1 q(y_{i}, α) is likelihood of experiment # 2

• and clearly the combined likelihood is

L = L_{1}L_{2}

• This final L is correct expression of results of both experiments.

• i.e., one should multiply likelihoods when combining experiments.

• We will now apply this to case when each L_{i} is Gaussian in theoretical parameter α.

3.2 Overview of χ^{2} Method

• This is standard and (often) fully correct (best) method for combining results from several experiments.

• It is, in fact, used (because it is easy conceptually and for computer) even when it is inapplicable,

• e.g., if there are correlations between supposed “independent” random variables

• or if the basic experiments have too few events to be Gaussian.

• The problem to be solved is: given N experimental measurements x_{i} and errors σ_{i} plus estimates ξ_{i} of the x_{i}, which are functions of n_{T} theoretical parameters.

• Find best values of theoretical parameters.

Statistics for Physics  3: Estimation of Parameters

Derivation of χ^{2} from Maximum Likelihood

• Then the probability of each observation is

p_{i} = \frac{1}{\sqrt{2π}σ_{i}} exp \{\frac{−(x_{i} − ξ_{i})^{2}}{2σ^{2}}\}

• The likelihood

L = N ; ∏ ; i=1 p_{i} ∝ exp \{− N ; ∑ ; i=1 \frac{(x_{i} − ξ_{i})^{2}}{2σ_{i}^{2}}\}

• so, the principle of maximizing the likelihood is equivalent to minimizing χ^{2} where

χ^{2} = N ; ∑ ; i=1 \frac{(x_{i} − ξ_{i})^{2}}{σ_{i}^{2}}  (19)

• with respect to n_{T} theoretical parameters in ξ_{i}.

CPS713 STAT G.C.Fox, Feb 96  122  gcf@npac.syr.edu

Mean and Standard Deviation of χ^{2}

• Notice that χ^{2} is the sum N ; ∑ ; i=1v_{i}^{2} where v_{i} are independent Gaussian random variables of unit standard deviation and zero mean.

• We can calculate the χ^{2} distribution exactly (Mathews and Walker pages 393–395) but in practice this is rarely necessary.

• Thus, we can easily show that—for fixed ξ_{i}—that

⟨χ^{2}⟩ = N σ_{χ^{2}}^{2} = 2N (20)

and this is sufficient for large N as central limit theorem says χ^{2} will be Gaussian with above mean and standard deviation.

• Now Equation (20) is not quite correct because if we determine the theoretical parameters to minimize χ^{2}, then implicitly ξ_{i} are also random variables that are functions of x_{j}.

• It is not too hard to show that in this case one should just replace N by N − n_{T} (number of degrees of freedom) in (20).

• In any case, (20) is very important because it allows an easy criterion for the goodness of fit.

Interpretation of Value of χ^{2}

• Thus, if value of χ^{2} is many (say > 3) standard deviations away from its mean, then the theoretical form is in doubt.

• The usual maximum likelihood method does not have this advantage: remember L was unnormalized and so its value had no significance for the goodness of fit.

• This is quite a difficult problem and an elementary discussion is given in MacDonald notes, page 14 onwards.

• The only way to get a feeling for χ^{2} is to use it: so I won’t dwell here on the important techniques used to minimize χ^{2} to get good theoretical parameters and estimate errors in the same. We will come back to this later in the notes.

• Mathews and Walker solve the simplest problem, i.e., ξ_{i} are linear in theoretical parameters.

• We do an example of this below.

• Bevington or similar book should be read before tackling a real problem.

Correlated Observations in χ^{2}

• Note that if x_{i} are calculated from the same population, e.g.,

x_{1} = \frac{1}{n} m ; ∑ ; k=1 y_{k}
x_{2} = \frac{1}{n} m ; ∑ ; k=1 y_{k}^{2}, etc.

• Then they are not independent,

• but we can use the theorem proved at the end of Section 2 to show that x_{i} will (for large m) be Gaussian.

• However, there is a correlation matrix, i.e., L ∝ exp \{−\frac{1}{2}(x − ξ)^{T}M(x − ξ)\} as discussed in Solmitz.

• By diagonalizing M we can (in principle) reduce this to the usual χ^{2} method with x_{i} replaced by eigenvectors of M.

Example of Multiple Moments

• This idea is important in the method of moments when one calculates several moments from the same distribution of events.

• In fact, one usually drops (forgets) off diagonal terms in M and forms χ^{2} from non-independent random variables.

• This is done in high energy physics when fitting the moments ρ_{00}, ρ_{1−1}, ρ_{10} of a spin 1^{−} particle decay.

• Of course, in this case, the off-diagonal elements are never given by the experimentalist recording his or her data: so one is forced to set them zero.

Linear fits to Cross Section Data

Example of Use of χ^{2}—Overview

• Consider the example in Mathews and Walker, page 390.

σ(θ) = a_{0} + a_{1} cos θ + a_{2} cos^{2} θ + a_{3} cos^{3} θ

• The measurements x_{i}’s are cross sections σ (θ_{i}) at particular scattering angles θ_{i}, and the theoretical predictions ξ_{i} are linear in unknown parameters a_{m}.

ξ_{i} = n_{T} ; ∑ ; m=1 C_{im}a_{m}

χ^{2} Formulation of Linear Fit

for Cross Section Data—I

• Minimize:

N ; ∑ ; i=1\frac{(x_{i} − ξ_{i})^{2}}{σ_{i}^{2}}
x_{i} = const. N_{i}
σ_{i} = const. \sqrt{N_{i}}
= \frac{x_{i}}{\sqrt{N_{i}}}

where i’th bin (i.e., range of θ observed, which are centered at θ_{i}) has N_{i} events.

• Minimizing χ^{2} with respect to unknown a_{m} gives

N ; ∑ ; i=1 \frac{(x_{i} − ξ_{i}^{∗})}{σ_{i}^{2}} \frac{∂ξ_{i}^{∗}}{∂a_{m}} = 0 or
_{l}^{∑} M_{ml} a_{l}^{∗} = b_{m}

• The superscript ^{∗} denotes χ^{2} or maximum likelihood estimate.

Statistics for Physics  3: Estimation of Parameters

χ^{2} Formulation of Linear Fit for Cross Section Data—II

• The “χ^{2}” matrix M_{ml} is given by:

M_{ml} = _{i}^{∑} \frac{C_{im}C_{il}}{σ_{i}^{2}}

• where M_{ml} is a positive definite symmetric matrix

= _{i}^{∑} \frac{(∂ξ_{i}/∂a_{m})(∂ξ_{i}/∂a_{l})}{σ_{i}^{2}}

• The right hand side

b_{m} = _{i}^{∑} \frac{x_{i}C_{im}}{σ_{i}^{2}} = _{i}^{∑} \frac{x_{i}}{σ_{i}^{2}} \frac{∂ξ_{i}}{∂a_{m}}

• Remember:

• ξ_{i} = σ_{theoretical} (θ_{i}) = a_{0} + a_{1} cos θ_{i} + a_{2} cos^{2} θ_{i} + a_{3} cos^{3} θ_{i}

• C_{im} = cos^{m} θ_{i} for m = 1, 2, 3

Statistics for Physics  3: Estimation of Parameters

Solution and Errors from Linear χ^{2}

• We can write the result in matrix form a^{∗} = M^{−1}b

• We can find errors in either of the two ways

• a) regard a^{∗} as a function of x_{i} and calculate standard deviation from central limit theorem.

• b) Use the shape of likelihood, L = exp (−\frac{1}{2}χ^{2}).

• We have—as far as shape of L is concerned: L = exp (independent of a+ no terms linear in (a − a^{∗}) due to “maximum” condition −\frac{1}{2}(a − a^{∗})^{T}M(a − a^{∗})),

• i.e., χ^{2} matrix M is indeed error matrix for a.

Estimate of Errors to be Used

\underline{in χ^{2} Fits}

• One should also mention that there are two versions of χ^{2} depending on whether error σ_{i} taken from “experiment” or “estimated theoretically.”

• In the latter case, it will be a function of the parameters to be varied.

• For example, suppose x_{i} is number of events in some counting experiment.

• The σ_{i} (experimental) = \sqrt{x_{i}}

• but σ_{i} (theory) = \sqrt{ξ_{i}}.

• Now x_{i} is only Gaussian for large x_{i}, but then \sqrt{x_{i}} and \sqrt{ξ_{i}} are essentially same.

• If x_{i} = 0, then theoretical σ is to be preferred!—although χ^{2} method itself is in difficulties as x_{i} not Gaussian.

• It is technically best to go back to maximum likelihood and put in true Poisson distribution for x_{i}.

3.3 Counter Physics—Introduction

and Poisson Process—I

• For historical reasons, I am forced to take my examples from high energy physics.

• Take M happy smiling counters—all arranged in a row: Suppose they count random events that fall in their province.

• Then each of them has a parameter λ_{i} associated with them such that distribution of number n_{i} of events in counter is a Poisson distribution with parameter λ_{i}t (t is time).

Counter Physics—Introduction

and Poisson Process—II

• This follows from our theory of Birth processes. In high energy physics, λ_{i} is proportional to _{∫_{BIN}}dt\frac{dσ}{dt} where BIN is t-bite covered by counter (t the momentum transfer is a particular kinematic range—not time!).

• The analysis is a straightforward application of maximum likelihood method:

p_{i} = e^{−λ_{i}t}\frac{(λ_{i}t)^{n_{i}}}{n_{i}!}  (21)

• (n_{i} observed counts, λ_{i} theory: t—now time!—may be known or not).

L = M ; ∏ ; i=1 p_{i}  or
ln L = − M ; ∑ ; i=1 ln n_{i}! − t M ; ∑ ; i=1 λ_{i} + M ; ∑ ; i=1 n_{i} ln(λ_{i}t) (22)

Unnormalized Counter Experiment—I

• Suppose the time t—which includes beam intensity, etc.—is unknown.

• Then we use Equation (22) to determine it (regard t as a theoretical parameter)

\frac{∂(ln L)}{∂t} = 0  or

t = \frac{M ; ∑ ; i=1n_{i}}{M ; ∑ ; i=1λ_{i}}

= \frac{N}{Λ}

where N = total Poisson parameter summed over all counters, and Λ = M ; ∑ ; i=1λ_{i}.

• Now if we eliminate t, we get

ln L = − M ; ∑ ; i=1 ln n_{i}! − N + N ln N
+ M ; ∑ ; i=1 n_{i} ln (\frac{λ_{i}}{Λ})

Statistics for Physics  3: Estimation of Parameters

Unnormalized Counter Experiment—II

• or dropping terms independent of λ_{i} :

ln L = M ; ∑ ; i=1 n_{i} ln (\frac{λ_{i}}{Λ})
= N ; ∑ ; j=1 ln p̃(C_{j})  (23)

where C_{j} is counter in which j’th event arrived and p̃(k) = \frac{λ_{k}}{Λ} is relative probability of landing in k’th counter.

• Note p̃ is normalized M ; ∑ ; k=1p̃(k) = 1.

• This is the discrete maximum likelihood method: • we have N observations:

• each has only M possible outcomes

• and the L = N ; ∏ ; j=1p̃(C_{j}) is product of probabilities—normalized properly.

Unnormalized Counter

Experiment—III

• If λ_{i} are functions of some theoretical parameters, we find them by maximizing Equation (23).

• In particular, if all values of μ_{i} = \frac{λ_{i}}{Λ} are free parameters, we maximize (23) subject to constant

M ; ∑ ; i=1μ_{i} = 1:

• this gives of course:

μ_{i} = \frac{n_{i}}{N}

Normalized Counter Experiment

• Now, we suppose t is no longer an unknown theoretical parameter, but a known number.

ln L = const. + \underbrace{\{N ln(tΛ) − tΛ\}}_{(a)} + \underbrace{M ; ∑ ; i=1 n_{i} ln(λ/Λ)}_{(b)}

• (b) is the form we found before.

• (a) can be regarded as ln p_{N} where p_{N} is probability of getting N (which should now be regarded as an experimental observable) total counts in time t

• (p_{N} = e^{−tΛ}(tΛ)^{N} which is up to a constant, the usual Poisson distribution).

• This requires some care to apply previous theorems as probabilities must be normalized in proofs in Equation (14). However, theorems do hold.

• Note Orear calls this normalized case as the extended maximum likelihood method.

Gaussian Limit of Counter Physics

• Of course, if each n_{i} is large

p_{i} ∝ exp [\frac{−\{n_{i} − λ_{i}t\}^{2}}{λ_{i}t}]

and we get back χ^{2} method.

• Note: σ_{i} is derived as theoretical value of error! (replace λ_{i}t in denominator by n_{i} to get experimental value).

• So we only get something different by using the full maximum likelihood method, if some n_{i} are small and Poisson distribution is not in regime where the distribution is essentially Gaussian.

3.4 Bubble Chamber

or Spectrometer Physics—I

• This is like Equation (16) except that every event (not every counter) has its own probability p_{i}.

• Typically, the events are labelled by continuous parameter(s), and so p_{i} is a continuous probability density.

• So far, we haven’t discussed how to apply maximum likelihood method to a counting problem with continuous probability distribution

• (lifetime example was continuous but not a counting experiment).

Bubble Chamber

or Spectrometer Physics—II

• What we do is divide range of continuous parameters {such as (0 ≤ θ ≤ 2π) in (1 + p cos θ)/2π polarization distribution} into M small bins.

• M should be so large that

• (a) probability of two events in a given experiment falling in same bin is zero.

• (b) the theoretical p_{i} is constant over each bin.

• then we have M observations—each with one of two results:

• success—1 event in bin

• failure—0 events in bin

Statistics for Physics  3: Estimation of Parameters

Bubble/Spectrometer—Physics Individual Events with Continuous

Distributions—III

• Let N be total number of events: then likelihood is

L = _{successes}^{∏} p(success) _{failures}^{∏} p(failure)
= N ; ∏ ; i=1 exp(−p_{i} t dx_{i})(p_{i} dx_{i} t)
⊗ M−N ; ∏ ; j=1 exp(−p_{j} dx_{j} t)

• writing out L in terms of Poisson probabilities for 0 or 1 events.

• Now M >> N and M bins dx_{j} cover space. Thus

L = exp [−∫ p(x) dx t] N ; ∏ ; i=1 p_{i} dx_{i} t

Bubble/Spectrometer—Physics

Individual Events with Continuous

Distributions—IV

• Now, fortunately, arbitrary dx_{i}, i = 1 … N are independent of theoretical parameters and so can be dropped (at further cost to the significance of normalization of L).

• So our fundamental formula for L in a continuous distribution is

L = exp [− ∫ p(x) dx t] N ; ∏ ; i=1 p(x_{i}, α)t

Bubble/Spectrometer—Physics

Individual Events with Continuous

Distributions—V

• Note that just as in Section 3.3—if time t is unknown we can find it from

\frac{∂L}{∂t} = 0  which gives
t = \frac{N}{∫ p(x) dx}

• when we get, with respect to remaining parameters α,

L(α) ∝ N ; ∏ ; i=1 \{\frac{p(x_{i}, α)}{∫ p(x, α)dx}\}  (24)

—the product of normalized probability distributions.

• If t is known, the discussion is again similar to this case in Section 3.3.

Individual Events Become

Binned Data—I

• If we have so many events that I can divide them into K groups of n_{k} (k = 1 … K) events where p has the same value for each member of the group, then:

L = exp {− ∫ p dx t} K ; ∏ ; k=1(pt)^{n_{k}}
∝ exp {− K ; ∑ ; k=1 λ_{k}t} K ; ∏ ; k=1(λ_{k}t)^{n_{k}}  (25)

where

λ_{k} = ∫_{bin for k’th group} p dx,

and we return exactly to formulae in Section 3.3.

• Further, if n_{k} is large for each k, our Poisson distributions become Gaussian and we can use χ^{2} method to analyze data.

• This is the only case where χ^{2} is correct in the sense of giving the best answer.

Individual Events Become

Binned Data—II

• However, even if N is not so large that p is constant in each group, then one can still bin data and throw away exact values of x_{i} for each event—just record what bin it is in.

• Obviously, if p is not constant over each bin, this throws away information.

• However, Equation (25) is now exact for L and one chooses bins so that n_{k} is large and χ^{2} applicable.

• This is used in practice because χ^{2} is so much easier and cheaper to use.

• Often one can bin data and minimize χ^{2}: this gives values for theoretical parameters which can be used as initial values for expensive full likelihood fit which will then need only a few iterations.

• One can carry piggy back further—and use a super low grade method (e.g., method of moments—see Section 3.5) to initialize parameters before χ^{2} fit.

• This is discussed by Solmitz (“histogram method”) and Eadie, et al.

Statistics for Physics  3: Estimation of Parameters

Method of Moments—I

• This is a simple and, in general, non-optimal method.

• Even here—for more than one parameter—there are two versions • (a) ignore correlations (used in practice but wrong) and

• (b) include correlations (best of a bad lot!). Of course, simplicity can lead to huge savings in computer time, and so it should not be sneered at.

• The basic idea of method of moments is simple:

• Consider a set of observations x_{i}, then central limit theorem

implies

\frac{1}{N} N ; ∑ ; i=1 f(x_{i}) → μ

• where

μ(α) = ∫ p(x, α)f(x)dx

• Here, μ is μ(α) where α is some unknown parameter. So putting

μ(α) = \frac{1}{N} N ; ∑ ; i=1 f(x_{i})

• gives one displaymath for α with solution α = α^{∗}:

• α^{∗} is the moment estimate of α.

Method of Moments—II

• Note:

• (a) One gets a different moment method for each choice of f:

• from maximum likelihood theorem, all of these are asymptotically worse or perhaps just not better than maximum likelihood method.

• One must choose moment which is simple and has small expected error.

• (b) If we have to estimate n theoretical parameters α, then we just get n equations by taking n moments (i.e., n different values of f).

Method of Moments—III

• (c) Error in μ is theoretically

\frac{1}{\sqrt{N}} \sqrt{∫ f^{2}p(x, α)dx − μ_{0}^{2}}

• and experimentally

\frac{1}{\sqrt{N}} \sqrt{\frac{1}{N} N ; ∑ ; i=1 f^{2}(x_{i}) − (\frac{1}{N} N ; ∑ ; i=1 f(x_{i}))^{2}}

• using μ_{0} to be true value of mean and

• μ = \frac{1}{N} N ; ∑ ; i=1f(x_{i}) to be estimate.

• By central limit theorem, for large N, experimental estimate will give correct error.

• However, if all x_{i} are the same (always true if N = 1!), then experimental error = 0, which is simply wrong.

• The moral is: use theoretical estimate if you only have a few events.

Statistics for Physics  3: Estimation of Parameters

Example of Moment Method

• Consider

p(x) = \frac{1}{2}(1 + ax) − 1 ≤ x ≤ 1

• Use Moment Method for any power m of x

⟨x^{m}⟩ = \frac{a}{2} +1 ; ∫ ; −1x^{m+1}dx
= \frac{1}{m + 2}  : m odd

• So if M_{m} = moment method using ⟨x^{m}⟩ and let M_{m} produce estimate a_{m}^{∗}.

• Then

a_{m}^{∗} = (m + 2) N ; ∑ ; i=1 \frac{x_{i}^{m}}{N}

• The error^{2} in estimate a_{m}^{∗} is

σ^{2}|_{m} = \frac{(m + 2)^{2}}{N} \{ +1 ; ∫ ; −1x^{2m}(1 + ax)dx − \frac{a^{2}}{(m + 2)^{2}}\}
= \frac{1}{N} \{\frac{(m + 2)^{2}}{2m + 1} − a^{2}\}

• which is a minimum for m = 1 when

σ^{2}|_{1} = \frac{1}{N} (3 − a^{2})

CPS713 STAT G.C.Fox, Feb 96  149  gcf@npac.syr.edu

Comparison of Moment Method

with Maximum Likelihood—I

• Now, Maximum Likelihood Method gives result with error estimate given by

\frac{1}{[Nσ^{2}]}|_{m.l.m.} = +1 ; ∫ ; −1 \frac{dx}{p} (\frac{∂p}{∂a})^{2}

• and from our general theory, this general formula gives in our general case:

\frac{1}{[Nσ^{2}]}|_{m.l.m.} = \frac{1}{2} +1 ; ∫ ; −1dx\frac{x^{2}}{(1 + ax)}
= \frac{1}{2a^{2}} \{ln [\frac{1 + a}{1 − a}] − 2a\}
= \frac{1}{3} + \frac{a^{2}}{5} + \frac{a^{3}}{7} + …
= \frac{1}{3} when a = 0

CPS713 STAT G.C.Fox, Feb 96  150  gcf@npac.syr.edu

Comparison of Moment Method

with Maximum Likelihood—II

• The value of \frac{1}{Nσ^{2}}|_{m.l.m.} is equal to \frac{1}{Nσ^{2}}|_{1} value when a = 0

• but

\frac{1}{[Nσ^{2}]}|_{m.l.m.} = ∞ when a = ±1.

• So

σ^{2}|_{m.l.m.} = σ^{2}|_{1} a = 0
< σ^{2}|_{1} a ≠ 0.

• i.e., maximum likelihood is always better, or as good as, moment method—in agreement with general theorem.

Several Parameters in

Method of Moments—I

• As we have already mentioned, the case of several parameters is treated by taking several moments.

• Now our moment errors are described by a full second moment matrix for off diagonal terms come from ∫ (f(x) − ⟨f⟩) (g(x) − ⟨g⟩) p(x) dx which is nonzero in general.

• It is conventional to ignore correlations:

• The expansion in terms of Y_{m}^{l}(θ, φ) of dσ/dt or an angular decay distribution is typical example.

• The expansion functions Y_{m}^{l} are orthonormal, so

W (θ, φ) = _{lm}^{∑} a_{lm}Y_{lm}
a_{lm} ∝ ⟨Y_{lm}^{∗}⟩

Several Parameters in

Method of Moments—II

• but still a_{lm} have correlated errors as ⟨Y_{lm}Y_{l′m′}⟩ ≠ ⟨Y_{lm}⟩⟨Y_{l′m′}⟩.

• For mean is weighted with probability distribution in θ and φ, and Y_{lm} are only orthogonal with respect to unit weight.

• These correlations in high energy physics example of dσ/dt express fact that backward peak (<< forward) comes from correlated cancellations.

Remarks on Counting Errors

in Moments—I

• In formula

⟨f⟩ = \frac{N ; ∑ ; i=1f(x_{i})}{N}  (26)

• I used to worry if you should include counting error (\sqrt{N}) in N in error estimate which otherwise only involves the Standard Deviation of f.

• If N is a priori fixed, clearly there is no \sqrt{N} error but suppose the events used in Equation (26) are those that fall in a particular subregion (counter) of a bigger sample space (and f(x) is defined only over this subregion).

• Then the estimate of probability of falling into counter certainly involves \sqrt{N} error! What about the mean of a quantity within subregion—does this involve this error?

• The problem is: what is error in a moment: should it include counting error?

• The answer is no—one uses identical formulae in case where number of events prescribed and not subject to counting error.

Statistics for Physics  3: Estimation of Parameters

Remarks on Counting Errors

in Moments—II

• We show this in a model problem in two variables:

• x is variable on which moment taken and y is binned.

• Choose:

P (x, y) = p(y)\frac{[1 + α(y)x]}{2}
Take a y-bin y_{1} ≤ y ≤ y_{2}
set S(y) = 1 : y_{1} ≤ y ≤ y_{2}
= 0 : otherwise.

• Assume α(y) is constant at α over bin and then method of moments is

α = 3 \frac{N_{T} ; ∑ ; i=1f(x_{i}, y_{i})}{N_{T} ; ∑ ; i=1g(x_{i}, y_{i})} = \frac{3F}{G}

• where

f(x, y) = xS(y)
g(x, y) = S(y)
CPS713 STAT G.C.Fox, Feb 96  155  gcf@npac.syr.edu

Remarks on Counting Errors

in Moments—III

• There are N_{T} total events and N fall into bin y_{1} ≤ y ≤ y_{2}.

• Now we can apply F/G theorem from end of Section 2,

• and it is not hard to show that f and g separately have counting errors \sqrt{N},

• but they are correlated and correlation cancels counting errors, so that α has identical error formula to standard moment method.

• (The estimate of p(y) will involve counting error, of course).

4: MINIMIZATION

Statistics for Physics  4: Minimization

4 Minimization

4.1 Statement of Problem

• The basic problem is to minimize f(α_{1} … α_{n})—a function of n variables α_{1} … α_{n}.

• Equivalently, we must find the zero of an arbitrary function f′(x_{1} … x_{n}) of n variables.

• This arises in previous sections when we wish either to minimize χ^{2} or maximize the likelihood L (equivalently to minimize—ln L).

• The methods divide according to

(i) Can calulate f only (ii) Can calculate f, f′, and f′′.

• (i) is extremely time consuming (i.e., it takes a lot of function evaluations to converge to desired minimum).

• In all problems I have faced, it has been simple to calculate derivatives and so I will only discuss (ii).

4.2 Derivative Formalism for χ^{2}—I

• One great advantage of the χ^{2} or maximum likelihood problem is that to calculate the second derivative of the function to be minimized only requires the first derivative of the theoretical functions.

• First consider this in the case of χ^{2}:

χ^{2} = N ; ∑ ; i=1 \frac{(t_{i} − e_{i})^{2}}{σ_{i}^{2}}  (27)

• where e_{i} is experimental measurement, σ_{i} is error, and t_{i} is theory—a function of n parameters.

• We take, without loss of generality, n = 1 and one parameter α in the following:

\frac{∂χ^{2}}{∂α} = N ; ∑ ; i=1 2(t_{i} − e_{i})\frac{∂t_{i}}{∂α}/σ_{i}^{2}
\frac{∂^{2}χ^{2}}{∂α^{2}} = S_{1} + S_{2}
S_{1} = N ; ∑ ; i=1 2(t_{1} − e_{i})\frac{∂^{2}t_{i}}{∂α^{2}}/σ_{i}^{2}
S_{2} = 2 N ; ∑ ; i=1 \frac{(\frac{∂t_{i}}{∂α} )2}{σ_{i}^{2}}

Statistics for Physics  4: Minimization

Derivative Formalism for χ^{2}—II

• Now, S_{1} is O(\sqrt{N}) and asymptotically negligible.

• To show this, regard it as a random variable which is a function of the N random variables e_{1} … e_{N}.

• By the Central Limit Theorem, S_{1} has mean

μ = N ; ∑ ; i=1 μ_{i}

• and Standard Deviation σ.

σ^{2} = N ; ∑ ; i=1 ω_{i}^{2}

• where

\frac{2(t_{i} − e_{i})}{σ_{i}^{2}} \frac{∂^{2}t_{i}}{∂α^{2}}

has mean μ_{i} and Standard Deviation ω_{i}.

CPS713 STAT G.C.Fox, Feb 96  160  gcf@npac.syr.edu

Derivative Formalism for χ^{2}—III

• Now if t_{i} was evaluated with true theoretical parameters μ_{i} = 0 and ω_{i} is some number.

• So σ^{2} is O(N) and S_{1} is O(\sqrt{N}) as desired.

• Meanwhile, S_{2} is O(N) as it is sum of N positive numbers.

• Whence—near solution

\frac{∂^{2}χ^{2}}{∂α^{2}} = S_{2} = 2 N ; ∑ ; i=1 (\frac{∂t_{i}}{∂α})^{2} /σ_{i}^{2}

• As claimed: second derivatives of χ^{2} only require evaluation of first derivatives of theory (only true when N large).

Derivative Calculation for

Maximum Likelihood—I

• Let’s do the same thing for maximum likelihood method.

− ln L = − N ; ∑ ; i=1 ln q_{i}

• where q_{i} are probability of our i = 1 … N events.

• Assume q_{i} are normalized ∫ q(α)dx = 1 as in Equation (25).

• This normalization is crucial but even when it does not seem to be true as in our subtle modifications in Section 4.3 with the “generalized” maximum likelihood method, one can easily modify the analysis, and show it is still valid.

• We get in normalized q_{i} case,

\frac{−∂^{2} ln L}{∂α^{2}} = S_{1} + S_{2}

• with

S_{1} = N ; ∑ ; i=1 \frac{1}{q_{i}} \frac{∂^{2}q_{i}}{∂α^{2}}

S_{2} = N ; ∑ ; i=1 \frac{1}{q_{i}^{2}} (\frac{∂q_{i}}{∂α})^{2}

CPS713 STAT G.C.Fox, Feb 96  162  gcf@npac.syr.edu

Derivative Calculation for

Maximum Likelihood—II

• Now, S_{1,2} are directly of form suitable to apply central limit theorem.

• Thus, they are sums over independent random variables taken from same distribution.

⟨S_{1}⟩ = N ⟨\frac{1}{q} \frac{∂^{2}q}{∂α^{2}} ⟩
= N ∫ q · \frac{1}{q} \frac{∂^{2}q}{∂α^{2}}dα
= N \frac{∂^{2}}{∂α^{2}} \{∫ q dα\}
= N \frac{∂^{2}}{∂α^{2}}(1) = 0

• So again S_{1} only has standard deviation term and is O(\sqrt{N}).

• Similarly S_{2} is O(N) and dominates.

• It follows that exactly same methods can be applied to χ^{2} and maximum likelihood fits.

• In the future we will only mention the former.

Statistics for Physics  4: Minimization

Derivative Calculation for

Maximum Likelihood—III

• As far as your computer program is concerned, one can handle χ^{2} and likelihood with same program.

• They just differ in calculational complexity N. • In χ^{2}, N is the number of “bins.”

• In likelihood, N is number of events.

• My old fitting programs • Manxcat—Marquardt’s Method

• “Feefoefum’s Well Nigh Perfect Fit” using eigenvalue method

• see later for methods

• were designed to interface with either χ^{2} or maximum likelihood.

4.3 The Derivative Method for χ^{2}—I

• Let’s first generalize the previous discussion to n parameters—getting a vector F_{j} of first and matrix M_{jk} of second derivatives

F_{j} = \frac{∂χ^{2}}{∂α_{j}} = N ; ∑ ; i=1 \frac{(t_{i} − e_{i})}{σ_{i}^{2}} \frac{∂t_{i}}{∂α_{j}}
M_{jk} = \frac{1}{2} \frac{∂^{2}χ^{2}}{∂α_{j}∂α_{k}} = N ; ∑ ; j,k=1 \frac{1}{σ_{i}^{2}} \frac{∂t_{i}}{∂α_{j}} \frac{∂t_{i}}{∂α_{k}} (28)

• Note that—in our approximation—M is positive semi-definite—which is also a technical advantage.

The Derivative Method for χ^{2}—II

• Now given any guess α = α_{0}, expand χ^{2}.

χ^{2}(α) = χ^{2}(α_{0})+_{j}^{∑}(α−α_{0})_{j} F_{j} +_{jk}^{∑}(α−α_{0})_{j} (α−α_{0})_{k} M_{jk}

• Adopt a vector notation with Δα_{j} = α_{j} − α_{0j}.

χ^{2}(\underline{α}) = χ^{2}(\underline{α}_{0}) + Δ\underline{α F} + Δ\underline{α}^{T}MΔ\underline{α}  (29)

• Minimum condition is

\frac{∂χ^{2}(\underline{α})}{∂α_{j}} = 0  or
MΔ\underline{α} = −\underline{F}  (30)

• So we have new guess

\underline{α} = Δ\underline{α} + \underline{α}_{0} = \underline{α}_{0} − M^{−1}\underline{F}

for true minimum.

• Note that as M positive definite in Equation (30), this new value of \underline{α} indeed gives a minimum (not a maximum or saddle point), as long as second order Taylor expansion valid.

Error Calculation from χ^{2} Expansion

• When we are at a minimum, \underline{F} = 0 and Equation (29) becomes

χ^{2}(\underline{α}) = χ^{2}(\underline{α}_{0}) + Δ\underline{α}^{T}MΔ\underline{α}

• or remembering that the likelihood

L = exp (−\frac{1}{2}χ^{2})
L = L(α_{0}) ⊗ exp (−\frac{1}{2}Δ\underline{α}^{T}MΔ\underline{α})  (31)

• Equation (31) is expected Gaussian form for L.

• According to Section 3, M is just inverse of error matrix, i.e.,

⟨(α^{∗} − α̃)_{j}(α^{∗} − α̃)_{k}⟩ = M_{jk}^{−1}

• where \underline{α}^{∗} = estimate from minimizing χ^{2},

• \underline{α̃} = true parameters.

• In particular, “error” σ_{j} in α_{j} is

σ_{j} = \sqrt{(M^{−1})_{jj}}  (32)

and the other elements in M^{−1} are the much forgotten correlation parameters.

• Note that Δα_{j} = σ_{j}, all other Δα_{j} = 0, produces a χ^{2} change of 1.

4.4 The Real World and χ^{2}

Minimization—I

• For a real problem, the iterative method described in Section 4.3 always diverges: one must invent various artistic devices to ensure procedure converges.

• The simplest reliable algorithm is called Marquardt’s method and this is intuitively easy to understand.

• Consider χ^{2} as a function of several parameters.

The Real World and χ^{2}

Minimization—II

• χ^{2} will have a bunch of hills and valleys and we are trying to discover the location of the lowest point in the valley.

• If our assumptions are correct, then χ^{2} is just a neat quadratic function, and we have found a step that exactly moves to the minimum.

• Unfortunately, this is not usually correct, as the step \underline{α}_{0} → \underline{α}_{0} + Δ\underline{α} maybe too large to support the assumption, that terms in \frac{∂^{k}(χ^{2})}{∂α^{k}}, k ≥ 3 are zero, to be correct.

The Real World and χ^{2} Minimization

Example of Two Variables—I

• Consider the case when there are two variables.

χ^{2} = χ_{0}^{2} + linear terms
+ A_{1} (α_{1} − α_{1}^{0})^{2} + A_{2} (α_{2} − α_{2}^{0})^{2}
+ B_{1} (α_{1} − α_{1}^{0})^{3} + B_{2} (α_{1} − α_{1}^{0})^{2} (α_{2} − α_{2}^{0}) + …
+ …

• where we have taken the special case where there is no term proportional to (α_{1} − α_{1}^{0})(α_{2} − α_{2}^{0});

• such a term only confuses the current discussion, and we will soon see how to put it back in.

• The difficulty (lack of convergence of iterative scheme) occurs when

A_{1} >> A_{2}

Statistics for Physics  4: Minimization

The Real World and χ^{2} Minimization Example of Two Variables—II Pictures of Two Variable Case

• so that as a function of α_{1}, we get

• and as a function of α_{2}, we get

with a valley with steep walls in α_{1} and gentle walls in α_{2} direction.

• The predicted shift from Equation (30), Δα_{2} ∼ \frac{1}{A_{2}} and is typically large compared to Δα_{1} ∼ \frac{1}{A_{1}}.

CPS713 STAT G.C.Fox, Feb 96  171  gcf@npac.syr.edu

The Real World and χ^{2} Minimization

Example of Two Variables—III

• Unfortunately, there is no reason to believe that this large shift is well estimated as assumptions that lead to expansion break down.

• Terms like B_{2}(Δα_{1})^{2}Δα_{2} are quite likely to be big.

• We see that estimate in Equation (30) is false and one will typically jump out of valley and find

χ^{2} (\underline{α}_{0} + Δ\underline{α}) >> χ^{2} (\underline{α}_{0})

• The method will diverge!

• In the above example, the estimate for Δα_{1} was valid and not spoiled by Δα_{2} estimate.

• In a real example, these will be mixing between α_{1} and α_{2} due to nonzero \frac{∂^{2}χ}{∂α_{1}∂α_{2}}.

• Both shifts Δα_{1}, Δα_{2} will be large and specious.

Statistics for Physics  4: Minimization

Eigenfunction Analysis

• There is an intuitively obvious solution of this.

• Consider the second order expansion term

(\underline{α} − \underline{α}_{0})^{T} M (\underline{α} − \underline{α}_{0})

• M is a symmetric positive semi-definite matrix.

• M can be diagonalized

• λ_{j} are eigenvalues of M.

• D is a diagonal matrix with diagonal entries λ_{j}.

• P is n × n matrix formed by eigenfunctions of M.

Transformation of M to Eigendirections

• Put \underline{β} = P\underline{α}.

• Then χ^{2} = χ_{0}^{2} + linear terms + λ_{1} (β_{1} − β_{01})^{2} + λ_{2} (β_{2} − β_{01})^{2} + . . ..

• Now the shifts Δβ_{j} can be calculated independently and do not affect each other.

• Assume that we order the eigenvalues

λ_{1} ≥ λ_{2} ≥ … ≥ λ_{n} ≥ 0

• Then clearly the first new shifts Δβ_{j} corresponding to λ_{1} … λ_{L}, say are reliable and the remainder j = L + 1 … n unreliable.

• Suppose we ignore the predictions for these last and set

Δβ_{j} = Predicted value j ≤ L
= 0 j > L

• We need a criterion for choosing L but subject to this, this method is intuitively reasonable.

What I Used to Do!

• In a program that I used for many years successfully, I chose L by the largest value such that

λ_{L} ≥ rλ_{1} and λ_{L+1} < rλ_{1}

• where initially the ratio r = 0.1 and r was increased or decreased based on experience,

• i.e., if the resultant χ^{2} (\underline{β}_{0} + Δβ) was indeed smaller (as predicted) or larger than χ^{2} (\underline{β}_{0}).

• This method is equivalent to replacing M by the regularized version.

Why is M Poorly Conditioned?

• I note that in using this method, I got to look at values of the λ_{j}. It was not uncommon for

\frac{λ_{max}}{λ_{min}} = \frac{λ_{1}}{λ_{n}} ∼ 10^{7} for n ∼ 30

• You can ask why one gets such large ratios.

• Intuitively, the eigenvectors of M with small eigenvalue correspond to linear combinations of α_{j} that are badly determined.

• So very small λ’s say that your theoretical formalism is redundant or alternatively the data inadequate.

• For instance, in fitting a cross section \frac{dσ}{d(cos θ)} to

n ; ∑ ; j=1α_{j}P_{j−1}(cos θ)

• then often the coefficients of the higher order Legendre functions are poorly determined, and

• if you choose n to be “too big,” then small eigenvalues will develop.

• In this case, one can cure the problem by reducing the top index n but in general no such simple solution exists.

CPS713 STAT G.C.Fox, Feb 96  176  gcf@npac.syr.edu

Marquardt’s Method

• Marquardt’s method is a simpler technique but with similar goals to the eigenvalue method.

• One replaces

M → M + QI

• where I is the identity matrix.

• So this means we replace M by

• Clearly, if λ_{L} ∼ Q, then the two formalism’s are roughly equivalent.

• Eigenvectors with eigenvalues > Q, are shifted by naive prediction while those with eigenvalue << Q, have their shift drastically reduced.

Estimate of Marquardt

Parameter Q

• One needs an artistic method to estimate Q which again is varied based on experience. Q is increased when predicted

χ^{2} (\underline{β}_{0} + Δβ) > χ^{2} (\underline{β}_{0})

and is otherwise decreased or left unchanged.

• The natural “scale” for Q is given by

\frac{1}{n} T r M = \frac{1}{n} n ; ∑ ; j=1 λ_{j}

the average eigenvalue.

• Marquardt’s method is less precise than the diagonalization method but much faster;

• finding eigenvalues and vectors is time consuming.

Experience with Marquardt’s Method

• I have used both successfully.

• I used Marquardt’s method in the analysis of some 10^{6} Fermilab events where on each event, I fitted known shower slopes to measured photon energy measurements in a segmented lead-scintillator sandwich.

• Here n ∼ 10 and the large number of fits (10^{6}) required a fast method.

• In these million fits, I saw no problem with divergence.

• This is case where it was necessary to find a fast reliable method.

• I used the eigenvalue method in analysis of experimental data in terms of theoretical models.

• Here, one is typically only doing a few (maybe one) fit and so speed of minimizer is irrelevant.

5: GOODNESS OF FIT

Statistics for Physics  5: Goodness of Fit

5. Goodness of Fit

5.1 The χ^{2} Distribution

• When the maximum likelihood method reduces to χ^{2}, it is relatively easy to judge goodness of fit.

χ^{2} = n ; ∑ ; i=1 \frac{(e_{i} − t_{i})^{2}}{σ_{i}^{2}}  (33)

where t_{i} depend on m theoretical parameters.

• Then the result of minimizing Equation (33) with respect to m theoretical parameters should be distributed according to χ^{2} distribution with n − m degrees of freedom.

• This function is plotted on p. 64 of EDJRS.

• If k is number of degrees of freedom.

p(χ^{2}) = \frac{(\frac{χ^{2}}{2} )^{\frac{k}{2}−1} exp (−\frac{χ^{2}}{2} )}{2Γ (\frac{k}{2} )}  (34)

• ⟨χ^{2}⟩ = k and Standard Deviation is \sqrt{2k}.

CPS713 STAT G.C.Fox, Feb 96  181  gcf@npac.syr.edu

Statistics for Physics  5: Goodness of Fit

χ^{2} Distribution for One Degree of Freedom

• Consider χ^{2} with one degree of freedom

χ^{2} = \frac{(x − t)^{2}}{σ^{2}} where
p(x) = \frac{1}{σ\sqrt{2π}} exp (−\frac{(x − t)^{2}}{2σ^{2}})
p(χ^{2}) = p(x)\frac{dx}{dχ^{2}} = \frac{p(x)σ^{2}}{2(x − t)} · 2
p(χ^{2}) = \frac{exp (−\frac{χ^{2}}{2} )}{\sqrt{2π} 2 \sqrt{χ^{2}}} · 2

• which is value in Equation (34) using Γ(1/2) = \sqrt{π} • (Here, factor of two comes as ±(x − t) given same χ^{2} value).

General χ^{2} Distribution

• One can derive the χ^{2} distribution for k degrees of freedom in two ways.

• 1.) Mathews and Walker

• Direct calculation as for k = 1 noting that k ; ∑ ; i=1x_{i}^{2}

with x_{i} of mean 0, Standard Deviation 1 is r^{2} where r is radial distance in k-dimensional space.

• Then factors multiplying exp (−\frac{χ^{2}}{2} ) in Equation (34) are just volume element in k dimensions.

• 2.) Eadie, et al. show characteristic function of χ^{2} with one degree of freedom is

+∞ ; ∫ ; −∞exp(i t χ^{2}) p(χ^{2}) d(χ^{2}) = \frac{1}{\sqrt{1 − 2 i t}}

• Therefore, characteristic function for χ^{2} with k degrees of freedom is (1 − 2 i t)^{−\frac{k}{2}}

• and this can be inverse fourier transformed to give Equation (34).

Statistics for Physics  5: Goodness of Fit

χ^{2} Characteristic Function

• Introduce the Characteristic Function F (χ^{2}) by

F (χ^{2}) = χ^{2} ; ∫ ; 0p(χ^{2})d(χ^{2})
= \frac{1}{2Γ (\frac{k}{2} )} χ^{2} ; ∫ ; 0(\frac{u}{2})^{(^{\frac{k}{2}})−1} e^{\frac{−u}{2}}du

• F has probability distribution w(F)

w(F) = p(χ^{2})\frac{d(χ^{2})}{dF} = 1

so F is uniformly distributed between 0 and 1.

5.2 Confidence Level

• We prefer to define C(χ^{2}) = 1 − F (χ^{2}) which is called the confidence level

C(χ^{2}) = 1 is χ^{2} = 0 “good”
= 0 is χ^{2} = ∞ “bad” or “low confidence”

• Suppose you conduct a series of experiments which gives you a sequence of χ^{2} values, χ^{2}|_{1,2...}.

• An example in high energy physics would come from analyzing a set of events each containing particle tracks—each track would have an individual χ^{2}.

• A given experiment would generate many million such χ^{2} values!

Typical Confidence Level Distributions

• If one can do the “errors” right, the resultant distribution in C(χ^{2}) should be flat

• But if you have “background,” i.e., events that do not satisfy claimed hypothesis, then one will have an excess of events at large χ^{2} >> (n − m) or small C(χ^{2}).

• C(χ^{2}) will look like

• The plot of C(χ^{2}) is a good way of estimating background.

Real World Confidence Levels

• Unfortunately, one often gets C(χ^{2}) plots that look like

• with no clear “flat” piece.

• This maybe because your “signal” is very small.

• It could be because you got your errors wrong (!), e.g., maybe the true σ_{i} are all larger than ones used in calculating χ^{2}.

• In this case, it makes sense to scale χ^{2} (χ_{new}^{2} = const. χ_{old}^{2}) until we get a nice flat distribution for confidence level C near 1 and analyze as above.

Error Estimate from χ^{2} Value

• If we expand χ^{2} near its minimum with—for convenience—just one theoretical parameter, λ.

χ^{2} = χ_{0}^{2} + α (λ − λ_{0})^{2}  (35)

• where λ_{0} is our estimate of λ.

• Then probability distribution of λ is just

α exp (−\frac{χ^{2}}{2})

α exp (−\frac{α(λ − λ_{0})^{2}}{2})

i.e., \frac{1}{\sqrt{α}} is Standard Deviation of estimate.

• Returning to Equation (35), we find that if we put λ = λ_{0} ± \frac{1}{\sqrt{α}} (i.e., at its one standard deviation value),

• χ^{2} is increased by 1.

• This χ^{2} increase “rule” is common method of finding “error” in estimate λ even when the Taylor expansion equation is not a good approximation.

5.3 Bartlett’s Function—Overview

• The previous analysis can be generalized to Bartlett’s S function which is defined as

S(λ) = \frac{1}{\sqrt{I(λ)}} \frac{∂ ln L}{∂λ}

• I(λ) is the information ⟨(\frac{∂}{∂λ} ln L)^{2}⟩

I(λ) = N ∫ [\frac{∂}{∂λ} ln P (x, λ)]^{2} P (x, λ)dx

• where N = Number of events and ∫ [\frac{∂}{∂λ} ln P (x, λ)]^{2} P (x, λ)dx can be calculated from theoretical forms.

Bartlett’s Function for χ^{2}

• From Equation (35) we have

\frac{∂}{∂λ} ln L = −\frac{2α(λ − λ_{0})}{2} for L = exp (−\frac{χ^{2}}{2})

• and using our earlier asymptotic result

α = I(λ) we find
S(λ) = \sqrt{I(λ)}(λ_{0} − λ)
= \sqrt{α}(λ_{0} − λ)

• Therefore, λ = λ_{0} ± one. Standard Deviation is just S(λ) = ∓1.

• In fact, for any number of parameters n,

• S(λ) has mean 0 and variance 1, and so the rule S(λ) = ±1 (±l for l Standard Deviations) is an appropriate generalization of the rule that one increases χ^{2} by 1(l).

5.4 Hypothesis Testing—I

• We wish to test a specific hypothesis H_{0} (referred to as null hypothesis) against a set of possible other hypotheses.

• In goodness of fit analysis, we test H_{0} (= hypothesis that data described by claimed theory with parameters found from maximizing L) against “all” other hypotheses.

• If we are trying to find type of a particle then

• H_{0} might be particle is an electron,

• H_{1}, that is a pion and these could be only possibilities.

Hypothesis Testing—II

• We are given our Standard N observations (x_{1} … x_{N}) on which to make our decision.

• We do this by constructing yet another statistic Y which is some function of x_{1} … x_{N}.

• We assume that Y has distribution P_{Y}^{(i)} if hypothesis i is true (i = 0 or 1 if two hypotheses).

• Y could be vector valued but we will take scalar case here.

• In the best of all possible worlds, one can find a Y such that the values taken by it under two hypotheses are disjoint, i.e.,

∫ P_{Y}^{(0)}(y)P_{Y}^{(1)}(y)dy = 0  (36)

• Unfortunately, Equation (36) is not always attainable

• but one should always try to minimize the probability overlap integral.

Hypothesis Testing—III

• We define r as region in Y space where we will deem hypothesis (1) to be true and (0) to be false.

• In the remaining region of Y space, R − r, (R is total Y range) we will take (0) to be true, (1) to be false.

• This scheme is characterized by two constants:

• 1. α = level of significance or size of test

= ∫ ; r P_{Y}^{(0)}(y)dy.

• α is chance that (0) will be declared false when it is in fact true.

• 2. β = ∫ ; R−rP_{Y}^{(1)}(y)dy is the contamination, i.e.,

• β is chance that (0) will be declared true when it is in fact false.

• Clearly, it depends on the application for the balance between efficiency 1 − α and contamination β.

5.5 Student’s t Distribution

for Hypothesis Testing—I

• Suppose that the x_{i} are drawn from a normal (Gaussian) distribution of unknown mean and Standard Deviation.

H^{(0)} is hypothesis that mean = μ^{(0)}
H^{(1)} is hypothesis that mean = μ^{(1)}

• Define

t = \frac{(\overline{x} − μ^{(0)}) \sqrt{n}}{\overline{σ}} where
\overline{x} = \frac{1}{n} n ; ∑ ; i=1x_{i}
\overline{σ}^{2} = \frac{1}{n − 1} n ; ∑ ; i=1 (x_{i} − \overline{x})^{2}

• Under hypothesis (0), t follows a so-called t-distribution with n − 1 degrees of freedom. (f = n − 1).

• t would be Gaussian with unit standard deviation if \overline{σ} taken from theory—not from experiment.

Statistics for Physics  5: Goodness of Fit

Student’s t Distribution

for Hypothesis Testing—II

• The probability distribution is

p(t) = \frac{Γ (\frac{f+1}{2})}{\sqrt{nπ} Γ (\frac{f}{2} )} \frac{1}{(1 + \frac{t^{2}}{f} )^{(f+1)/2}}

• f = 1 is Cauchy distribution.

• As f → ∞, the t distribution becomes Gaussian with mean 0 and Standard Deviation 1.

• This can be expected for as f → ∞, using data to find \overline{σ}, should be accurate.

• If μ_{1} > μ_{0}, choose 0 < t_{a} < (μ_{1} − μ_{0}), and

• one could say that 0 is true if t < t_{a} and 1 is true if t > t_{a}.

CPS713 STAT G.C.Fox, Feb 96  195  gcf@npac.syr.edu

Student’s t Distribution

for Hypothesis Testing—III

• The corresponding efficiency/contamination are trivially calculable from cumulative probability functions of the student’s t distribution.

• As discussed on page 168 of Brandt, this method can be adapted (use average \overline{σ}) to test if two sets of measurements have the same mean

• (e.g., Set 1: Smokers, Set 2: Non-Smokers)

• Above, we show the regions with α (level of significance) and β (contamination), defined earlier.

5.6 Robust Estimation–Motivation–I

• See Andrews, et al. “Robust Estimates of Location,” Princeton University Press;

• “Robustness in Statistics,” edited by R. L. Laurer, G. N. Wilkinson, Academic Press (1979).

• In the maximum likelihood method we found the method that had the smallest possible Standard Deviation (asymptotically at least).

• However, as we emphasized in the study of non-parametric estimates, it is quite possible (and indeed probable) that systematic error from incorrect theoretical form of probability is large.

• In many cases, one knows reasonably well the analytic form of the “interesting” effect but what one observes has other totally unknown contributions.

• One must choose between parameterizing the unknown or ignoring known form of known effect.

Robust Estimation—Motivation—II

• Consider an example from physics where one observes mass distribution with a Breit Wigner for a resonance, say, plus “small” (usually large!) contributions to the probability distribution (e.g., background which must be estimated as polynomial/phase space) which must be parameterized in an ad hoc fashion.

• It is hard to evaluate distortion due to an incorrect parameterization of these unknown contributions.

• Obviously, the correct choice of method will depend on relative size of signal (with known functional form), and background.

Robust Estimation—Motivation—III

• There is thus a great deal of interest in the study of so called robust techniques which are insensitive to errors in our knowledge about the underlying distribution.

• Statisticians are very interested in such techniques but (apart from non-parametric methods we have already discussed) there has been little use/experience of them in the type of hard problems that physicists are want to find.

• Further, the derivation of robust methods is rather intuitive and not as mathematically precise as those of classical statistics.

Robust Determination of Mean—I

• The first book cited (Andrews) has an instructive study of a very simple problem.

• Given a univariate sample {x^{(i)} i = 1, N} from some distribution,

• Find an estimate of the mean μ of the underlying probability distribution.

• The natural method is to estimate μ from the arithmetic mean

μ_{Am}^{∗} = \frac{1}{N} N ; ∑ ; i=1 x^{(i)}

• As long as the underlying distribution has a finite Standard Deviation, this will have an error ∝ \frac{1}{\sqrt{N}} which → 0 as N → ∞.

Statistics for Physics  5: Goodness of Fit

Robust Determination of Mean—II

• Further, the method is unbiased for any N, i.e., ⟨μ_{Am}^{∗}⟩ = μ

• If the x^{(i)} are distributed according to either the • Gaussian

p(x) ∝ exp (−\frac{(x − μ)^{2}}{2σ^{2}})

• or exponential form

p(x) ∝ exp (−\frac{x}{μ})

• then we know from maximum likelihood method that the arithmetic mean is the best method.

• However, for other distributions, the arithmetic mean is unlikely to be best.

Robust Determination of Mean—III

• For a trivial example, suppose I measured y = x^{2} where x had an exponential distribution.

• Then the best estimate of the mean of y is obviously

y_{clever}^{∗} = (\frac{1}{N} N ; ∑ ; i=1 (y^{(i)})^{1/2})^{2}

and not

y_{Am}^{∗} = \frac{1}{N} N ; ∑ ; i=1 y^{(i)}

• The problem of estimating μ is hardest for distributions p(x) that have long tails.

• Then with a small probability, we will generate events from the tail, but these will be way off the means,

• and if included in calculation of mean, will distort results in a finite sample.

Robust Determination of Mean

Trimmed Mean

• One method that is more robust than the arithmetic mean, is the trimmed mean

μ_{trim}^{∗} = \frac{N(1−p) ; ∑ ; i=pN+1x^{(i)}}{N(1 − 2p)}

where pN events are left off each end.

• For the special case of distributions symmetric about their mean, the trimmed μ^{∗} will also have mean μ and error α\frac{1}{\sqrt{N}}.

• The variance of μ_{trim}^{∗} will be larger than μ_{Am}^{∗} for Gaussianly distributed x^{(i)}.

• But as you investigate distributions broader than Gaussians you will find the variance of μ_{trim}^{∗} becoming smaller than that of μ_{Am}^{∗}.

• A special case of μ_{trim}^{∗} is the median (i.e., middle value of x^{(i)} when they are numerically ordered) gotten with p = 0.5.

Likelihood Based Robust

Estimation of Means—I

• An interesting class of robust estimators are based on the modification of the maximum likelihood equations.

• Thus, let

z^{(i)} = \frac{(x^{(i)} − μ)}{σ}

and suppose that p(x, μ) is just a function of z. (True for Gaussian case!).

• The likelihood equation for μ is

\frac{∂}{∂μ} ln L = 0

or

N ; ∑ ; i=1 \frac{∂}{∂μ} ln p (\frac{(x^{(i)} − μ)}{σ}) = 0

CPS713 STAT G.C.Fox, Feb 96  204  gcf@npac.syr.edu

Likelihood Based Robust

Estimation of Means—II

• Now, as we are assuming p to be an even function of z, i.e., symmetric about the mean, we write the above as

N ; ∑ ; i=1 Ψ (\frac{(x^{(i)} − μ)}{σ}) = 0  (37)

• where Ψ is an odd function of its argument.

• For the case of Gaussian p,

• Ψ is just (x^{(i)} − μ)/σ^{2} and Equation (37) is

\frac{1}{σ^{2}} N ; ∑ ; i=1 (x^{(i)} − μ_{L}^{∗}) = 0

• which gives the arithmetic mean as the maximum likelihood estimate of the mean.

Statistics for Physics  5: Goodness of Fit

Particular Robust Estimates for Likelihood—Huber

• Now we invent a function Ψ which is less sensitive to x^{(i)}’s far from the mean than (x^{(i)} − μ).

• In Huber’s method, one sets

Ψ(z) = z for |z| ≤ z_{c}
and Ψ(z) = z_{c} for z > z_{c}
= −z_{c} for z > −z_{c}

Statistics for Physics  5: Goodness of Fit

Particular Robust Estimates for

Likelihood—Andrews

• Many forms for Ψ have been proposed.

• The most successful was that due to Andrews with

\begin{matrix}Ψ(z) = sin(z/2.1) & |z| < 2.1π \\ = 0 & |z| ≥ 2.1π \end{matrix}

• Note z is \frac{(x−μ)}{σ} where σ supplies natural scale for problem.

• We do not need a very good estimate for σ but it should be robust!

• A simple choice is σ = mean|x^{(i)} − μ| which formula is applied iteratively as value of μ is refined.

CPS713 STAT G.C.Fox, Feb 96  207  gcf@npac.syr.edu

Comparison of Methods—I

• A large number of tests are described in first reference to this section. Choose three distributions.

• Let N = Gaussian Distribution with σ = 1.

• C = Cauchy.

• M = Gaussian ÷ uniform random variable between 0 and 1.

• Here is a sample of their results for 10 observations

• The above table gives variances in units of which I am not certain; but the relative values within a column are independent of this.

Comparison of Methods—II

• We see very clearly that the arithmetic mean is hopeless for broad distributions.

• The median and Andrews methods are very good in these cases with the Andrews being better in Gaussian case (σ only 4% more than best estimate) but less robust in other examples.

• Note that all the methods except the arithmetic mean require all the x^{(i)} being in memory so that we can order them.

• So the computer complexity is much higher for the robust estimates.

• The book of Andrews gives computer programs for all the robust estimators.

Statistics for Physics  6: Generation of Random Numbers

6: GENERATION OF RANDOM NUMBERS 6 Overview of Generation of

Random Numbers

• Knuth is a reasonable pedagogical description although it is a little out of date.

• Chiu (C3P-526) describes newer algorithms which are important in some cases.

• Press, et al. (Numerical Recipes) implements algorithms at about the level of Knuth.

• The National High Performance Software Exchange has a review of parallel random number generators.

Generation of Random

Numbers—Genuinely Random

• One needs specialized hardware to generate true random numbers,

• e.g.,

for successive decays (Λ → pπ^{−}) of a bunch of Λ’s.

• There are a lot of difficulties here.

• Not only must you measure θ in an unbiased fashion (impossible),

• but also to find a random number with 32 bits of precision would require several separate θ values.

• Maybe each gives four bits of information and so eight measurements of θ are needed to get one random number.

Congruential Pseudorandom Numbers

• Rather, one usually employs a deterministic method which gives “essentially” random numbers. These are called pseudorandom.

• Most “standard” library random number generators use the so called linear congruential generators which generate a sequence of integers I_{j} by

I_{j+1} = [aI_{j} + c] mod(m)

• Here, m is called the modulus, a is called the multiplier, and c is called the increment.

• This is deterministic—given any seed I_{0}, one generates an identical sequence I_{0}, I_{1}, I_{2} …

• There are most m possible random numbers, but there may be less. One could generate a cycle if

I_{j_{1}} = I_{j_{2}} j_{1} > j_{2}

• Then I_{j_{1}+1} = I_{j_{2}+1}, etc. Also, I_{j_{1}−1} = I_{j_{2}−1} and in particular I_{j_{1}−j_{2}} = I_{0}.

Congruential Random Numbers

in Monte Carlo

• If the cycle has a length m, then imagine our calculation of π

π = +1 ; ∫ ; −1 +1 ; ∫ ; −1dx dy θ [1 − x^{2} − y^{2}]

• Then random numbers after the m’th are no longer independent and so our use of central limit theorem breaks down.

• Our error is ∝ max [\frac{1}{\sqrt{N}} , \frac{1}{\sqrt{m}}] where N = number of events and m is cycle length.

• The references give various formal criteria to guide the choice of a, c, and m.

• Unfortunately, I do not see any of these references as being very good as they do not test and compare some key choices.

Statistics for Physics  6: Generation of Random Numbers

A Practical Congruential

Pseudorandom Number Generator

• In my group, we have used in the past

m = 2^{31}
a = 1103515245
c = 12345

• with C code which is essentially

• newran = [a * oldran + c]& MASK,

• where MASK has lower 31 bits = 1. All higher bits = 0.

• We get a floating point number between 0 and 1 by:

floatran = newran/2147483648.0

• This is very simple to code, especially in C which allows you to use machine capabilities effectively.

Shift Register Pseudorandom

Numbers—Overview

• The shift register methods are not discussed (correctly) in Knuth, and may well be better, as they are as easy to calculate and have longer period.

• Chiu has implemented a particularly simple method.

I_{j} = I_{j−q} ⊕ I_{j−p}

where ⊕ = Exclusive OR

1 ⊕ 1 = 0
1 ⊕ 0 = 0 ⊕ 1 = 1
0 ⊕ 0 = 0

This ⊕ operates independently on each bit of the numbers I_{j−q}, I_{j−p}.

• The period of this could be as large as 2^{p} − 1, which is big for Chiu’s choice p = 250, q = 147.

Initialization of

Pseudorandom Numbers

• The shift register method is unpopular because it is newer and because it is harder to initialize. However, Chiu describes a reasonable way using a congruential generator to get first random numbers.

• Typically, a random number package includes a call to a routine with a name like ranset.

• A call ranset (seed)

• sets I_{0} = seed for the congruential method.

• If you want a truly random start, a standard “dodge” is to set the seed to the binary image (remember I_{0} is just a (32 bit) integer) of the system clock.

• One may wish to redo a calculation with the same initial seed so use

clock (seed) print ( in binary form, seed) ranset (seed)

Correlations

• Often poorly chosen random number generators exhibit correlations between neighboring random numbers. This is well illustrated by tests by Chiu of integer generator supplied by Borland with Turbo Pascal V. 3-0.

• Other generators are better!

Decorrelation of

Pseudorandom Numbers

• Numerical Recipes has a program that removes any such problems. Its algorithm is essentially as follows:

• Set up a list R(1 … L) of L random numbers.

• Then choose randomly j where 1 ≤ j ≤ L.

• Set desired random number = R(j).

• Replace R(j) by a new random number.

• These correlations are particularly important in accurate calculations of condensed matter phase transitions

• as you are measuring correlations between local spins and so these must have uncorrelated random generators.

Generation of Gaussian

Random Numbers—I

• Gaussian Random Numbers can be generated by:

g = \frac{N ; ∑ ; i=1u_{i} − \frac{N}{2}}{σ\sqrt{N}}

• where σ^{2} = 1/12 and u_{i} are uniformly distributed between 0 and 1, is Gaussianly distributed for large N.

• This, for N ∼ 10 is sometimes used in computers.

• More convenient is often

g_{1} = (−2 ln u_{1})^{\frac{1}{2}} cos 2πu_{2}
g_{2} = (−2 ln u_{1})^{\frac{1}{2}} sin 2πu_{2}

where u_{i} are uniform in [0,1] and then g_{i} are independent random numbers which are Gaussianly distributed with zero mean and unit variance.

Statistics for Physics  6: Generation of Random Numbers

Generation of Gaussian Random Numbers—II

• You can prove this by taking polar coordinates in

exp [−\frac{1}{2} (g_{1}^{2} + g_{2}^{2})] dg_{1} dg_{2}

• with g_{1} = r cos θ, g_{2} = r sin θ.

• Then

exp (−\frac{1}{2} [g_{1}^{2} + g_{2}^{2}]) dg_{1} dg_{2}

= exp (−\frac{1}{2}r^{2}) r dr dθ
= exp (−\frac{1}{2}r^{2}) d (\frac{1}{2}r^{2}) dθ
= d (exp (−\frac{1}{2}r^{2})) dθ

• Therefore, the choice

u_{1} = exp (−\frac{1}{2}r^{2}) u_{2} = θ/2π

as uniformly distributed in [0,1], implies that g_{1}, g_{2} are Gaussian.

CPS713 STAT G.C.Fox, Feb 96  221  gcf@npac.syr.edu

Generation of General

Random Number Distributions

• Example: Exponential Distribution

• Suppose x has distribution e^{−x} for x = 0 to ∞.

e^{−x}dx = d(e^{−x})
y = e^{−x} is uniform

• Therefore, x = − ln y where y is a uniform random variable in range 0 to 1.

• More generally put,

p(x)dx = d(F (x))

• Then w(x) = F (x), the characteristic function is uniformly distributed

w = x ; ∫ ; −∞ p(x)dx

• (could be different lower limit if as above p(x) = 0 where x ≤ 0).

• You must be able to invert integral to be able to do this.

Generation of Exponentially

Distributed Random Numbers

• For the choice

p(x) = e^{−x}

• we get

w(x) = 1 − e^{−x}

• This is essentially what we got on the previous foil because

• e^{−x} (our intuitive choice) and w(x) = 1 − e^{−x} (our general case) are both uniform in [0,1] together.

Statistics for Physics  6: Generation of Random Numbers

Interpolation to Find General Distributions—I

• Suppose p(x) is given by a table

• in say equally spaced intervals x_{k} − x_{k−1} = δx.

Statistics for Physics  6: Generation of Random Numbers

Interpolation to Find General Distributions—II

• Then we can form the cumulative distribution

x_{i} 0 = F (x_{1})
x_{1} + δx p(x_{1})  = F (x_{1} + δx)
x_{1} + 2δx p(x_{1}) + p(x_{2}) = F (x_{1} + 2δx)

⋮

x_{m+1} m ; ∑ ; i=1p(x_{i}) = 1 if normalized = F (x_{m+1})

• Then generate a random number r uniformly with 0 ≤ r ≤ 1

• Then use inverse interpolation to find x F (x) = r

• Then x is distributed according to p(x)

Statistics for Physics  6: Generation of Random Numbers

Accept/Reject Method for General

Probability Distributions

• Let p_{max} = \begin{matrix}max \\ x \end{matrix} [p(x)] for x_{1} ≤ x ≤ x_{2}.

• Let w_{max} be any number > p_{max}.

• Generate x uniformly in [x_{1}, x_{2}].

• Generate r uniformly in [0, 1].

• If r > \frac{p(x)}{w_{max}}, ignore x and start again.

• If r ≤ \frac{p(x)}{w_{max}}, accept x.

• One can combine this with technique described on previous foils.

• Suppose f(x)dx = d(F (x)) and we generate x in [x_{1}, x_{2}] distributed according to f(x) using exact method described first.

• This assumes we choose f(x) so that it can be integrated to find F (x).

• Now use to accept-reject method or interpolation method with the probability distribution p′(x) = p(x)/f(x)

• We choose f(x) so that it is roughly p(x) so that p′(x) is more slowly varying than p(x).

Accept-Reject Methods

for Event Generation—I

• Accept-Reject is sometimes an important technique in event generation.

• Suppose we generate a set N of events i with weight w_{i} representing the production probability.

• Then if we want to calculate some simple quantity q_{i}, the canonical method is to take

⟨Q⟩ = \frac{N ; ∑ ; i=1q_{i}w_{i}}{N ; ∑ ; i=1w_{i}}  (38)

CPS713 STAT G.C.Fox, Feb 96  227  gcf@npac.syr.edu

Accept-Reject Methods

for Event Generation—II

• An Alternative is to use acceptance-rejection to find a set M of events j with “weight 1.” Then

⟨Q⟩ = \frac{1}{M} M ; ∑ ; j=1 q_{j}  (39)

• When is Equation (39) superior to Equation (38)?

• If you are just calculating ⟨Q⟩, then (38) contains all the information you have and is “best.”

• However, if you are doing a lot of work to calculate q_{i},

• or maybe there are many different quantities q_{i},

• the accept-reject technique is best.

• Thus, it is clearly a waste to spend a lot of time calculating q_{i} for an event whose weight w_{i} is small, e.g., 0.001 w_{max}!