DMelt:Statistics/4 Statistical Limits

Table of contents

Statistical limits

Here we will shows a few examples of how to set statistical limits.

Let us build 2 experiments with observed data, and calculate probability for the given limit. The user can iterate towards whatever limit type they choose, e.g., for a 90% upper limit they adjust the limit guess until the probability is shown as 0.10; for a 90% lower limit they would aim for 0.90.

# Build 2 experiments with observed data, and calculate
# probability for the given limit.  The user can  iterate towards 
# whatever limit type they choose, e.g., for a 90% upper 
# limit they adjust the limit guess until the probability is shown as 0.10; 
# for a 90% lower limit they would aim for 0.90.

from jhpro.stat.limit import *

sensitivity=1;
sigma=0;
limitguess=90;
MCevents=10000;
s=StatConfidence(sensitivity,sigma,limitguess,MCevents);

# build an experimental data
# - 200 observed events
# - 140 expected background
# - 20  error on expected background 
exp1=ExpData(200,140,20)
s.addData(exp1)

# second experiment with 150 observed events
exp2=ExpData(190,190,20)
s.addData(exp2)

print "-> Get probabilities for upper limits:"
s.runUpper()
print "Prob (Cousins+Highland)=",s.getProbabilityCH()   
print "Prob (BaBar SWG)=",s.getProbabilitySWG()
print "Prob (Jeffreys)=",s.getProbabilityJ()

# run for lower limit
print "\n->  Get probabilities for lower limits:"
s.runLower()
print "Prob (Cousins+Highland)=",s.getProbabilityCH()              
print "Prob (BaBar SWG)=",s.getProbabilitySWG()
print "Prob (Jeffreys)=",s.getProbabilityJ()

The output is:

Get probabilities for lower limits:
Prob (Cousins+Highland)= 0.6716
Prob (BaBar SWG)= 0.6719
Prob (Jeffreys)= 0.6671

Signal exclusion limits on a smooth background

This program demonstrates the computation of 95 % C.L. limits with correct treatment of statistical errors. Signal hypothesis is excluded at the 95% CL if CLs = 0.05 and at more than the 95% CL if CLs < 0.05, assuming that signal is present.

This algorithm uses the Likelihood ratio semi-Bayesian method. This method is described in this paper (PDF). It takes signal, background and data histograms and runs a set of Monte Carlo experiments (>50000 in this example ) in order to compute the limits.

Let's give some definitions: [math]\displaystyle{ CL_{s+b} =P_{s+b} (X \lt X_{obs}) }[/math] is the confidence level for excluding the possibility of simultaneous presence of new particle production and background (the s + b hypothesis). This corresponds to the probability that the test statistic would be less than or equal to that observed in the data, assuming the presence of both signal and background at their hypothesized levels.

The confidence level [math]\displaystyle{ (1-CL_{s+b}) }[/math] may be used to quote exclusion limits although it has the disturbing property that if too few candidates are ob- served to account for the estimated background, then any signal, and even the background itself, may be excluded at a high confidence level. [math]\displaystyle{ CL_{b} =P_{b} (X \lt X_{obs}) }[/math] is the confidence level for the background alone. It expresses the probability that background processes would give fewer than or equal to the number of candidates observed. Then [math]\displaystyle{ 1-CL_{b} }[/math] gives the probability that the background could have fluctuated to produce a distribution of candidates at least as signal-like as those observed in the data. [math]\displaystyle{ CL_{s} = CL_{s+b}/CL_{b} }[/math] is the Modified Frequentist confidence level. It is used to calculate 95% CL upper limits (CLs less than 0.05) on the signal

Furthermore, the PDFs of X in the signal+background and background hypotheses allow computation of the expected confidence levels <CLs+b>, <CLb>, and <CLs>, assuming the presence only of background. These are indications of how well an experiment would do on average in excluding a signal if the signal truly is not present, and are the important figures of merit when optimizing an analysis for exclusion.

"""
Theory: Signal hypothesis is excluded at the 95% CL if CLs = 0.05
and at more than the 95% CL if CLs < 0.05, assuming that signal is present

Authors:  Christophe.Delaere@cern.ch on 21.08.02 and Sergei Chekanov
"""

from java.awt import Color,Font
from java.util import Random
from jhplot import *
from jhpro.stat import *

c1 = HPlot("Canvas",600,400)
c1.setGTitle("Computation of 95 % C.L. limits")
c1.visible()
c1.setRange(-4,4,0.0,120)
c1.setNameX("Variable")
c1.setNameY("Events")


# set
background = H1D("Background",30,-4.0,4.0)
background.setColor(Color.green)
background.setFill(1)
background.setFillColor(Color.green)
background.setErrAll(0)

signal     = H1D("Signal",30,-4.0,4.0)
signal.setFill(1)
signal.setFillColor(Color.red)
signal.setColor(Color.red)

data     = H1D("Data",30,-4.0,4.0)
data.setColor(Color.black)
data.setStyle("p")

r=Random()
for i in range(25000):
      background.fill(r.nextGaussian(),0.02)
      signal.fill(1+0.2*r.nextGaussian(),0.001)
for i in range(500):
      data.fill(r.nextGaussian(),1.0)    

sigback=background.oper(signal,"Signal+Background","+")
sigback.setErrAll(0)

c1.draw(signal)
c1.draw(background)
c1.draw(sigback)
c1.draw(data)

datasource = DataSource(signal, background,data)
climit = CLimits(datasource, 100000)
confidence = climit.getLimit()
print "CLs    : " ,confidence.getCLs()
print "CLb    : " ,confidence.getCLb()
print "CLsb    : " ,confidence.getCLsb()
print "expected <CLs>    : " ,confidence.getExpectedCLs_b()
print "expected <CLb>    : " ,confidence.getExpectedCLb_b()
print "expected <CLsb>    : " ,confidence.getExpectedCLb_b()

print "Signal hypothesis is excluded at level (%) ", (1-confidence.getCLs())*100.

# export to some image (png,eps,pdf,jpeg...)
# c1.export(Editor.DocMasterName()+".png");
# edit the image
# IEditor(Editor.DocMasterName()+".png");

The output is:

CLs    :  0.00129900435575
CLb    :  0.07939
CLsb    :  0.000103127955803
expected <CLs>    :  0.018271281974
expected <CLb>    :  0.50001
expected <CLsb>    :  0.50001
Signal hypothesis is excluded at level (%)  99.8700995644

Calculating exclusion limits for a new theory in science

In many fields of science, it is important to understand the relevance of new theories or hypotheses in a description of experimental data, assuming that such data are already well represented by predictions of some well-accepted theory. A popular statistical method for setting upper limits (also called exclusion limits) on model parameters of a new theory is based on the CLs_method (particle physics).

Let us show how to set an exclusion limit in case of a counting experiment. We will use particle physics, as an example of hardcore science, where a new theory is considered to be excluded at the 95% confidence level (CL) if CLs = 0.05, and at more than the 95% CL if CLs< 0.05, assuming the presence of a signal that represents a new theory.

In this small tutorial, we will set the exclusion limits at the 95% CL on empirical distributions ("histograms"). We will use 3 histograms: (1) first histogram represents counts of actual measurements, (2) second histogram is the prediction of a well-established theory, (3) third distribution is a prediction ("signal") from some new theory.

Let us generate these histograms using one-tailed Gaussian distributions. The first histogram includes 10,000 events distributed according to one-tailed Gaussian distribution with the standard deviation 1. The second histogram represents our expectation for a known theory ("background"). The latter histogram has a similar shape and the number of events. The third histogram ("signal") corresponds to the prediction of a new theory. We assume that the new theory contributes to the tail of the experimental data (near x=2). We want to find the maximum number of events of signal (which is a parameter of this new theory) which can be reliable excluded, assuming that the given experimental data is already well described by the well-established theory.

The Python/Jython code below creates all the required histograms, scans through the signal strength, and calculates the CLs exclusion limits for a given signal strength:

from java.awt import Color
from java.util import Random
from jhplot import *
from jhpro.stat import *

data=H1D("Fake data",30,0.0,4.0)
data.setColor(Color.black); data.setStyle("p")

backg=H1D("Backg",30,0.0,4.0)
backg.setColor(Color.green); backg.setFill(1); backg.setErrAll(0)
backg.setFillColor(Color.green)

signal= H1D("Signal",30,0.0,4.0)
signal.setFill(1); signal.setFillColor(Color.red)
signal.setColor(Color.red)

r=Random(2L) # make random numbers reproducible
Max=10000    # events in empirical data
strength=10  # initial signal strength at start of the scan

for i in xrange(Max): # fill data and prediction
    data.fill(abs(r.nextGaussian()))
    backg.fill(abs(r.nextGaussian()))

for ev in range(100):  # scan through the parameter strength
   strength=strength+5 # increment Nr of signal events by 10
   signal.clear()
   for i in xrange(strength):
        signal.fill(abs(2+0.3*r.nextGaussian()))
   events=signal.sumAllBinHeights()
   datasource = DataSource(signal, backg, data)
   print "Calculating ..  Nr of signal events:",events
   climit = CLimits(datasource,50000)
   conf = climit.getLimit()
   valCLs=conf.getCLs(); print "  CLs = ",valCLs
   if  valCLs == float("inf"): continue
   exc=int((1-valCLs)*100)
   print "  Signal is excluded at level ",exc,"%"  
   del climit
   if (valCLs<0.05):
       print "Stop! Nr ",events," needed  for ",exc,"% CL"
       break

c1 = HPlot()  # show the results
c1.visible(); c1.setRange(0,3,0.0,1500)
sigback=backg.oper(signal,"Signal+Backg","+")
sigback.setColor(Color.blue)
sigback.setErrAll(0)
c1.draw([backg,signal,sigback,data])
c1.export("limit.pdf")

The above code determines how much of the signal should by allowed (in terms of the event counts) so it can be excluded at 95% confidence level, making this new theory irrelevant to the empirical data. We start from 10 initial events for the signal, calculate the exclusion limit using the class CLimit that uses a Monte Carlo method for settings statistical limits, and then increment the count number of the signal by 5. We stop the calculation when the signal is excluded at 95% CL. This corresponds to about 100 signal events. One can decrease the step size and increase the number of Monte Carlo experiments in order to get a better precision.

The created figure ("limit.pdf") shows the hypothetical "experimental" data overplayed with the blue histogram that includes both the well-established theory and the new theory (excluded at 95% CL). The signal distribution excluded at 95% CL is shown in red.