A Unified Approach to the Classical Statistical Analysis of Small Signals
MATLAB by Alex Read, March 2021, University of Oslo, Department of Physics
Introduction
To find a unified confidence interval, as defined by Feldman and Cousins, for a gaussian measurement with the mean parameter 
we need to map out the acceptance region based on the likelihood ratio ordering principle.
we need to map out the acceptance region based on the likelihood ratio ordering principle.The likelihood is 
. The log-likelihood is 
. If we consider only then the likelihood is not defined for 
. The log-likelihood is minimized, under the constraint , by 
. The minimized log-likelihood for the best-fit is 
for 
and 
for 
. The log-likelihood ratio is 
for and 
for .
. The log-likelihood is . If we consider only then the likelihood is not defined for . The log-likelihood is minimized, under the constraint , by . The minimized log-likelihood for the best-fit is for and for . The log-likelihood ratio is for and for .Below we see how the log-likelihood ratio looks for various small values of μ. For 
the ratio is parabolic in an acceptance region of 90% or more, therefore the 90% CL F&C confidence intervals will coincide with the intervals defined by the log-likelihood contour defined by 
, where
. For 
the lower end of the acceptance region will be below 0, in the region where the log-likelihood is linear, and the computation of the interval is more complicated. We find numerically the threshold of 
, and the corresponding interval, that yields the specified acceptance (CL). For small enough x it won't be possible to define a 2-sided acceptance region, and there will only be an upper limit. This is the unifying feature of the F&C intervals - the automatic switch from 1-sided to 2-sided confidence intervals, always maintaining explicit coverage for all values of . Coverage means that regardles of the observed x, the derived interval will have a probability of 
of including the true value μ.
the ratio is parabolic in an acceptance region of 90% or more, therefore the 90% CL F&C confidence intervals will coincide with the intervals defined by the log-likelihood contour defined by , where. For the lower end of the acceptance region will be below 0, in the region where the log-likelihood is linear, and the computation of the interval is more complicated. We find numerically the threshold of , and the corresponding interval, that yields the specified acceptance (CL). For small enough x it won't be possible to define a 2-sided acceptance region, and there will only be an upper limit. This is the unifying feature of the F&C intervals - the automatic switch from 1-sided to 2-sided confidence intervals, always maintaining explicit coverage for all values of . Coverage means that regardles of the observed x, the derived interval will have a probability of of including the true value μ.xm = linspace(-3,0,60);
xp = linspace(0,6,60);
figure
mu = 0;
plot(xm,mu^2-2*mu*xm,'k--','LineWidth',2), hold on
plot(xp,(xp-mu).^2,'k-')
mu = 0.5;
plot(xm,mu^2-2*mu*xm,'r--')
plot(xp,(xp-mu).^2,'r-')
mu = 1;
plot(xm,mu^2-2*mu*xm,'g--')
plot(xp,(xp-mu).^2,'g-')
mu = 2;
plot(xm,mu^2-2*mu*xm,'m--')
plot(xp,(xp-mu).^2,'m-')
xlabel('Observed mean $x$','Interpreter','latex')
ylabel('$-2\ln\Lambda(x|\mu)$','Interpreter','latex')
grid on
title('Log-likelihood ratio for $\mu\ge 0$','Interpreter','latex')
ca = gca;
ca.YLim = [0,9];
cs.Xlim = [-3,6];
yline(2.7055);
xline(0);
legend('$\mu=0, x<0$','$\mu=0, x\ge 0$','$\mu=0.5, x<0$','$\mu=0.5, x\ge 0$','$\mu=1, x<0$','$\mu=1, x\ge 0$','$\mu=2, x<0$','$\mu=2, x\ge 0$','Interpreter','latex')
Form the confidence belt
Loop over true value μ for a given 
. Check if the acceptance region α has both sides . If so, we can find the acceptance region by inverting the cdf of the normal distribution since the region is defined by 
. If the lower side has , the F&C prescription is to add regions of equal likelihood (or probability), i.e. we have to solve 
for y, where 
for and 
for , in other words 
and 
.
. Check if the acceptance region α has both sides . If so, we can find the acceptance region by inverting the cdf of the normal distribution since the region is defined by . If the lower side has , the F&C prescription is to add regions of equal likelihood (or probability), i.e. we have to solve for y, where for and for , in other words and .mu = linspace(0,6,100);
mu(1) = 1E-6; % Avoid divide-by-zero problems with mu=0
alpha = 0.90
alpha = 0.9000
fun = @(y,mu,alpha) F(y,mu,alpha);
x1 = zeros(size(mu));
x2 = zeros(size(mu));
for i=1:length(mu)
p = normcdf(2*mu(i),mu(i))- normcdf(0,mu(i));
if p > alpha
% Lower side of acceptance region has x>=0
x1(i) = norminv((1-alpha)/2,mu(i));
x2(i) = norminv(1-(1-alpha)/2,mu(i));
else
% Lower side of acceptance region has x<0
fun = @(y,mu,alpha) F(y,mu,alpha);
y = fzero(@(y) fun(y,mu(i),alpha),5);
x1(i) = max(-5,(mu(i)^2-y)/(2*mu(i)));
x2(i) = mu(i) + sqrt(y);
ilast = i;
end
end
Plot the confidence belt
The part of the belt that coincides with log-likelihood contours is shown in green. The part in black is affected by the constraint .
.figure
plot(x1,mu,'k-'), hold on
plot(x2,mu,'k-')
plot(x1(ilast+1:end),mu(ilast+1:end),'g-')
plot(x2(ilast+1:end),mu(ilast+1:end),'g-')
ca = gca;
ca.XLim = [-2,4];
ca.YLim = [0,6];
xlabel('Measured mean $x$',"Interpreter","latex")
ylabel('Mean \mu')
grid on
title('Feldman and Cousins 90\% CL w/$$\mu\ge 0$$','Interpreter',"latex")
Find limits directly
Find the upper limits as a function of the observed x without referring to tabulated confidence belts. Superimpose some points on the confidence belt to be sure it is done correctly.
xo = -2:0.1:4;
for i=1:length(xo)
muup(i) = find_mu_up(xo(i),alpha);
end
plot(xo,muup,'r+')
Find the lower limits as a function of the observed x without referring to tabulated confidence belts. Superimpose some points on the confidence belt to be sure it is done correctly.
for i=1:length(xo)
mudown(i) = find_mu_down(xo(i),alpha);
end
plot(xo,mudown,'ko')
Superimpose the 
95% and 90% CL upper limits
95% and 90% CL upper limitsx = -2:0.1:4;
ul = iCLsUpperLimit(x,0.95);
plot(x,ul,'b-')
ul = iCLsUpperLimit(x,0.90);
plot(x,ul,'m-')
legend('Belt left','Belt right','Belt left','Belt right','F&C UL','F&C LL','CLs 95% CL UL','CLs 90% CL UL','Location','best')
Reproduce Table X in the F&C article
Extend it to larger 
and 
CL
and CLalphaTab=[0.6827,0.90,0.95,0.99,chi2cdf(5^2,1)];
fprintf('xo 68.27%% CL 90%% CL 95%% CL 99%% CL 5 sigma CL')
for xo = -3:0.1:10
for i=1:length(alphaTab)
mudown(i) = find_mu_down(xo,alphaTab(i));
muup(i) = find_mu_up(xo,alphaTab(i));
end
fprintf('%4.1f %4.2f,%4.2f %4.2f,%4.2f %4.2f,%4.2f %4.2f,%4.2f %4.2f,%4.2f\n',...
xo,mudown(1),muup(1),mudown(2),muup(2),mudown(3),muup(3),mudown(4),muup(4),mudown(5),muup(5))
end
Functions
Solving 
for y allows us to find the confidence belt between and for the case 
. This could also be solved with a inline function.
for y allows us to find the confidence belt between and for the case . This could also be solved with a inline function.function [value] = F(y,mu,alpha)
% Finding the zero of this function allows
% us to find x1(y) and x2(y) for a given alpha.
%
x1 = max(-5,(mu^2-y)/(2*mu));
x2 = mu + sqrt(y);
value = normcdf(x2,mu) - normcdf(x1,mu) - alpha;
end
Calculate the upper limit as directly as possible, i.e. without pre-tabulating the confidence belt results above. The upper limit of the confidence interval is given by the μ that places the observation at the left edge of the confidence belt. If 
then both ends of the confidence interval are within the parabolic part of the log-likelihood ratio and we can solve for 
. For 
we use that 
and and solve 
for μ. Make sure the interval that we search for the solution is more than wide enough to cover the solution! The upper limit will be in the neighborhood of 
plus the number of standard deviations corresponding to the confidence level (alpha).
at the left edge of the confidence belt. If then both ends of the confidence interval are within the parabolic part of the log-likelihood ratio and we can solve for . For we use that and and solve for μ. Make sure the interval that we search for the solution is more than wide enough to cover the solution! The upper limit will be in the neighborhood of plus the number of standard deviations corresponding to the confidence level (alpha).function [mu] = find_mu_up(x1,alpha)
interval = [0,max(0,x1)+1.1*sqrt(chi2inv(alpha,1))];
if x1>=0
fun = @(mu,x1,alpha) normcdf(x1,mu)-(1-alpha)/2;
mu = fzero(@(mu) fun(mu,x1,alpha),interval);
else
fun2 = @(mu,x1,alpha) normcdf(mu+sqrt(mu^2-2*mu*x1),mu)-normcdf(x1,mu)-alpha;
mu = fzero(@(mu) fun2(mu,x1,alpha),interval);
end
end
Calculate the lower limit as directly as possible, i.e. without pre-tabulating the confidence belt results above. First of all, if 
, then we know the lower limit must be 0. Well above 
, if 
, then we know that both ends of the confidence interval are within the parabolic part of the likelihood ratio and we can solve 
for μ. In the remaining intermediate region we must solve for 
for μ.
, then we know the lower limit must be 0. Well above , if , then we know that both ends of the confidence interval are within the parabolic part of the likelihood ratio and we can solve for μ. In the remaining intermediate region we must solve for for μ.function [mu] = find_mu_down(x2,alpha)
p = normcdf(x2,0);
if p<alpha
mu = 0;
else
interval = [0,max(1,x2)];
p = normcdf(x2,x2/2)-normcdf(0,x2/2);
if p>=alpha
fun = @(mu,x2,alpha) normcdf(x2,mu,'upper')-(1-alpha)/2;
mu = fzero(@(mu) fun(mu,x2,alpha),interval);
else
fun = @(mu,x2,alpha) normcdf(x2,mu)-normcdf((2*x2*mu-x2^2)/(2*mu),mu)-alpha;
mu = fzero(@(mu) fun(mu,x2,alpha),interval);
end
end
end
Calculate the upper limit for the unit Gaussian pdf (use this internal copy to make the live script self-contained).
upper limit for the unit Gaussian pdf (use this internal copy to make the live script self-contained).function [uplim,pmu,ompb] = iCLsUpperLimit(x,CL)
%iCLsUpperLimit Compute CLs upper limit.
% [uplim,pmu,ompb] = CLsUpperLimit(x,CL) finds the upper limit
% at the CL confidence level using CLs=pmu/(1-pb) for the unit
% gaussian measurement x with the null (background-only) hypothesis
% at 0 and requiring the signal to be positive. The p-value for the
% signal plus background hypothesis, pmu, and ompb, the p-value for the
% null hypothesis, often called 1-p_b, are also returned.
%
% If x is a list of observations, lists of upper limits and
% p-values are returned.
ompb = normcdf(x,0);
for i=1:length(x)
uplim(i) = fzero(@(mu) normcdf(x(i),mu)/ompb(i)-(1-CL),x(i));
pmu(i) = normcdf(x(i),uplim(i));
end
end