View source for Physics:PP/11/Statistical Significance and Discovery

{{sidebar box|[[Physics:PP/11|<< Back]]}}

{{Notice | Use [[Software:Jas4pp|Jas4pp]] or [[Software:DataMelt|DataMelt]] to process examples of this section.}}

== Statistical Significance and Discovery ==

Experimental measurements often  have the disadvantage that they might have been caused by pure statistical accidents.
In such cases we need  to estimate the level of statistical significance (credibility) of observations,
to ensure that a likelihood of such observations not happening by chance is
small. One can approach this question from the different direction by
estimating  
the probability that the observed results have  occurred due to 
statistical accident. Then, a large value of this probability 
would represent a small level of statistical significance (and vice versa).
For a good recent review see <cite>signif_book</cite>.

The probability we are talking about  is defined by the number called p-value.
Mathematically speaking, it  tells about the probability of a 
result being observed, given that a certain statement (the null hypothesis) is true. 
The p-value is close to 1 when results observed in measurements
could have occurred by chance. 
Observations are statistically significant when the p-value is small, i.e. 
when one can demonstrate that the probability of obtaining 
observations by chance only is relatively low. 

Let us give an example. Historical records on sunspot (Sunspots are relatively dark areas on the surface (photosphere) of the Sun.) 
activity suggest that solar activity is at the level of 80 sunspots per year.
The population of sunspots is normally distributed with the standard deviation of <math>\sigma=20</math>,  and 
the records are large enough to assume the Central Limit Theorem.
During one year, however, 130 sunspots were observed. Is this finding significant? Does it indicate
a beginning of unusual sun activity? 

To estimate the likelihood of the observation of 130 (or more) sunspots in a given year, 
we will calculate the area under the Normal (Gaussian) probability density function 
from minus infinity to 130. The Normal distribution should have the mean 80 and <math>\sigma=20</math>.
Then this  value should be subtracted from 1 in order to calculate the probability for events with <math>>130</math> sunspots. 

<jcode>
from org.apache.commons.math3.distribution import *
no=NormalDistribution(80, 20)
print "p-value=",1-no.cumulativeProbability(130)
</jcode>

The output probability is 0.0062. It tells that the significance of the observation is very high: 
Only 6 years out of 1000  have occurrences with 130 sunspots due to statistical nature.
This value was obtained by using an analytic description of the Normal distribution as 
implemented in the Apache Common Math library. In reality,
the underlying distribution (and asked questions on significance) can be rather complex, therefore, one needs numeric methods 
to estimate statistical significance. 

One can run an experiment where many events are generated 
using a Normal distribution with the mean 80 and <math>\sigma=20</math>, and counting a fraction of events with equal or more 130 values.
Here is the code that numerically estimates the probability and plots a histogram with the null hypothesis.

Here is the numeric calculation of the p-value:

<jcode>
from cern.jet.random.engine import *
from cern.jet.random import *
from jhplot import *

nr=Normal(80,20,MersenneTwister())
h1=H1D("Normal",100,0,200)
sum=0.; total=100000
for i in range(total):
    va=nr.nextDouble()
    h1.fill(va)
    if (va>=130): sum +=1.0
c1=HPlot()
c1.visible(); c1.setAutoRange()
c1.draw(h1)
print "p value=",sum/total
</jcode>

The code prints the expected p-value of <math>\simeq 0.0062</math>. 
The actual precision will depend on the number of random numbers generated.
One should increase the number of random events to reduce statistical uncertainty
on this estimate. 

Now you know the main trick in the numeric estimation of the p-values.
If you need to calculate the p-value numerically, first one should
make multiple observations and determine the underlying distribution function, or your "null hypothesis".
Then run numeric tests and  estimate how many occurrences are expected for a total (large) number of observations,
similar to those in the above example.  The only difficult part is to find the correct null hypothesis.

You may ask the oppose question: How many events we should see when we require a p-value less than  0.05, which is 
rather common question in many fields of science (biology, psychology, physics, etc.). 
We should calculate the value, x, for which the area under the Gaussian probability density function 
(integrated from minus infinity to x) is equal to the argument:

<jcode>
from org.apache.commons.math3.distribution import *
no=NormalDistribution(80, 20)
print "p=0.05:",no.inverseCumulativeProbability(0.95)
</jcode>

This example returns 112.89. 

Now let us consider one special case.
If events are rare, i.e. when  zero is a very likely number to observe, 
we should use a Poisson distribution instead of the Normal distribution.
In the example with the sunspots, one can ask a question of how many sunspots occur with a certain size or type.
For example, if we are talking about rare events with well known mean 0.5, then an 
observation of 3 events can already be quite significant.
Let us calculate the p-value in this situation:

<jcode>
from org.apache.commons.math3.distribution import *
no=PoissonDistribution(0.5)  # mean to 0.5 
print "p-value=",1-no.cumulativeProbability(3) # 3 events
</jcode>

The code returns rather small  probability of 0.0018.

In real life,  expectations for experimental outcomes have some uncertainties.  
Let us calculate a statistical significance for an experiment with 7 events,
while the expectation is  <math>3.1\pm 0.3 </math>, where 0.3 is uncertainty.  In many situations, the uncertainty can be treated 
in the frequentist approach, i.e. we will assume that 0.3 is the standard deviation of a Gaussian distribution that has
the mean value 3.1 (which is often true if the prediction is obtained using different probabilistic techniques). 
In order to account for uncertainty in the prediction, the mean value of the Poisson distribution should be
obtained from a Gaussian distribution with the mean and width of the prediction <cite>Highland</cite>.

<jcode>
from org.apache.commons.math3.random import *

r=RandomDataGenerator()
observed=7.        # observed number of events
bkg=3.1            # expected value 
err=0.3            # expected error
tot=1000000; sum1=0;
for i in range(tot):
   mean=r.nextGaussian(bkg, err)
   value=r.nextPoisson(mean)
   if (value>=observed): sum1 +=1.0
print "p=", sum1/tot
</jcode>

This code returns 0.041. Note that the above example is not very abstract - it was taken from the real-life situation of the top-quark discovery
when the prediction was called "expected background rate", while the number 7 was the total number of observed top-quark candidates <cite>pekka</cite>. 
 
In the case if you are interested in more conservative  estimate of the p-value, and have doubts  in the
statistical nature of the uncertainty, remove the Gaussian distribution and shift the mean value of the Poisson distribution by the
uncertainty,  i.e. evaluate the Poisson distribution with the mean 3.4.

== Discovery sensitivity ==

It should be said that, in many ares of science,  events in a counting experiments 
are considered to be statistically significant
when the p-value is less than 0.05. 
This means that an observation has more than 95% chance of being true, i.e. not caused by pure chance. 
In some scientific fields, such as particle physics, a discovery is announced only when the p-value
is less than  <math>3\cdot 10^{-7} </math> (assuming one-tailed p-value).
This value corresponds to the five-sigma ( <math>5\sigma </math>) observation, i.e.
this probability can be obtained by integrating the tail of the Normal distribution
above 5 standard deviations.
Putting in other words,   
the observation would occur by chance alone as rarely as a value sampled from a 
Gaussian distribution would be five standard deviations from the mean. 
The Gaussian significance, which will be denoted as Z, is related to
the p-value as  <math>Z=\Phi^{-1}(1-p)</math>, where  <math>\Phi^{-1}</math> is the inverse of the cumulative distribution.

Let us see how we can covert the p-value to the Z value in the expression  <math>Z\cdot\sigma </math>:

<jcode>
from org.apache.commons.math3.distribution import *
p=NormalDistribution()
print "Z=",-p.inverseCumulativeProbability(3E-07)
</jcode>

You will see a value Z close to 5, i.e. we can say  <math>5\sigma</math> observation.

If the p-value is smaller than  0.003 (which corresponds to
three sigma,  <math>3\sigma </math>), one can report an evidence for a new phenomena (or a particle)  being observed.
The reason for such a strong requirement by particle physicists  is that new laws of physics should be discovered with much stronger evidence than
what is usually required for the p-value of 0.05 which still lives a good room for statistical fluctuations and various instrumental anomalies.
In addition, the five-sigma rule is a good protection against problems in estimating the p-value in real life
when the background model (or the null hypothesis) is not known well, and measurements can be effected by 
uncertainties that are hard to foreseen. Such uncertainties do not have statistical nature, and they 
are usually called systematic uncertainties.

Now let us come back to the previous example and calculate the p-value that corresponds 
to three sigma ( <math>3\sigma</math>) and five sigma ( <math>5\sigma </math>).
Then we estimate the expected numbers of events we should observe to claim 3 and 5  <math>\sigma </math> discoveries: 

<jcode>
from org.apache.commons.math3.distribution import *
mean=80; sd=20
no=NormalDistribution(mean, sd)
p1=no.cumulativeProbability(mean+3*sd)
p2=no.cumulativeProbability(mean+5*sd)
print "p-value for 3 sigma:",1-p1
print "Events:",int(no.inverseCumulativeProbability(p1))
print "p-value for 5 sigma:",1-p2
print "Events:",int(no.inverseCumulativeProbability(p2))
</jcode>


This prints the answer:
 
<jcode>
p-value for 3 sigma: 0.00135
Events: 139
p-value for 5 sigma: 2.86651E-07
Events: 180
</jcode>
 
As you can see,  the claim of a new phenomena that causes high sun activity
requires the observation of 180 sunspots in a single year assuming the <math>5\sigma</math>  approach. 
Observing 140 sunspots in a single year is good enough for an evidence of an unusual sun activity 
assuming the  <math>3\sigma </math> criteria. 
As explained before, simple situations in  which the null-hypothesis
is well established and systematical variations (i.e those that do not occur by chance) are well understood, 
do not require the  <math>5\sigma</math> extreme, i.e. values of p=0.05 and p=0.01 should be usually sufficient.


This section is discussed in more detail in the book "Numeric Computation and Statistical Data Analysis on the Java Platform" <cite>chekanovbook</cite> by S.Chekanov.

{{Author|Sergei Chekanov|Schekanov}}

==References==
{{Reflist}}