Topics in Statistical Data Analysis for HEP
Lecture 2: Multivariate Methods

CERN-JINR European School

of High Energy Physics

Bautzen, 14–27 June 2009

Glen Cowan

Physics Department

Royal Holloway, University of London

g.cowan@rhul.ac.uk

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

Outline

Lecture #1: An introduction to Bayesian statistical methods

Role of probability in data analysis (Frequentist, Bayesian)

A simple fitting problem : Frequentist vs. Bayesian solution

Bayesian computation, Markov Chain Monte Carlo

Setting limits / making a discovery

Lecture #2: Multivariate methods for HEP

Event selection as a statistical test

Neyman-Pearson lemma and likelihood ratio test

Some multivariate classifiers:

Boosted Decision Trees

Support Vector Machines

Resources on multivariate methods

C.M. Bishop, Pattern Recognition and Machine Learning, Springer, 2006

T. Hastie, R. Tibshirani, J. Friedman, The Elements of Statistical Learning, Springer, 2001

R. Duda, P. Hart, D. Stork, Pattern Classification, 2nd ed., Wiley, 2001

A. Webb, Statistical Pattern Recognition, 2nd ed., Wiley, 2002

Materials from some recent meetings:

PHYSTAT conference series (2002, 2003, 2005, 2007,...) see	www.phystat.org

Caltech workshop on multivariate analysis, 11 February, 2008

indico.cern.ch/conferenceDisplay.py?confId=27385

SLAC Lectures on Machine Learning by Ilya Narsky (2006)

www-group.slac.stanford.edu/sluo/Lectures/Stat2006_Lectures.html

Books:

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

Linear decision boundaries

A linear decision boundary is only

optimal when both classes follow

multivariate Gaussians with equal

covariances and different means.

x1

x2

For some other cases a linear

boundary is almost useless.

x1

x2

Nonlinear transformation of inputs

x1

x2

We can try to find a transformation,

so that the transformed “feature space” variables can be separated

better by a linear boundary:

Here, guess fixed

basis functions

(no free parameters)

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 +1 for points on one side of the hyperplane and -1 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 function 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 margin, and xi = 0 at or on the right side.

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 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

Summary

Boosted Decision Trees and Support Vector Machines are two examples

of relatively modern developments in Machine Learning that are only

recently attracting attention in HEP.

There are now many multivariate methods on the market and it is

difficult to make general statements about performance; this is often

very specific to the problem.

Expect advanced multivariate methods to have a major impact in areas

where one struggles for statistical significance, not in precision

measurements.

A simpler (e.g. “cut-based”) analysis may be considered more robust,

but e.g. a 5s signal from an SVM supported by 4s from cuts may win.

Fortunately tools to investigate these methods are now widely available.

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

Boosted decision tree example

First use of boosted decision trees in HEP was for particle identification

for the MiniBoone neutrino oscillation experiment.

H.J.Yang, B.P. Roe, J. Zhu, “Studies of Boosted Decision Trees for MiniBooNE Particle Identification”, Physics/0508045, Nucl. Instum. & Meth. A 555(2005) 370-385.

B.P. Roe, H.J. Yang, J. Zhu, Y. Liu, I. Stancu, G. McGregor, ”Boosted decision trees as an alternative to artificial neural networks for particle identification”, physics/0408124, NIMA 543 (2005) 577-584.

Particle i.d. in MiniBooNE

Search for nm to ne oscillations

required particle i.d. using

information from Cherenkov

detector.

Large number (~200) input

variables measured for each

event.

H.J. Yang, MiniBooNE PID, DNP06

MiniBoone boosted decision tree

Here performance stable after a few hundred trees

MiniBooNE Decision tree performance

~ eb / es

Comparison of boosting algorithms

A number of boosting algorithms on the market; differ in the

update rule for the tree weight.

AdaBoost study with linear classifier

J. Sochman, J. Matas,

cmp.felk.cvut.cz

Start with a problem for which a linear classifier is weak:

AdaBoost study with linear classifier

J. Sochman, J. Matas,

cmp.felk.cvut.cz

Imperfect pdf estimation

What if the approximation we use (e.g., parametric form, assumption

of variable independence, etc.) to estimate p(x) is wrong?

If we use poor estimates to construct the test variable

then the discrimination between the event classes will not be optimal.

But can this cause us e.g. to make a false discovery?

Even if the estimate of p(x) used in the discriminating variable are

imperfect, this will not affect the accuracy of the distributions f(y|H0),

f(y|H1); this only depends on the reliability of the training data.

Using the classifier output for discovery

y

f(y)

y

N(y)

Normalized to unity

Normalized to expected

number of events

excess?

signal

background

background

search

region

Discovery = number of events found in search region incompatible

with background-only hypothesis. Maximize the probability of this

happening by setting ycut for maximum s/√b (roughly true).

ycut

Controlling false discovery

So for a reliable discovery what we depend on is an accurate estimate

of the expected number of background events, and this accuracy only

depends on the quality of the training data; works for any function y(x).

But we do not blindly rely on the MC model for the training data for

background; we need to test it by comparing to real data in control

samples where no signal is expected.

The ability to perform these tests will depend on on the complexity of

the analysis methods.

2008-05-12T14:01:21