Is there evidence for a peak in this
data?

1

Is there evidence for a peak in thisdata?

“Observa6on of an Exo6c S=+1
Baryon in Exclusive Photoproduc6on from the Deuteron”
S. Stepanyan et al, CLAS Collab, Phys.Rev.LeK. 91 (2003) 252001“The sta6s6cal signiﬁcance of the peak is 5.2 ± 0.6 σ”

2Is there evidence for a peak in thisdata?

“Observa6on of an Exo6c S=+1
Baryon in Exclusive Photoproduc6on from the Deuteron”
S. Stepanyan et al, CLAS Collab, Phys.Rev.LeK. 91 (2003) 252001“The sta6s6cal signiﬁcance of the peak is 5.2 ± 0.6 σ”“A Bayesian analysis of pentaquark signals from CLAS data”
D. G. Ireland et al, CLAS Collab, Phys. Rev. LeK. 100, 052001 (2008)“The ln(RE) value for g2a (-0.408) indicates weak evidence in favourof the data model without a peak in the spectrum.”
Comment on “Bayesian Analysis of Pentaquark Signals from CLAS
3Data”
Bob Cousins, hKp://arxiv.org/abs/0807.1330

Sta6s6cal Issues in Searches for
New Physics
Louis Lyons
Imperial College, London
andOxfordKyoto
May 2016

4•PhyStat-v• INTERNATIONAL WORKSHOP ON STATISTICAL ISSUESIN NEUTRINO PHYSICS•THE KAVLI INSTITUTE FOR THE PHYSICS AND MATHEMATICS OF THE UNIVERSE, KASHIWA,
JAPAN

•

30 MAY – 1 JUNE 2016••

Local Organising Committee: Mark HARTZ/Christophe BRONNER/Richard
CALLAND/Yoshinari HAYATO/Yasuhiro NISHIMURA/Kimihiro OKUMURA

•
•
•
•

Scientific Organising Committee: Yoshi UCHIDA/Jun CAO/Daniel
CHERDACK/Robert COUSINS/David VAN DYK/Mark HARTZ/Pilar
HERNANDEZ/Joe FORMAGGIO/Thomas JUNK/Asher KABOTH/LouisLYONS/Shun SAITO/Subir SARKAR/Elizabeth WORCESTER/Kai ZUBER.
PhyStatNu_LOC@ipmu.jp conference.ipmu.jp/PhyStat-nu PhyStat-nu@imperial.ac.uk( Ask Google for PHYSTAT-nu )5

Theme:Using data to make judgements about H1 (New Physics) versus
H0 (S.M. with nothing new)Why?
Experiments are expensive and time-consumingsoWorth investing effort in statistical analysis
à better information from data
Topics:Blind Analysis
Why 5σ for discovery?
Significance
P(A|B) ≠ P(B|A)
Meaning of p-values
Wilks’ Theorem
LEE = Look Elsewhere EffectBackground Systematics
Coverage
Example of misleading inference
p0 v p1 plots
Higgs search: Discovery and spin
(N.B. Several of these topics have no unique solutions from Statisticians)
Conclusions

6

Choosing between 2 hypotheses
Hypothesis testing: New particle or statistical fluctuation?
H0 = b
H1 = b + s

7Choosing between 2 hypotheses
Possible methods:Δχ2p-value of statistic à
lnL–ratio
Bayesian:
Posterior odds
Bayes factor
Bayes information criterion (BIC)
Akaike ……..
(AIC)
Minimise “cost”
See ‘Comparing two hypotheses’
hKp://www.physics.ox.ac.uk/users/lyons/H0H1_A~1.pdf

8

(a)

(b)

H0

tobs
With 2 hypotheses,
each with own pdf,
p-values aredefined as tailareas, pointing in
towards each other

H1

tobs

t
p1

t
p0

(c)
H0

H1

tobs

t

9

Procedure for choosing between 2 hypotheses
1) No sensitivity
H02) Maybe

3) Easy separation

H1

t

β

tcrit αProcedure: Obtain expected distributions for data statistic (e.g. L-ratio) for H0 and H1
Choose α (e.g. 95%, 3σ, 5σ ?) and CL for p1 (e.g. 95%)
Given b, α determines tcrit
b+s defines β.

For s > smin, separation of curves à discovery or excln

1-β = Power of test
Now data:

If tobs ≥ tcrit (i.e. p0 ≤ α), discovery at level α
If tobs < tcrit, no discovery.

If p1 < 1– CL, exclude H1

10BLIND ANALYSES
Why blind analysis?

Data sta6s6c, selec6ons, correc6ons, method

Methods of blinding

Add random number to result *Study procedure with simula6on only
Look at only ﬁrst frac6on of data
Keep the signal box closed
Keep MC parameters hidden
Keep unknown frac6on visible for each bin

Disadvantages
Takes longer 6me
Usually not available for searches for unknown

Aner analysis is unblinded, don’t change anything unless ……..

* Luis Alvarez sugges6on re “discovery” of free quarks
11Why 5σ for Discovery?
Sta6s6cians ridicule our belief in extreme tails (esp. for systema6cs)
Our reasons:1) Past history (Many 3σ and 4σ eﬀects have gone away)
2) LEE (see later)
3) Worries about underes6mated systema6cs
4) Subconscious Bayes calcula6onp(H1|x) = p(x|H1) * π(H1)p(H0|x)p(x|H0) π(H0)Posterior Likelihood Priors
prob
ra6o
“Extraordinary claims require extraordinary evidence”
N.B. Points 2), 3) and 4) are experiment-dependent
Alterna6ve sugges6on:
L.L. “Discovering the signiﬁcance of 5σ”

hKp://arxiv.org/abs/1310.1284

12

How many σ’s for discovery?
SEARCHSURPRISE

IMPACTLEE

SYSTEMATICS No. σ

Higgs searchMediumVery highM

Medium5

Single top

No

Low

No

No

3

SUSY

Yes

Very highVery largeYes

7

Bs oscilla@ons

Medium/Low Medium

Δm

No

4

Neutrino osc

MediumHigh

sin22ϑ, Δm2

No4

Bsà μ μNo

Low/Medium

No

Medium3

Pentaquark

Yes

High/V. highM, decaymode

Medium7(g-2)μ anom

Yes

High

No

Yes

4

H spin ≠ 0

Yes

High

No

Medium5

4th gen q, l, ν

Yes

High

M, modeNo

6

Dark energy

Yes

Very highStrength

Yes

5

Grav Waves

No

High

EnormousYes

8

Sugges6ons to provoke discussion, rather than `delivered on Mt. Sinai’
Bob Cousins: “2 independent expts each with 3.5σ beKer than one expt with 5σ”

13SigniﬁcanceSignificance = S/√B

or similar ?

Potential Problems:•Uncertainty in B
•Non-Gaussian behaviour of Poisson, especially in tail
•Number of bins in histogram, no. of other histograms [LEE]•Choice of cuts, bins

(Blind analyses)

For future experiments:• Optimising: Could give S =0.1, B = 10-4, S/√B =10
14

P(A|B) ≠ P(B|A)
Remind Lab or University media contact person that:
Prob[data, given H0] is very small
does not imply that
Prob[H0, given data] is also very small.
e.g. Prob{data | speed of ν ≤ c}= very small
does not imply
Prob{speed of ν≤c | data} = very smallorProb{speed of ν>c | data} ~ 1
Everyday situation, 2nd most convincing example:
Pack of playing cardsp(spade|king) = 1/4
p(king|spade) = 1/13

15

P(A|B) ≠ P(B|A)
Remind Lab or University media contact person that:
Prob[data, given H0] is very small
does not imply that
Prob[H0, given data] is also very small.
e.g. Prob{data | speed of ν ≤ c}= very small
does not imply
Prob{speed of ν≤c | data} = very smallorProb{speed of ν>c | data} ~ 1
Everyday examplep(pregnant|female) ~ 3%p(female|pregnant) >> 3%

What p-values are (and are not)H0 pdfReject H0 if t > tcrit (p < α )p-value = prob that t ≥ tobs
tcrit
Small p à data and theory have poor compa6bility
Small p-value does NOT automa6cally imply that theory is unlikely
Bayes prob(Theory|data) related to prob(data|Theory) = Likelihoodby Bayes Th, including Bayesian priorp0 = αtp-values are misunderstood. e.g. An6-HEP jibe:“Par6cle Physicists don’t know what they are doing, because half theirp ˂ 0.05 exclusions turn out to be wrong”
Demonstrates lack of understanding of p-values[All results rejec6ng energy conserva6on with p ˂α =.05 cut will turn out tobe ‘wrong’]

17

Combining diﬀerent p-values
Several results quote independent p-values for same eﬀect:
p1, p2, p3…..
e.g. 0.9, 0.001, 0.3 ……..
What is combined signiﬁcance?
Not just p1*p2*p3…..
If 10 expts each have p ~ 0.5, product ~ 0.001 and is clearly NOT correct
combined p
n −1
S = z * (-ln z)j j! , z = p1p2p3…….
j =0(e.g. For 2 measurements, S = z * (1 - lnz) ≥ z )Problems:
1) Recipe is not unique (Uniform dist in n-D hypercube à uniform in 1-D)
2) Formula is not associa@ve
Combining {{p1 and p2}, and then p3} gives diﬀerent answer
from {{p3 and p2}, and then p1} , or all together
Due to diﬀerent op@ons for “more extreme than x1, x2, x3”.
3) Small p’s due to diﬀerent discrepancies∑/******* BeKer to combine data ************
18

Wilks’ Theorem
Data = some distribution e.g. mass histogramFor H0 and H1, calculate best fit weighted sum of squares S0 and S1Examples: 1) H0 = polynomial of degree 3
H1 = polynomial of degree 52) H0 = background only
H1 = bgd+peak with free M0 and cross-section
3) H0 = normal neutrino hierarchy
H1 = inverted hierarchy
If H0 true, S0 distributed as χ2 with ndf = ν0
If H1 true, S1 distributed as χ2 with ndf = ν1
If H0 true, what is distribution of ΔS = S0 – S1? Expect not large.
Wilks’ Theorem:
ΔS distributed as χ2 with ndf = ν0 – ν1 provided:a) H0 is trueb) H0 and H1 are nestedc) Params for H1à H0 are well defined, and not on boundaryd) Data is asymptotic

Is it χ2?

19

Wilks’ Theorem, contd
Examples: Does Wilks’ Th apply?
1) H0 = polynomial of degree 3
H1 = polynomial of degree 5
YES: ΔS distributed as χ2 with ndf = (d-4) – (d-6) = 2
2) H0 = background only
H1 = bgd + peak with free M0 and cross-section
NO: H0 and H1 nested, but M0 undefined when H1à H0. ΔS≠χ2(but not too serious for fixed M)
3) H0 = normal neutrino hierarchy*********
H1 = inverted hierarchy*********
NO: Not nested. ΔS≠χ2(e.g. can have Δχ2 negative)
N.B. 1: Even when W. Th. does not apply, it does not mean that ΔS
is irrelevant, but you cannot use W. Th. for its expected distribution.
N.B. 2: For large ndf, better to use ΔS, rather than S1 and S0 separatelyIs diﬀerence in S distributed as χ2 ?

What is peak at zero?
Why not half the entries?

Demortier:
H0 = quadratic bgd
H1 = ……………… +Gaussian of fixed width,variable location & ampl

Protassov, van Dyk, Connors, ….
H0 = continuum(a) H1 = narrow emission line(b) H1 = wider emission line(c) H1 = absorption line
Nominal significance level = 5%21Is difference in S distributed as χ2 ?, contd.

So need to determine the ΔS distribution by Monte Carlo
N.B.
1) For mass spectrum, determining ΔS for hypothesis H1
when data is generated according to H0 is not trivial,
because there will be lots of local minima
2) If we are interested in 5σ significance level, needs lots ofMC simulations (or intelligent MC generation)3) Asymptotic formulae may be useful (see K. Cranmer, G. Cowan,
E. Gross and O. Vitells, 'Asymptotic formulae for likelihood-based tests of new
physics',

22

Look Elsewhere Eﬀect (LEE)
Prob of bgd ﬂuctua6on at that place = local p-value
Prob of bgd ﬂuctua6on ‘anywhere’ = global p-value
Global p > Local p
Where is `anywhere’?
a) Any loca6on in this histogram in sensible rangeb) Any loca6on in this histogramc) Also in histogram produced with diﬀerent cuts, binning, etc.
d) Also in other plausible histograms for this analysise) Also in other searches in this PHYSICS group (e.g. SUSY at CMS)f) In any search in this experiment (e.g. CMS)g) In all CERN expts (e.g. LHC expts + NA62 + OPERA + ASACUSA + ….)h) In all HEP expts
etc.
d) relevant for graduate student doing analysisf) relevant for experiment’s Spokesperson
INFORMAL CONSENSUS:
Quote local p, and global p according to a) above.
Explain which global p

23

Example of LEE: Stonehenge

24

25

Are alignments signiﬁcant?
•Atkinson replied with his ar6cle "Moonshine on Stonehenge"
in An1quity in 1966, poin6ng out that some of the pits which ….. had used
for his sight lines were more likely to have been natural depressions, andthat he had allowed a margin of error of up to 2 degrees in his alignments.
Atkinson found that the probability of so many alignments being visible
from 165 points to be close to 0.5 rather that the "one in a million"
possibility which ….. had claimed.

•….. had been examining stone circles since the 1950s in search ofastronomical alignments and the megalithic yard. It was not un6l 1973
that he turned his aKen6on to Stonehenge. He chose to ignore alignments
between features within the monument, considering them to be too close
together to be reliable. He looked for landscape features that could havemarked lunar and solar events. However, one of …..'s key sites, Peter's
Mound, turned out to be a twen6eth-century rubbish dump.
26

Background systema6cs

27

Background systema6cs, contd
Signif from comparing χ2’s for H0 (bgd only) and for H1 (bgd + signal)
Typically, bgd = func6onal form fa with free params
e.g. 4th order polynomial
Uncertain6es in params included in signif calcula6on
But what if func6onal form is diﬀerent ? e.g. fb
Typical approach:If fb best ﬁt is bad, not relevant for systema6cs
If fb best ﬁt is ~comparable to fa ﬁt, include contribu6on to systema6cs
But what is ‘~comparable’?
Other approaches:
Proﬁle likelihood over diﬀerent bgd parametric forms
hKp://arxiv.org/pdf/1408.6865v1.pdf?
Background subtrac6on
sPlotsNon-parametric background
Bayes
etc
No common consensus yet among experiments on best approach{Spectra with mul6ple peaks are more diﬃcult}28

“Handling uncertain@es in backgroundshapes: the discrete proﬁling method”
Dauncey, Kenzie, Wardle and Davies (Imperial College, CMS)
arXiv:1408.6865v1 [physics.data-an]
Has been used in CMS analysis of Hàγγ
Problem with ‘Typical approach’: Alterna6ve func6onal
forms do or don’t contribute to systema6cs by hard cut, so
systema6cs can change discon6nuously wrt ∆χ2
Method is like proﬁle L for con6nuous nuisance params
Here ‘proﬁle’ over discrete func6onal forms
29Reminder of Proﬁle LυStat uncertainty on s from widthof L ﬁxed at υbestsContours of lnL(s,υ)
s = physics paramυ = nuisance param

Total uncertainty on s from widthof L(s,υprof(s)) = Lprof
υprof(s) is best value of υ at that s
υprof(s) as fn of s lies on green line
Total uncert ≥ stat uncertainty
30

-2lnL∆s
31

Red curve: Best value of nuisance param υBlue curves: Other values of υ
Horizontal line: Intersec6on with red curveà
sta6s6cal uncertainty‘Typical approach’: Decide which blue curves have small enough ∆Systema6c is largest change in minima wrt red curves’.
Proﬁle L: Envelope of lots of blue curves
Wider than red curve, because of systema6cs (υ)
For L = multi-D Gaussian, agrees with ‘Typical approach’
Dauncey et al use envelope of finite number of functional forms
32Point of controversy!
Two types of ‘other func6ons’:a) Diﬀerent func6on types e.g.
Σai xi versus Σai/xib) Given fn form but diﬀerent number of terms
DDKW deal with b) by -2lnL à -2lnL + kn
n = number of extra free params wrt best
k = 1, as in AIC (= Akaike Informa6on Criterion)
Opposi6on claim choice k=1 is arbitrary.
DDKW agree but have studied diﬀerent values, and say k=1 is op6mal for them.
Also, any parametric method needs to make such a choice33

Example of misleading inference
Ofer Vitells, Weizmann Ins6tute PhD thesis (2014)
On-oﬀ problem (signal + bgd, bgd only)
e.g. non = 10, moﬀ = 0
i.e. convincing evidence for signal
Now, to improve analysis, look at spectra of events (e.g. in mass) in “on” and“oﬀ” regions
e.g. Use 100 narrow bins à ni = 1 for 10 bins, mi = 0 for all bins
Assume bins are chosen so that signal expecta6on si is uniform in all bins
but
bgd bi is unknown

34

Likelihood: L(s,bi) = e–Ks e–(1+τ)Σbi Πj(s+bj)
K = number of bins
(e.g. 100)
τ = scale factor for bgd(e.g. 1)
j = “on” bins with event(e.g. 1….. 10)
Profile over background nuisance params bi
Lprof(s) has largest value at
s=0
if non < K/(1+τ)
s=non/K
if non ≥ K/(1+τ)
Similar result for Bayesian marginalisation of L(s,bi) over backgrounds bi
i.e. With many bins, profile (or marginalised) L has largest value at s=0,
even though non = 10 and moff=0
BUT when mass distribution ignored (i.e. just counting experiment),
signal+bgd is favoured over just bgd
35

WHY?
Background given greater freedom with large number K of nuisanceparameters
Compare:Neyman and ScoK, “Consistent es6mates based on par6ally
consistent observa6ons”, Econometrica 16: 1-32 (1948)
Data = n pairs X1i = G(µi, σ2)X2i = G(µi, σ2)Param of interest = σ2
Nuisance params = µi. Number increases with n
Proﬁle L es6mate of σ2 are biassed E = σ2/2
and inconsistent (bias does not tend to 0 as n à ∞)MORAL: Beware!
36

p0 v p1 plots
Preprint by Luc Demor6er and LL,“Tes6ng Hypotheses in Par6cle Physics:
Plots of p0 versus p1”
hKp://arxiv.org/abs/1408.6123
For hypotheses H0 and H1, p0 and p1
are the tail probabili6es for data
sta6s6c t
Provide insights on:
CLs for exclusion
Punzi deﬁni6on of sensi6vity
Rela6on of p-values and Likelihoods
Probability of misleading evidence
Sampling to foregone conclusion
Jeﬀreys-Lindley paradox
37CLs = p1/(1-p0) à diagonal line
Provides protec6on against excluding H1 when liKle or no sensi6vity
Punzi deﬁni6on of sensi6vity:
Enough separa6on of pdf’s for no chance of ambiguity
Δµ
H0

H1

t
Can read oﬀ power of test
e.g. If H0 is true, what isprob of rejec6ng H1?
N.B. p0 = tail towards H1
p1 = tail towards H038

Why p ≠ Likelihood ra6o
Measure diﬀerent things:
p0 refers just to H0; L01 compares H0 and H1

Depends on amount of data:
e.g. Poisson coun6ng expt liKle data:
For H0, μ0 = 1.0. For H1, μ1 =10.0
Observe n = 10 p0 ~ 10-7 L01 ~10-5
Now with 100 6mes as much data, μ0 = 100.0 μ1 =1000.0
Observe n = 160 p0 ~ 10-7 L01 ~10+14
39

Jeﬀreys-Lindley Paradox
H0 = simple,
H1 has µ freep0 can favour H1, while B01 can favour H0
B01 = L0 / ∫L1(s) π(s) ds

Likelihood ra6o depends on signal :e.g. Poisson coun6ng expt small signal s:
For H0, μ0 = 1.0. For H1, μ1 =10.0
Observe n = 10 p0 ~ 10-7
L01 ~10-5 and favours H1
Now with 100 6mes as much signal s, μ0 = 100.0 μ1 =1000.0
Observe n = 160 p0 ~ 10-7
L01 ~10+14 and favours H0
B01 involves intergra6on over s in denominator, so a wide enough range
will result in favouring H0
However, for B01 to favour H0 when p0 is equivalent to 5σ, integra6on
range for s has to be O(106) 6mes Gaussian widths
40

WHY LIMITS?
Michelson-Morley experiment à death of aether
HEP experiments: If UL on expected rate for new
par6cle < expected, exclude par6cle
CERN CLW (Jan 2000)
FNAL CLW (March 2000)
Heinrich, PHYSTAT-LHC, “Review of Banﬀ Challenge”
41Methods (no systema6cs)
Bayes (needs priors e.g. const, 1/μ, 1/√μ, μ, …..)
Frequen6st (needs ordering rule,possible empty intervals, F-C)
Likelihood (DON’T integrate your L)
χ2 (σ2 =μ)χ2(σ2 = n)
Recommenda6on 7 from CERN CLW: “Show your L”
1) Not always prac6cal
2) Not suﬃcient for frequen6st methods
42Ilya Narsky, FNAL CLW 2000

43

DESIRABLE PROPERTIES• Coverage• Interval length• Behaviour when n < b• Limit increases as σb increases• Uniﬁed with discovery and interval es6ma6on

44

90% Classical interval for Gaussian
σ=1

Xobs = 3
Xobs = 1
Xobs =-2

μ≥0e.g. m2(νe)

Two-sided range
Upper limit
No region for µ
45

FELDMAN - COUSINS
Wants to avoid empty classical intervals àUses “L-ra6o ordering principle” to resolve
ambiguity about “which 90% region?”[Neyman + Pearson say L-ra6o is best for
hypothesis tes6ng]
Uniﬁed à No ‘Flip-Flop’ problem
46Classical (Neyman) Confidence Intervals
Uses only P(data|theory)

TheoreticalParameterµData xExample:Param = Temp at centre of Sun
Data = Est. flux of solar neutrinosµ≥0

Prob(µl<µ<µu) = α

No prior for µ47

Classical (Neyman) Confidence Intervals
Uses only P(data|theory)

TheoreticalParameterµData x<1.5
1.5 – 2.2
>2.2

µ range
Empty
Upper limit
2-sidedData xExample:Param = Temp at centre of Sun
Data = est. flux of solar neutrinosµ≥0

No prior for µ48

Feldman-Cousins
90% conf intervals

Xobs = -2 now gives upper limit

49

Search for Higgs:
Hà γ γ: low S/B, high statistics

50

HàZ Z à 4 l: high S/B, low statistics

51

p-value for ‘No Higgs’ versus mH

52

Mass of Higgs:
Likelihood versus mass

53

Comparing 0+ versus 0- for Higgs(like Neutrino Mass Hierarchy)

http://cms.web.cern.ch/news/highlights-cms-results-presented-hcp

54

Conclusions
Resources:
Software exists:
e.g. RooStats
Books exist: Barlow, Cowan, James, Lyons, Roe,…..
New: `Data Analysis in HEP: A Practical Guide to
Statistical Methods’ , Behnke et al.
PDG sections on Prob, Statistics, Monte Carlo
CMS and ATLAS have Statistics Committees (and BaBar and CDF
earlier) – see their websites

Before re-inventing the wheel, try to see if Statisticians have alreadyfound a solution to your statistics analysis problem.
Don’t use a square wheel if a circular one already exists.
“Good luck”

55