Feldman and Cousins Confidence Intervals

Feldman and Cousins solved the problem on how to choose confidence intervals in a unified way (that is without basing the choice on the data) [Feldman1998]. The functions gammapy.stats.fc_* give you access to a numerical solution to the Feldman Cousins algorithm. It can be used for any type of statistical distribution and is not limited to a Poisson process with background or a Gaussian with a boundary at the origin.

The basic ingredient to fc_construct_acceptance_intervals_pdfs is a matrix of \(P(X|\\mu)\) (see e.g. equation (3.1) and (3.2) in [Feldman1998]). Every row is a probability density function (PDF) of x and the columns are built up by varying the signal strength \(\\mu\). The other parameter is the desired confidence level (C.L.). The function will return another matrix of acceptance intervals where 1 means the point is part of the acceptance interval and 0 means it is outside. This can be easily converted to upper and lower limits (fc_get_limits), which simply connect the outside 1s for different \(\\mu\) values. An upper or lower limit is obtained by drawing a vertical line at the measured value and calculating the intersection (fc_find_limit).

Examples

Assume you have a Poisson background with known mean 3.0. We generate the matrix of \(P(X|\\mu)\) like this

import gammapy.stats as gstats
import numpy as np
from scipy import stats

x_bins = np.arange(0, 50)
mu_bins = np.linspace(0, 15, 15 / 0.005 + 1, endpoint=True)
matrix = [stats.poisson(mu + 3.0).pmf(x_bins) for mu in mu_bins]

Now we generate the 90% acceptance intervals and construct the lower and upper limit from them:

acceptance_intervals = gstats.fc_construct_acceptance_intervals_pdfs(matrix, 0.9)
LowerLimitNum, UpperLimitNum, _ = gstats.fc_get_limits(mu_bins, x_bins, acceptance_intervals)

Let’s say you measured x = 1, then the 90% upper limit would be:

>>> gstats.fc_find_limit(1, UpperLimitNum, mu_bins)
1.875

The following plot shows the confidence belt based on the Feldman and Cousins principle for a 90% confidence level for the unknown Poisson signal mean \(\\mu\). It is a reproduction of Fig. 7 from [Feldman1998]. It should be noted that the plot in the paper is inconsistent with Table IV from the same paper: the lower limit is off by one bin to the left.

()

../_images/plot_fc_poisson.png

Assume you have an experiment where the observable x is simply the measured value of \(\\mu\) in an experiment with a Gaussian resolution with known width \(\\sigma\). The following plot shows the confidence belt based on the Feldman and Cousins principle for a 90% confidence level for the mean of the Gaussian \(\\mu\), constrained to be non-negative. it reproduces Fig. 10 from [Feldman1998].

()

../_images/plot_fc_gauss.png

Acceptance Interval Fixing

Feldman and Cousins point out that sometimes the set of intersected horizontal lines is not simply connected. These artefacts are corrected with fc_fix_limits.

>>> gstats.fc_fix_limits(LowerLimitNum, UpperLimitNum)

The following script in the examples directory demonstrates the problem: example_fc_demonstrate_artefact.py

For mu = 0.745 the 90% acceptance interval is [0,8] and for mu = 0.750 it is [1,8]. A lot of the fast algorithms that do not compute the full confidence belt will come to the conclusion that the 90% confidence interval is [0, 0.745] and thus the upper limit when zero is measured should be 0.745 (one example is TFeldmanCousins that comes with ROOT, but is has the additional bug of making the confidence interval one mu bin to big, thus reporting 0.75 as upper limit).

For mu = 1.035 the 90% acceptance interval is [0,8] again and only starting mu = 1.060 will 0 no longer be in the 90% acceptance interval. Thus the correct upper limit according to the procedure described in [Feldman1998] should be 1.055, which is also the value given in the paper (rounded to 1.06).

Sensitivity

[Feldman1998] also defines experimental sensitivity as the average upper limit that would be obtained by an ensemble of experiments with the expected background and no true signal. It can be calculated using fc_find_average_upper_limit.

>>> gstats.fc_find_average_upper_limit(x_bins, matrix, UpperLimitNum, mu_bins)
4.41

General Case

In the more general case, one may not know the underlying PDF of \(P(X|\\mu)\). One way would be to generate \(P(X|\\mu)\) from Monte Carlo simulation. With a dictionary of mu values and lists of X values from Monte Carlo one can use fc_construct_acceptance_intervals to construct the confidence belts.

Here is an example, where the X values are generated from Monte Carlo (seed is fixed here, so the result is known):

import gammapy.stats as gstats
import numpy as np
from scipy import stats

x_bins = np.linspace(-10, 10, 100, endpoint=True)
mu_bins = np.linspace(0, 8, 8 / 0.05 + 1, endpoint=True)

np.random.seed(seed=1)

distribution_dict = dict((mu, [stats.norm.rvs(loc=mu, scale=1, size=5000)]) for mu in mu_bins)

acceptance_intervals = gstats.fc_construct_acceptance_intervals(distribution_dict, x_bins, 0.6827)

LowerLimitNum, UpperLimitNum, _ = gstats.fc_get_limits(mu_bins, x_bins, acceptance_intervals)

mu_upper_limit = gstats.fc_find_limit(1.7, UpperLimitNum, mu_bins)
# mu_upper_limit == 2.7

Verification

To verify that the numerical solution is working, the example plots can also be produced using the analytical solution. They look consistent. The scripts for the analytical solution are given in the examples directory: example_fc_poisson_analytical.py example_fc_gauss_analytical.py