%% A Unified Approach to the Classical Statistical Analysis of Small Signals
% Gary J. Feldman, Robert D. Cousins, Phys.Rev.D57:3873-3889,1998, arXiv:physics/9711021
% 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 \mu\ge 0 we need to map out the acceptance region based on the likelihood ratio ordering principle.
% The likelihood is L(x|\mu) = {1\over\sqrt{2\pi}} e^{-(x-\mu)^2/2}. The log-likelihood is -2\ln L=(x-\mu)^2+\ln 2\pi. If we consider only \mu\ge 0 then the likelihood is not defined for \mu<0. The log-likelihood is minimized, under the constraint \mu\ge 0, by \hat\mu = \max(x,0). The minimized log-likelihood for the best-fit is \ln 2\pi for x\ge 0 and \ln 2\pi +x^2 for x<0. The log-likelihood ratio is -2\ln\Lambda(\mu) = (x-\mu)^2 for x\ge 0 and (x-\mu)^2-x^2=\mu^2-2\mu x for x<0.
% Below we see how the log-likelihood ratio looks for various small values of \mu. For x>\sim 2 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 -2\ln\Lambda(x|\mu)<2.7055, whereP(\chi^2_1>2.7055) = 0.90. For x<\sim 2 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 -2\ln\Lambda(x|\mu), 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 \mu\ge 0. Coverage means that regardles of the observed x, the derived interval will have a probability of CL of including the true value \mu.

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 \mu for a given CL=\alpha. Check if the acceptance region \alpha has both sides x\ge 0. If so, we can find the acceptance region by inverting the cdf of the normal distribution since the region is defined by \mu\pm k(\alpha). If the lower side has x<0, the F&C prescription is to add regions of equal likelihood (or probability), i.e. we have to solve \int_{x_1(y)}^{x_2(y)}P(x|\mu)dx = \alpha for y, where y=(x_2-\mu)^2 for x\ge 0 and y=(x-\mu)^2-x^2 = \mu^2-2\mu x_1 for x<0, in other words x_1=(\mu^2-y)/(2\mu) and x_2=\mu+\sqrt{y}.


mu = linspace(0,6,100);
mu(1) = 1E-6; % Avoid divide-by-zero problems with mu=0
alpha = 0.90

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 \mu\ge 0.

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 CL_s 95% and 90% CL upper limits

x = -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 x_o and 5\sigma CL

alphaTab=[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 \int^{\mu+\sqrt{y}}_{(\mu^2-y)/(2\mu)}P(x|\mu)dx=\alpha for y allows us to find the confidence belt between x_1=(\mu^2-y)/(2\mu) and x_2=\mu+\sqrt{y} for the case x_1<0. 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 \mu that places the observation x_o at the left edge of the confidence belt. If x_o\ge 0 then both ends of the confidence interval are within the parabolic part of the log-likelihood ratio and we can solve for \int_{-\infty}^{x_o}P(x|\mu)dx=(1-\alpha)/2. For x_o<0 we use that x_o = (\mu^2-y)/(2\mu) and x_2 = \mu+\sqrt{y} and solve \int_{x_o}^{\mu+\sqrt{\mu^2-2\mu x_o}}P(x|\mu)dx=\alpha for \mu. 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 \hat\mu\ge 0 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 \int_{-\infty}^{x_o}P(x|\mu=0)dx<\alpha, then we know the lower limit must be 0. Well above x_o=0, if \int_{0}^{x_o}P(x|\mu=x_0/2)dx>\alpha, then we know that both ends of the confidence interval are within the parabolic part of the likelihood ratio and we can solve \int_{x_o}^\infty P(x|mu)dx=(1-\alpha)/2 for \mu. In the remaining intermediate region we must solve for \int^{x_o}_{(2x_o\mu-x_o^2)/(2\mu)}P(x|\mu)dx = \alpha for \mu.

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 CL_s 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