Cowan - Aachen14 - 1 - Statistical Methods For Particle Physics

Statistical Methods for Particle Physics
Lecture 1: probability, random variables, MC

www.pp.rhul.ac.uk/~cowan/stat_aachen.html
Graduierten-Kolleg
RWTH Aachen
10-14 February 2014
Glen Cowan
Physics Department
Royal Holloway, University of London
g.cowan@rhul.ac.uk
www.pp.rhul.ac.uk/~cowan
G. Cowan
Aachen 2014 / Statistics for Particle Physics, Lecture 1
Outline
1 Probability
Definition, Bayes theorem, probability densities
and their properties, catalogue of pdfs, Monte Carlo
2 Statistical tests
general concepts, test statistics, multivariate methods,
goodness-of-fit tests
3 Parameter estimation
general concepts, maximum likelihood, variance of
estimators, least squares
4 Hypothesis tests for discovery and exclusion
discovery significance, sensitivity, setting limits
5 Further topics
systematic errors, Bayesian methods, MCMC
G. Cowan
Some statistics books, papers, etc.

G. Cowan, Statistical Data Analysis, Clarendon, Oxford, 1998
R.J. Barlow, Statistics: A Guide to the Use of Statistical Methods in
the Physical Sciences, Wiley, 1989
Ilya Narsky and Frank C. Porter, Statistical Analysis Techniques in
Particle Physics, Wiley, 2014.
L. Lyons, Statistics for Nuclear and Particle Physics, CUP, 1986
F. James., Statistical and Computational Methods in Experimental
Physics, 2nd ed., World Scientific, 2006
S. Brandt, Statistical and Computational Methods in Data
Analysis, Springer, New York, 1998 (with program library on CD)
J. Beringer et al. (Particle Data Group), Review of Particle Physics,
Phys. Rev. D86, 010001 (2012) ; see also pdg.lbl.gov sections on
probability, statistics, Monte Carlo
G. Cowan
Data analysis in particle physics

Observe events of a certain type
Measure characteristics of each event (particle momenta,

number of muons, energy of jets,...)
Theories (e.g. SM) predict distributions of these properties
up to free parameters, e.g., , GF, MZ, s, mH, ...
Some tasks of data analysis:
Estimate (measure) the parameters;
Quantify the uncertainty of the parameter estimates;
Test the extent to which the predictions of a theory
are in agreement with the data.
G. Cowan
Dealing with uncertainty

In particle physics there are various elements of uncertainty:
final theory not known,
thats why we search further
theory is not deterministic,
quantum mechanics
random measurement errors,
present even without quantum effects
things we could know in principle but dont,
e.g. from limitations of cost, time, ...
We can quantify the uncertainty using PROBABILITY
G. Cowan
A definition of probability
Consider a set S with subsets A, B, ...
Kolmogorov
axioms (1933)
From these axioms we can derive further properties, e.g.
G. Cowan
Conditional probability, independence

Also define conditional probability of A given B (with P(B) 0):
E.g. rolling dice:

Subsets A, B independent if:
If A, B independent,
N.B. do not confuse with disjoint subsets, i.e.,
G. Cowan
Interpretation of probability
I. Relative frequency
A, B, ... are outcomes of a repeatable experiment
cf. quantum mechanics, particle scattering, radioactive decay...

II. Subjective probability
A, B, ... are hypotheses (statements that are true or false)
Both interpretations consistent with Kolmogorov axioms.

In particle physics frequency interpretation often most useful,
but subjective probability can provide more natural treatment of
non-repeatable phenomena:
systematic uncertainties, probability that Higgs boson exists,...
G. Cowan
Bayes theorem
From the definition of conditional probability we have,
and
but
, so
Bayes theorem
First published (posthumously) by the

Reverend Thomas Bayes (17021761)
An essay towards solving a problem in the
doctrine of chances, Philos. Trans. R. Soc. 53
(1763) 370; reprinted in Biometrika, 45 (1958) 293.
G. Cowan
The law of total probability

Consider a subset B of
the sample space S,
divided into disjoint subsets Ai

such that i Ai = S,
Ai
B Ai
law of total probability
Bayes theorem becomes
G. Cowan
10
An example using Bayes theorem

Suppose the probability (for anyone) to have AIDS is:
prior probabilities, i.e.,
before any test carried out
Consider an AIDS test: result is + or -

probabilities to (in)correctly
identify an infected person
probabilities to (in)correctly
identify an uninfected person
Suppose your result is +. How worried should you be?

G. Cowan
11
Bayes theorem example (cont.)

The probability to have AIDS given a + result is
posterior probability
i.e. youre probably OK!
Your viewpoint: my degree of belief that I have AIDS is 3.2%
Your doctors viewpoint: 3.2% of people like this will have AIDS

G. Cowan
12
Frequentist Statistics general philosophy

In frequentist statistics, probabilities are associated only with
the data, i.e., outcomes of repeatable observations (shorthand:
).
Probability = limiting frequency

Probabilities such as
P (Higgs boson exists),
P (0.117 < s < 0.121),
etc. are either 0 or 1, but we dont know which.
The tools of frequentist statistics tell us what to expect, under
the assumption of certain probabilities, about hypothetical
repeated observations.
The preferred theories (models, hypotheses, ...) are those for
which our observations would be considered usual.
G. Cowan
13
Bayesian Statistics general philosophy

In Bayesian statistics, use subjective probability for hypotheses:
probability of the data assuming
hypothesis H (the likelihood)
posterior probability, i.e.,

after seeing the data
prior probability, i.e.,

before seeing the data
normalization involves sum

over all possible hypotheses
Bayes theorem has an if-then character: If your prior

probabilities were (H), then it says how these probabilities
should change in the light of the data.
No general prescription for priors (subjective!)
G. Cowan
14
Random variables and probability density functions

A random variable is a numerical characteristic assigned to an
element of the sample space; can be discrete or continuous.
Suppose outcome of experiment is continuous value x
f(x) = probability density function (pdf)
x must be somewhere
Or for discrete outcome xi with e.g. i = 1, 2, ... we have
probability mass function
x must take on one of its possible values
G. Cowan
15
Cumulative distribution function

Probability to have outcome less than or equal to x is
cumulative distribution function
Alternatively define pdf with

G. Cowan
16
Other types of probability densities

Outcome of experiment characterized by several values,
e.g. an n-component vector, (x1, ... xn)
joint pdf
Sometimes we want only pdf of some (or one) of the components
marginal pdf
x1, x2 independent if
Sometimes we want to consider some components as constant
conditional pdf
G. Cowan
17
Expectation values
Consider continuous r.v. x with pdf f (x).
Define expectation (mean) value as
Notation (often):
~ centre of gravity of pdf.
For a function y(x) with pdf g(y),

(equivalent)
Variance:
Notation:
Standard deviation:
~ width of pdf, same units as x.

G. Cowan
18
Covariance and correlation

Define covariance cov[x,y] (also use matrix notation Vxy) as
Correlation coefficient (dimensionless) defined as
If x, y, independent, i.e.,
, then
x and y, uncorrelated
N.B. converse not always true.

G. Cowan
19
Correlation (cont.)
G. Cowan
20
Error propagation
Suppose we measure a set of values
and we have the covariances
which quantify the measurement errors in the xi.
Now consider a function
What is the variance of
The hard way: use joint pdf
to find the pdf
then from g(y) find V[y] = E[y2] - (E[y])2.

Often not practical,
G. Cowan
may not even be fully known.
21
Error propagation (2)

Suppose we had
in practice only estimates given by the measured
Expand
to 1st order in a Taylor series about
To find V[y] we need E[y2] and E[y].

since
G. Cowan
22
Putting the ingredients together gives the variance of
G. Cowan
23

If the xi are uncorrelated, i.e.,
then this becomes
Similar for a set of m functions
or in matrix notation
G. Cowan
where
24

The error propagation formulae tell us the
covariances of a set of functions
in terms of
the covariances of the original variables.
Limitations: exact only if
linear.
Approximation breaks down if function
nonlinear over a region comparable
in size to the i.
y(x)
y
x
y(x)
?
N.B. We have said nothing about the exact pdf of the xi,
e.g., it doesnt have to be Gaussian.
G. Cowan
25
Error propagation special cases
That is, if the xi are uncorrelated:

add errors quadratically for the sum (or difference),
add relative errors quadratically for product (or ratio).
But correlations can change this completely...
G. Cowan
26
Error propagation special cases (2)

Consider
with
Now suppose = 1. Then
i.e. for 100% correlation, error in difference 0.
G. Cowan
27
Some distributions
Distribution/pdf
Binomial
Multinomial
Poisson
Uniform
Exponential
Gaussian
Chi-square
Cauchy
Landau
Beta
Gamma
Students t
G. Cowan
Example use in HEP

Branching ratio
Histogram with fixed N
Number of events found
Monte Carlo method
Decay time
Measurement error
Goodness-of-fit
Mass of resonance
Ionization energy loss
Prior pdf for efficiency
Sum of exponential variables
Resolution function with adjustable tails
28
Binomial distribution
Consider N independent experiments (Bernoulli trials):
outcome of each is success or failure,
probability of success on any given trial is p.
Define discrete r.v. n = number of successes (0 n N).
Probability of a specific outcome (in order), e.g. ssfsf is
But order not important; there are

ways (permutations) to get n successes in N trials, total
probability for n is sum of probabilities for each permutation.
G. Cowan
29
Binomial distribution (2)

The binomial distribution is therefore
random
variable
parameters
For the expectation value and variance we find:
G. Cowan
30
Binomial distribution (3)

Binomial distribution for several values of the parameters:
Example: observe N decays of W, the number n of which are

W is a binomial r.v., p = branching ratio.
G. Cowan
31
Multinomial distribution
Like binomial but now m outcomes instead of two, probabilities are
For N trials we want the probability to obtain:

n1 of outcome 1,
n2 of outcome 2,
nm of outcome m.
This is the multinomial distribution for
G. Cowan
32
Multinomial distribution (2)

Now consider outcome i as success, all others as failure.
all ni individually binomial with parameters N, pi
for all i
One can also find the covariance to be
Example:
represents a histogram
with m bins, N total entries, all entries independent.
G. Cowan
33
Poisson distribution
Consider binomial n in the limit
n follows the Poisson distribution:
Example: number of scattering events

n with cross section found for a fixed
integrated luminosity, with
G. Cowan
34
Uniform distribution
Consider a continuous r.v. x with - < x < . Uniform pdf is:
N.B. For any r.v. x with cumulative distribution F(x),

y = F(x) is uniform in [0,1].
Example: for 0 , E is uniform in [Emin, Emax], with
G. Cowan
35
Exponential distribution
The exponential pdf for the continuous r.v. x is defined by:
Example: proper decay time t of an unstable particle

( = mean lifetime)
Lack of memory (unique to exponential):
G. Cowan
36
Gaussian distribution
The Gaussian (normal) pdf for a continuous r.v. x is defined by:
(N.B. often , 2 denote

mean, variance of any
r.v., not only Gaussian.)
Special case: = 0, 2 = 1 (standard Gaussian):
If y ~ Gaussian with , 2, then x = (y - ) / follows (x).

G. Cowan
37
Gaussian pdf and the Central Limit Theorem

The Gaussian pdf is so useful because almost any random
variable that is a sum of a large number of small contributions
follows it. This follows from the Central Limit Theorem:
For n independent r.v.s xi with finite variances i2, otherwise
arbitrary pdfs, consider the sum
In the limit n , y is a Gaussian r.v. with
Measurement errors are often the sum of many contributions, so

frequently measured values can be treated as Gaussian r.v.s.
G. Cowan
38
Central Limit Theorem (2)

The CLT can be proved using characteristic functions (Fourier
transforms), see, e.g., SDA Chapter 10.
For finite n, the theorem is approximately valid to the
extent that the fluctuation of the sum is not dominated by
one (or few) terms.
Beware of measurement errors with non-Gaussian tails.
Good example: velocity component vx of air molecules.
OK example: total deflection due to multiple Coulomb scattering.
(Rare large angle deflections give non-Gaussian tail.)
Bad example: energy loss of charged particle traversing thin
gas layer. (Rare collisions make up large fraction of energy loss,
cf. Landau pdf.)
G. Cowan
39
Multivariate Gaussian distribution

Multivariate Gaussian pdf for the vector
are column vectors,
are transpose (row) vectors,
For n = 2 this is
where = cov[x1, x2]/(12) is the correlation coefficient.

G. Cowan
40
Chi-square (2) distribution

The chi-square pdf for the continuous r.v. z (z 0) is defined by
n = 1, 2, ... = number of degrees of

freedom (dof)
For independent Gaussian xi, i = 1, ..., n, means i, variances i2,

follows 2 pdf with n dof.
Example: goodness-of-fit test variable especially in conjunction
with method of least squares.
G. Cowan
41
Cauchy (Breit-Wigner) distribution

The Breit-Wigner pdf for the continuous r.v. x is defined by
( = 2, x0 = 0 is the Cauchy pdf.)

E[x] not well defined, V[x] .
x0 = mode (most probable value)
= full width at half maximum
Example: mass of resonance particle, e.g. , K*, 0, ...
= decay rate (inverse of mean lifetime)
G. Cowan
42
Landau distribution
For a charged particle with = v /c traversing a layer of matter
of thickness d, the energy loss follows the Landau pdf:
+ - + -

- + - +

L. Landau, J. Phys. USSR 8 (1944) 201; see also

W. Allison and J. Cobb, Ann. Rev. Nucl. Part. Sci. 30 (1980) 253.
G. Cowan
43
Landau distribution (2)

Long Landau tail
all moments
Mode (most probable

value) sensitive to ,
particle i.d.
G. Cowan
44
Beta distribution
Often used to represent pdf

of continuous r.v. nonzero only
between finite limits.
G. Cowan
45
Gamma distribution
Often used to represent pdf

of continuous r.v. nonzero only
in [0,].
Also e.g. sum of n exponential
r.v.s or time until nth event
in Poisson process ~ Gamma
G. Cowan
46
Student's t distribution
= number of degrees of freedom

(not necessarily integer)
= 1 gives Cauchy,
gives Gaussian.
G. Cowan
47
Student's t distribution (2)

If x ~ Gaussian with = 0, 2 = 1, and
z ~ 2 with n degrees of freedom, then
t = x / (z/n)1/2 follows Student's t with = n.
This arises in problems where one forms the ratio of a sample
mean to the sample standard deviation of Gaussian r.v.s.
The Student's t provides a bell-shaped pdf with adjustable
tails, ranging from those of a Gaussian, which fall off very
quickly, ( , but in fact already very Gauss-like for
= two dozen), to the very long-tailed Cauchy ( = 1).
Developed in 1908 by William Gosset, who worked under
the pseudonym "Student" for the Guinness Brewery.
G. Cowan
48
The Monte Carlo method

What it is: a numerical technique for calculating probabilities
and related quantities using sequences of random numbers.
The usual steps:
(1) Generate sequence r1, r2, ..., rm uniform in [0, 1].
(2) Use this to produce another sequence x1, x2, ..., xn
distributed according to some pdf f (x) in which
were interested (x can be a vector).
(3) Use the x values to estimate some property of f (x), e.g.,
fraction of x values with a < x < b gives
MC calculation = integration (at least formally)
MC generated values = simulated data
use for testing statistical procedures
G. Cowan
49
Random number generators

Goal: generate uniformly distributed values in [0, 1].
Toss coin for e.g. 32 bit number... (too tiring).
random number generator
= computer algorithm to generate r1, r2, ..., rn.
Example: multiplicative linear congruential generator (MLCG)
ni+1 = (a ni) mod m , where
ni = integer
a = multiplier
m = modulus
n0 = seed (initial value)
N.B. mod = modulus (remainder), e.g. 27 mod 5 = 2.
This rule produces a sequence of numbers n0, n1, ...
G. Cowan
50
Random number generators (2)

The sequence is (unfortunately) periodic!
Example (see Brandt Ch 4): a = 3, m = 7, n0 = 1
sequence repeats
Choose a, m to obtain long period (maximum = m - 1); m usually
close to the largest integer that can represented in the computer.
Only use a subset of a single period of the sequence.
G. Cowan
51
Random number generators (3)

are in [0, 1] but are they random?
Choose a, m so that the ri pass various tests of randomness:
uniform distribution in [0, 1],
all values independent (no correlations between pairs),
e.g. LEcuyer, Commun. ACM 31 (1988) 742 suggests
a = 40692
m = 2147483399
Far better generators available, e.g. TRandom3, based on Mersenne

twister algorithm, period = 219937 - 1 (a Mersenne prime).
See F. James, Comp. Phys. Comm. 60 (1990) 111; Brandt Ch. 4
G. Cowan
52
The transformation method

Given r1, r2,..., rn uniform in [0, 1], find x1, x2,..., xn
that follow f (x) by finding a suitable transformation x (r).
Require:
i.e.
That is,
G. Cowan
set
and solve for x (r).
53
Example of the transformation method

Exponential pdf:
Set
and solve for x (r).

works too.)
G. Cowan
54
The acceptance-rejection method
Enclose the pdf in a box:
(1) Generate a random number x, uniform in [xmin, xmax], i.e.

r1 is uniform in [0,1].
(2) Generate a 2nd independent random number u uniformly
distributed between 0 and fmax, i.e.
(3) If u < f (x), then accept x. If not, reject x and repeat.
G. Cowan
55
Example with acceptance-rejection method
If dot below curve, use

x value in histogram.
G. Cowan
56
Improving efficiency of the

acceptance-rejection method
The fraction of accepted points is equal to the fraction of
the boxs area under the curve.
For very peaked distributions, this may be very low and
thus the algorithm may be slow.
Improve by enclosing the pdf f(x) in a curve C h(x) that conforms
to f(x) more closely, where h(x) is a pdf from which we can
generate random values and C is a constant.
Generate points uniformly
over C h(x).
If point is below f(x),
accept x.
G. Cowan
57
Monte Carlo event generators

Simple example: e+e- +-

Generate cos and :
Less simple: event generators for a variety of reactions:

e+e- +-, hadrons, ...
pp hadrons, D-Y, SUSY,...
e.g. PYTHIA, HERWIG, ISAJET...
Output = events, i.e., for each event we get a list of
generated particles and their momentum vectors, types, etc.
G. Cowan
58
A simulated event
PYTHIA Monte Carlo

pp gluino-gluino
G. Cowan
59
Monte Carlo detector simulation

Takes as input the particle list and momenta from generator.
Simulates detector response:
multiple Coulomb scattering (generate scattering angle),
particle decays (generate lifetime),
ionization energy loss (generate ),
electromagnetic, hadronic showers,
production of signals, electronics response, ...
Output = simulated raw data input to reconstruction software:
track finding, fitting, etc.
Predict what you should see at detector level given a certain
hypothesis for generator level. Compare with the real data.
Estimate efficiencies = #events found / # events generated.
Programming package: GEANT
G. Cowan
60
Next time...
Today we have focused on probabilities, how they are
defined, interpreted, quantified, manipulated, etc.
In the following lecture we will begin talking about statistics,
i.e., how to make inferences about probabilities (e.g., probabilistic
models or hypotheses) given a sample of data.
G. Cowan
61
Extra slides
G. Cowan
62
Multivariate distributions
Outcome of experiment characterized by several values, e.g. an
n-component vector, (x1, ... xn)
joint pdf
Normalization:
G. Cowan
63
Marginal pdf
Sometimes we want only pdf of
some (or one) of the components:
marginal pdf
x1, x2 independent if
G. Cowan
64
Marginal pdf (2)
Marginal pdf ~
projection of joint pdf
onto individual axes.
G. Cowan
65
Conditional pdf
Sometimes we want to consider some components of joint pdf as
constant. Recall conditional probability:
conditional pdfs:
Bayes theorem becomes:

Recall A, B independent if
x, y independent if
G. Cowan
66
Conditional pdfs (2)

E.g. joint pdf f(x,y) used to find conditional pdfs h(y|x1), h(y|x2):
Basically treat some of the r.v.s as constant, then divide the joint
pdf by the marginal pdf of those variables being held constant so
that what is left has correct normalization, e.g.,
G. Cowan
67
Functions of a random variable

A function of a random variable is itself a random variable.
Suppose x follows a pdf f(x), consider a function a(x).
What is the pdf g(a)?
dS = region of x space for which

a is in [a, a+da].
For one-variable case with unique
inverse this is simply
G. Cowan
68
Functions without unique inverse

If inverse of a(x) not unique,
include all dx intervals in dS
which correspond to da:
Example:
G. Cowan
69
Functions of more than one r.v.

Consider r.v.s
and a function
dS = region of x-space between (hyper)surfaces defined by
G. Cowan
70
Functions of more than one r.v. (2)

Example: r.v.s x, y > 0 follow joint pdf f(x,y),
consider the function z = xy. What is g(z)?
(Mellin convolution)
G. Cowan
71
More on transformation of variables

Consider a random vector
with joint pdf
Form n linearly independent functions

for which the inverse functions
exist.
Then the joint pdf of the vector of functions is

where J is the
Jacobian determinant:
For e.g.
G. Cowan
integrate
over the unwanted components.
72

Cowan - Aachen14 - 1 - Statistical Methods For Particle Physics

Uploaded by

Copyright:

Available Formats

Cowan - Aachen14 - 1 - Statistical Methods For Particle Physics

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Cowan - Aachen14 - 1 - Statistical Methods For Particle Physics

Uploaded by

Copyright:

Available Formats

Statistical Methods for Particle Physics

Lecture 1: probability, random variables, MC

Aachen 2014 / Statistics for Particle Physics, Lecture 1

Aachen 2014 / Statistics for Particle Physics, Lecture 1

Some statistics books, papers, etc.

Aachen 2014 / Statistics for Particle Physics, Lecture 1

Data analysis in particle physics

Measure characteristics of each event (particle momenta,

Aachen 2014 / Statistics for Particle Physics, Lecture 1

Dealing with uncertainty

Aachen 2014 / Statistics for Particle Physics, Lecture 1

Aachen 2014 / Statistics for Particle Physics, Lecture 1

Conditional probability, independence

E.g. rolling dice:

Aachen 2014 / Statistics for Particle Physics, Lecture 1

cf. quantum mechanics, particle scattering, radioactive decay...

Both interpretations consistent with Kolmogorov axioms.

Aachen 2014 / Statistics for Particle Physics, Lecture 1

First published (posthumously) by the

Aachen 2014 / Statistics for Particle Physics, Lecture 1

The law of total probability

divided into disjoint subsets Ai

law of total probability

Bayes theorem becomes

Aachen 2014 / Statistics for Particle Physics, Lecture 1

An example using Bayes theorem

Aachen 2014 / Statistics for Particle Physics, Lecture 1

Bayes theorem example (cont.)

Aachen 2014 / Statistics for Particle Physics, Lecture 1

Frequentist Statistics general philosophy

Probability = limiting frequency

Aachen 2014 / Statistics for Particle Physics, Lecture 1

Bayesian Statistics general philosophy

posterior probability, i.e.,

prior probability, i.e.,

normalization involves sum

Bayes theorem has an if-then character: If your prior

Aachen 2014 / Statistics for Particle Physics, Lecture 1

Random variables and probability density functions

Aachen 2014 / Statistics for Particle Physics, Lecture 1

Cumulative distribution function

Alternatively define pdf with

Aachen 2014 / Statistics for Particle Physics, Lecture 1

Other types of probability densities

Aachen 2014 / Statistics for Particle Physics, Lecture 1

~ centre of gravity of pdf.

For a function y(x) with pdf g(y),

~ width of pdf, same units as x.

Aachen 2014 / Statistics for Particle Physics, Lecture 1

Covariance and correlation

Correlation coefficient (dimensionless) defined as

N.B. converse not always true.

Aachen 2014 / Statistics for Particle Physics, Lecture 1

Aachen 2014 / Statistics for Particle Physics, Lecture 1

to find the pdf

then from g(y) find V[y] = E[y2] - (E[y])2.

may not even be fully known.

Aachen 2014 / Statistics for Particle Physics, Lecture 1

Error propagation (2)