Statistical Methods for LHC Physics

Glen Cowan

RHUL Physics

www.pp.rhul.ac.uk/~cowan

Particle Physics Seminar

University of Glasgow

26 February, 2009

Outline

Quick overview of physics at the Large Hadron Collider (LHC)

New multivariate methods for event selection

Boosted Decision Trees

Support Vector Machines

Some thoughts on Bayesian methods

Statistical combination of search channels

Outlook for data analysis at the LHC

Data analysis at the LHC

The LHC experiments are expensive

~ $1010 (accelerator and experiments)

the competition is intense

(ATLAS vs. CMS) vs. Tevatron

and the stakes are high:

4 sigma effect

5 sigma effect

So there is a strong motivation to extract all possible information

from the data.

The Large Hadron Collider

Counter-rotating proton beams

in 27 km circumference ring

pp centre-of-mass energy 14 TeV

( 10 TeV for 2009/10)

Detectors at 4 pp collision points:

ATLAS

CMS

LHCb (b physics)

ALICE (heavy ion physics)

general purpose

The ATLAS detector

2100 physicists

37 countries

167 universities/labs

25 m diameter

46 m length

7000 tonnes

~108 electronic channels

The Standard Model of particle physics

Matter...

+ gauge bosons...

photon (g), W±, Z, gluon (g)

+ relativity + quantum mechanics + symmetries... = Standard Model

25 free parameters (masses, coupling strengths,...).

Includes Higgs boson (not yet seen).

Almost certainly incomplete (e.g. no gravity).

Agrees with all experimental observations so far.

Many candidate extensions to SM (supersymmetry, extra dimensions,...)

A simulated SUSY event in ATLAS

high pT

muons

high pT jets

of hadrons

missing transverse energy

p

p

Background events

This event from Standard

Model ttbar production also

has high pT jets and muons,

and some missing transverse

energy.

 can easily mimic a SUSY event.

LHC data

At LHC, ~109 pp collision events per second, mostly uninteresting

do quick sifting, record ~200 events/sec

single event ~ 1 Mbyte

1 “year”  107 s, 1016 pp collisions / year

2  109 events recorded / year (~2 Pbyte / year)

For new/rare processes, rates at LHC can be vanishingly small

e.g. Higgs bosons detectable per year could be ~103

 'needle in a haystack'

For Standard Model and (many) non-SM processes we can generate

simulated data with Monte Carlo programs (including simulation

of the detector).

A simulated event

PYTHIA Monte Carlo

pp  gluino-gluino

.

.

.

Multivariate event selection

Suppose for each event we measure a set of numbers

x1 = jet pT

x2 = missing energy

x3 = particle i.d. measure, ...

follows some n-dimensional joint probability density, which

depends on the type of event produced, i.e., was it

We want to separate (classify) the

event types in a way that exploits the

information carried in many variables.

E.g. hypotheses (class labels) H0, H1, ...

Often simply “signal”, “background”

Finding an optimal decision boundary

Maybe select events with “cuts”:

xi < ci

xj < cj

Or maybe use some other type of decision boundary:

Goal of multivariate analysis is to do this in an “optimal” way.

H0

H0

H0

H1

H1

H1

Test statistics

The decision boundary is a surface in the n-dimensional space of

input variables, e.g.,

We can treat the y(x) as a scalar test statistic or discriminating

function, and try to define this function so that its distribution has the

maximum possible separation between the event types:

The decision boundary

is now effectively a single

cut on y(x), dividing

x-space into two

regions:

R0 (accept H0)

R1 (reject H0)

Constructing a test statistic

The Neyman-Pearson lemma states: to obtain the highest background

rejection for a given signal efficiency (highest power for a given

significance level), choose the acceptance region for signal such that

where c is a constant that determines the signal efficiency.

Equivalently, the optimal discriminating function is given by the

likelihood ratio:

N.B. any monotonic function of this is just as good.

Neyman-Pearson doesn't always help

The problem is that we usually don't have explicit formulae for the pdfs p(x|s), p(x|b), so for a given x we can't evaluate the likelihood ratio.

Instead we have Monte Carlo models for signal and background

processes, so we can produce simulated data:

generate

generate

“training data”

events of known type

Naive try: enter each (s,b) event into an n-dimensional histogram,

use e.g. M bins for each of the n dimensions, total of Mn cells.

n is potentially large  prohibitively large number of cells to populate,

can't generate enough training data.

General considerations

In all multivariate analyses we must consider e.g.

Choice of variables to use

Functional form of decision boundary (type of classifier)

Computational issues

Trade-off between sensitivity and complexity

Trade-off between statistical and systematic uncertainty

Our choices can depend on goals of the analysis, e.g.,

Event selection for further study

Searches for new event types

Decision boundary flexibility

overtraining

boundary too rigid

good trade-off

The decision boundary will be defined by some free parameters that

we adjust using training data (of known type) to achieve the best

separation between the event types.

Goal is to determine the boundary using a finite amount of training data

so as to best separate between the event types for an unseen data sample.

Some “standard” multivariate methods

Place cuts on individual variables

Simple, intuitive, in general not optimal

Linear discriminant (e.g. Fisher)

Simple, optimal if the event types are Gaussian distributed with

equal covariance, otherwise not optimal.

Probability Density Estimation based methods

Try to estimate p(x|s), p(x|b) then use

In principle best, difficult to estimate p(x) for high dimension.

Neural networks

Can produce arbitrary decision boundary (in principle optimal),

but can be difficult to train, result non-intuitive.

Decision trees

In a decision tree repeated cuts are made on a single variable until some stop criterion is reached.

The decision as to which variable is used is

based on best achieved improvement

in signal purity:

Example by MiniBooNE experiment,

B. Roe et al., NIM 543 (2005) 577

where wi. is the weight of the ith event.

Iterate until stop criterion reached, based e.g. on purity and minimum number of events in a node.

Decision trees (2)

The terminal nodes (leaves) are classified as signal or background

depending on majority vote (or e.g. signal fraction greater than a

specified threshold).

This classifies every point in input-variable space as either signal

or background, a decision tree classifier, with the discriminant function

Decision trees tend to be very sensitive to statistical fluctuations in

the training sample.

Methods such as boosting can be used to stabilize the tree.

Boosting

Boosting is a general method of creating a set of classifiers

which can be combined to achieve a new classifier that is more stable

and has a smaller error than any individual one.

Often applied to decision trees but, can be applied to any classifier.

Suppose we have a training sample T consisting of N events with

x1,...., xN event data vectors (each x multivariate)

y1,...., yN true class labels, +1 for signal, -1 for background

w1,...., wN event weights

Now define a rule to create from this an ensemble of training samples

T1, T2, .... , derive a classifier from each and average them.

AdaBoost

A successful boosting algorithm is AdaBoost (Freund & Schapire, 1997).

First initialize the training sample T1 using the original

x1,...., xN event data vectors

y1,...., yN true class labels (+1 or -1)

w1(1),...., wN(1) event weights

with the weights equal and normalized such that

Train the classifier f1(x) (e.g. a decision tree) using the weights w(1)

so as to minimize the classification error rate,

where I(X) = 1 if X is true and is zero otherwise.

Updating the event weights (AdaBoost)

Define the training sample for step k+1 from that of k by updating

the event weights according to

i = event index

k = training sample index

Assign a score to the kth classifier based on its error rate:

Normalize so that

Iterate K times, final classifier is

BDT example from MiniBooNE

~200 input variables for each event (n interaction producing e, m or p).

Each individual tree is relatively weak, with a misclassification

error rate ~ 0.4 – 0.45

B. Roe et al., NIM 543 (2005) 577

Monitoring overtraining

From MiniBooNE

example

Boosted decision tree summary

Advantage of boosted decision tree is it can handle a large number

of inputs. Those that provide little/no separation are rarely used as tree splitters are effectively ignored.

Easy to deal with inputs of mixed types (real, integer, categorical...).

If a tree has only a few leaves it is easy to visualize (but rarely use only a single tree).

There are a number of boosting algorithms, which differ primarily in the

rule for updating the weights (e-Boost, LogitBoost,...)

Other ways of combining weaker classifiers: Bagging (Boostrap-

Aggregating), generates the ensemble of classifiers by random sampling

with replacement from the full training sample.

Support Vector Machines

Support Vector Machines (SVMs) are an example of a kernel-based

classifier, which exploits a nonlinear mapping of the input variables

onto a higher dimensional feature space.

The SVM finds a linear decision boundary in the higher dimensional space.

But thanks to the “kernel trick” one does not every have to write down

explicitly the feature space transformation.

Some references for kernel methods and SVMs:

The books mentioned in www.pp.rhul.ac.uk/~cowan/mainz_lectures.html

C. Burges, A Tutorial on Support Vector Machines for Pattern Recognition,

research.microsoft.com/~cburges/papers/SVMTutorial.pdf

N. Cristianini and J.Shawe-Taylor. An Introduction to Support Vector Machines

and other kernel-based learning methods. Cambridge University Press, 2000.

The TMVA manual (!)

Linear SVMs

Consider a training data set consisting of

x1,...., xN event data vectors

y1,...., yN true class labels (+1 or -1)

Suppose the classes can be separated by a hyperplane defined by

a normal vector w and scalar offset b (the “bias”). We have

for all i

for all yi = +1

for all yi = -1

or equivalently

Bishop Ch. 7

Margin and support vectors

The distance between the hyperplanes defined by y(x)  x·w + b = +1

and y(x) = -1 is called the margin, which is:

If the training data are perfectly separated then this means there are

no points inside the margin.

Suppose there are points on the margin (this is equivalent to defining

the scale of w). These points are called support vectors.

Linear SVM classifier

The best classifier should have a large margin, so to maximize

we can minimize

subject to the constraints

for all i

We can define the classifier using

which is >0 for points on one side of the hyperplane and <0 on the other.

Lagrangian formulation

This constrained minimization problem can be reformulated using

a Lagrangian

positive Lagrange multipliers ai

We need to minimize L with respect to w and b and maximize

with respect to ai.

There is an ai for every training point. Those that lie on the margin

(the support vectors) have ai > 0, all others have ai = 0. The solution can be written

(sum only contains

support vectors)

Dual formulation

It can be shown that one finds the same solution a by minimizing

the dual Lagrangian

So this means that both the classifier function and the Lagrangian

only involve dot products of vectors in the input variable space.

The classifier output is thus

Nonseparable data

If the training data points cannot be separated by a hyperplane,

one can redefine the constraints by adding slack variables xi:

for all i

Thus the training point xi is allowed to be up to a distance xi on the wrong side of the boundary, and xi = 0 at or on the right side of the boundary.

For an error to occur we have xi > 1, so

is an upper bound on the number of training errors.

Cost function for nonseparable case

To limit the magnitudes of the xi we can define the error function that

we minimize to determine w to be

where C is a cost parameter we must choose that limits the amount

of misclassification. It turns out that for k=1 or 2 this is a quadratic

programming problem and furthermore for k=1 it corresponds to

minimizing the same dual Lagrangian

where the constraints on the ai become

Nonlinear SVM

So far we have only reformulated a way to determine a linear

classifier, which we know is useful only in limited circumstances.

But the important extension to nonlinear classifiers comes from first

transforming the input variables to feature space:

These will behave just as our new “input variables”. Everything

about the mathematical formulation of the SVM will look the same

as before except with f(x) appearing in the place of x.

Only dot products

Recall the SVM problem was formulated entirely in terms of dot

products of the input variables, e.g., the classifier output is

so in the feature space this becomes

The Kernel trick

How do the dot products help? It turns on that a broad class of

kernel functions can be written in the form:

Functions having this property must satisfy Mercer's condition

for any function g where

is finite.

So we don't even need to find explicitly the feature space transformation

f(x), we only need a kernel.

Finding kernels

There are a number of techniques for finding kernels, e.g., constructing

new ones from known ones according to certain rules (cf. Bishop Ch 6).

Frequently used kernels to construct classifiers are e.g.

polynomial

Gaussian

sigmoidal

Using an SVM

To use an SVM the user must as a minimum choose

a kernel function (e.g. Gaussian)

any free parameters in the kernel (e.g. the s of the Gaussian)

the cost parameter C (plays role of regularization parameter)

The training is relatively straightforward because, in contrast to neural

networks, the function to be minimized has a single global minimum.

Furthermore evaluating the classifier only requires that one retain

and sum over the support vectors, a relatively small number of points.

SVM in HEP

SVMs are very popular in the Machine Learning community but have

yet to find wide application in HEP. Here is an early example from

a CDF top quark anlaysis (A. Vaiciulis, contribution to PHYSTAT02).

signal

eff.

Software for multivariate analysis

TMVA, Höcker, Stelzer, Tegenfeldt, Voss, Voss, physics/0703039

StatPatternRecognition, I. Narsky, physics/0507143

Further info from www.hep.caltech.edu/~narsky/spr.html

Also wide variety of methods, many complementary to TMVA

Currently appears project no longer to be supported

From tmva.sourceforge.net, also distributed with ROOT

Variety of classifiers

Good manual

Bayesian vs. frequentist methods

Two schools of statistics use different interpretations of probability:

I. Relative frequency (frequentist statistics):

II. Subjective probability (Bayesian statistics):

In particle physics frequency interpretation most used, but subjective

probability can be more natural for non-repeatable phenomena:

systematic uncertainties, probability that Higgs boson exists...

Statistical vs. systematic errors

Statistical errors:

How much would the result fluctuate upon repetition	of the measurement?

Implies some set of assumptions to define	probability of outcome of the measurement.

Systematic errors:

What is the uncertainty in my result due to

uncertainty in my assumptions, e.g.,

model (theoretical) uncertainty;

modelling of measurement apparatus.

The sources of error do not vary upon repetition of the	measurement. Often result from uncertain

value of, e.g., calibration constants, efficiencies, etc.



Although numerical values of answer here same as in frequentist

case, interpretation is different (sometimes unimportant?)

Example: posterior pdf from MCMC

Sample the posterior pdf from previous example with MCMC:

Summarize pdf of parameter of

interest with, e.g., mean, median,

standard deviation, etc.

Statistical combination of Higgs channels in ATLAS

Status of Higgs search recently published in ATLAS "CSC Book"

Expected Performance of the ATLAS Experiment: Detector,

Trigger and Physics, arXiv:0901.0512, CERN-OPEN-2008-20.

Several statistical methods under study for combination of channels.

Here show method using “profile likelihood”

For now, concentrate on combining "discovery modes":

H → 

H → WW→ e

H → ZZ(*) → 4l (l = e, )

H → → ll, lh

Some modes not yet included, e.g., WW → ee, , lqq;

→ sensitivity will improve

Higgs production at the LHC

Statistical combination of channels

Treat systematics by means of profile likelihood method.

Considers fixed-mass hypothesis; "look-elsewhere-effect" must

be studied separately.

NB  and WW analyses also study floating-mass fits.

Approximate methods for discovery/exclusion significance.

Valid for L > 2 fb. For lower L need MC methods.

Other developments ongoing (Bayesian, CLs, look-elsewhere effect,...)

Statistical model

Bin i of a given channel has ni events, expectation value is

Expected signal and background are:

 is global strength parameter, common to all channels.

= 0 means background only, = 1 is SM hypothesis.

btot, s, b are

nuisance parameters

The likelihood function

 represents all nuisance parameters, e.g., background rate, shapes

here signal rate is

only parameter

of interest

The single-channel likelihood function uses Poisson model

for events in signal and control histograms:

The full likelihood function is

There is a likelihood Li(,i) for each channel, i = 1, …, N.

data in signal histogram

data in control

histogram

Profile likelihood ratio

To test hypothesized value of , construct profile likelihood ratio:

Maximized L for given 

Maximized L

Equivalently use q =  2 ln ():

data agree well with hypothesized  → q small

data disagree with hypothesized  → q large

Distribution of qunder assumption of  related to chi-square

(Wilk's theorem, approximation valid for roughly L > 2 fb):

p-value / significance of hypothesized 

Test hypothesized  by giving

p-value, probability to see data

with ≤ compatibility with 

compared to data observed:

Equivalently use significance,

Z, defined as equivalent number

of sigmas for a Gaussian

fluctuation in one direction:

Sensitivity

Discovery:

Generate data under s+b (= 1) hypothesis;

Test hypothesis  = 0 → p-value → Z.

Exclusion:

Generate data under background-only (= 0) hypothesis;

Test hypothesis .

If = 1 has p-value < 0.05 exclude mH at 95% CL.

Estimate median significance by setting data equal to expectation

values (Asimov data) or by using MC.

For median, can combine significances of individual channels.

For significance of e.g. a real data set, need global fit.

ATLAS combined discovery significance

ATLAS combined 95% CL exclusion limits

Summary of ATLAS Higgs combination

Study shown only uses Higgs channels most relevant to discovery

some modes not yet included → sensitivity will improve

Systematic uncertainties treated using profile likelihood.

Other developments ongoing (Bayesian, CLs, look-elsewhere effect,...)

Due to approximations used, combined results valid for L > 2 fb

For L = 2fb for ATLAS:

expected discovery 5 or more in 143 < mH < 179 GeV,

expected upper limit on mH at 95% CL is 115 GeV

Discussions underway on ATLAS/CMS combinations.

Software for global fit and full combination under development as

part of the RooStats project (joint ATLAS/CMS/RooFit/Root)

available as part of ROOT 5.22

Outlook for data analysis at the LHC

Recent developments from Machine Learning provide some new tools

for event classification with a number of advantages over those methods

commonly used in HEP, e.g.,

Boosted decision trees

Support Vector Machines

Bayesian methods can allow for a more natural treatment of

non-repeatable phenomena, e.g., model uncertainties.

MCMC can be used to marginalize posterior probabilities

Active development of combination methods for searches

Profile likelihood, Bayesian, “LEP”-style, ...

Software for these methods now much more easily available, expect

rapid development as the LHC begins to produce real results.

Quotes I like

“If you believe in something

you don't understand, you suffer,...”

– Stevie Wonder

“Keep it simple.

As simple as possible.

Not any simpler.”

– A. Einstein

Extra slides

Probability – quick review

Frequentist (A = outcome of

repeatable observation):

Subjective (A = hypothesis):

Conditional probability:

Bayes' theorem:

MCMC basics: Metropolis-Hastings algorithm

Goal: given an n-dimensional pdf

generate a sequence of points

1) Start at some point

2) Generate

Proposal density

e.g. Gaussian centred

about

3) Form Hastings test ratio

4) Generate

5) If

else

move to proposed point

old point repeated

6) Iterate

Metropolis-Hastings (continued)

This rule produces a correlated sequence of points (note how

each new point depends on the previous one).

For our purposes this correlation is not fatal, but statistical

errors larger than naive

The proposal density can be (almost) anything, but choose

so as to minimize autocorrelation. Often take proposal

density symmetric:

Test ratio is (Metropolis-Hastings):

I.e. if the proposed step is to a point of higher , take it;

if not, only take the step with probability

If proposed step rejected, hop in place.

Metropolis-Hastings caveats

Actually one can only prove that the sequence of points follows

the desired pdf in the limit where it runs forever.

There may be a “burn-in” period where the sequence does

not initially follow

Unfortunately there are few useful theorems to tell us when the

sequence has converged.

Look at trace plots, autocorrelation.

Check result with different proposal density.

Try using a number of different starting points and look

at spread in results.

The "look-elsewhere effect"

Look for Higgs at many mH values -- probability of seeing a large

fluctuation for some mH increased.

Combined significance shown here relates to fixed mH.

False discovery prob enhanced by ~ mass region explored / m

For H→ and H→WW, studied by allowing mH to float in fit:

→ 

Expected discovery significance

Expected p-value of SM vs mH

Example from validation exercise: ZZ(*) → 4l

5 level

(One minus)

cumulative

distributions.

Band gives 68%

CL limits.

½2

½2

½2

½2

Distributions of q0 for 2, 10 fb from MC compared to ½2

~

post-

erior

2008-05-12T14:01:21