DRAFT 0.3

Glen Cowan RHUL Physics 15 January, 2012

Asymptotic distribution of the ratio of two profile likelihoods

1 Introduction

In Ref. [1], distributions for test statistics based on likelihood ratios were given, using approximations that become exact in the limit of a large data sample. The main results were given for test statistics based on a ratio of the profile likelihood (defined below) to the maximized likelihood. Another type of likelihood ratio, namely, a ratio of two profile likelihoods corresponding to two hypotheses has also been widely used in HEP searches, especially at the Tevatron. A special case of such a likelihood ratio was also discussed in Sec. 3.8 of [1], but it was restricted to the case where the distributions of the corresponding likelihood-ratio statistic (defined below) was approximately Gaussian. The main purpose of the present note is to extend the result of [1] to cases where the Gaussian approximation for the test statistic is not valid. The more general result is needed when the standard deviation of the estimator of the signal strength depends significantly on the hypothesized signal strength.

Consider a model containing a rate parameter μ as well as nuisance parameters θ = (θ_{1}, …, θ_{m}), described by a likelihood function L(μ, θ). The tests of hypothesized parameter value μ were based on the ratio of the profile to maximized likelihoods:

λ(μ) = \frac{L(μ, ˆ θ̂(μ))}{L(μ̂, θ̂)} .  (1)

Here μ̂, and θ̂ are the values of the parameters that maximize the likelihood and ˆ θ̂(μ) denotes the values of the θ that maximize the likelihood for a given value of μ. That is, the numerator of (1) is the profile likelihood and the denominator is the maximized likelihood.

In this note, we will assume one measures a continuous variable x that follows a Gaussian distribution with mean μ and standard deviation σ_{μ}. This is equivalent to the assumption that the asymptotic approximations of Wald and Wilks hold, as described in Ref. [1]. There, the Maximum Likelihood (ML) estimator for μ, μ̂, played the role of the measured quantity x. We assume that there are two specific hypothesis of interest:

H_{0} : x ∼ Gauss(μ_{0}, σ_{0}) ,  (2)
H_{1} : x ∼ Gauss(μ_{1}, σ_{1}) .  (3)

In a search for a new process, H_{0} typically represents the hypothesis that the data sample consists entirely of events from known (background) processes, and and H_{1} is the hypothesis that includes a sought after signal process.

Often, for example, the measurement is based on counting a number of events n that are Poisson distributed with a mean value μs + b, where s and b represent the expected number of events from signal and background processes, respectively, and μ is a strength parameter defined such that μ = 0 gives the background-only hypothesis H_{0} and μ = 1 corresponds to the nominal alternative H_{1}. In this case the ML estimator for μ is

μ̂ = \frac{n − b}{s} .  (4)

Since n is Poisson distributed with mean μs + b, the variance of n is V [n] = μs + b and therefore the standard deviation of μ̂ is

σ_{μ} = \frac{\sqrt{μs + b}}{s} .  (5)

Here we will assume that specification of the hypotheses H_{0} and H_{1} provides known values for the two mean values μ_{0} and μ_{1} as well as the standard deviations σ_{0} and σ_{1}. Here for illustrative purposes we will take the parameter μ under the two hypotheses H_{0} and H_{1} to be μ_{0} = 0 and μ_{1} = 1, with the corresponding standard deviations

σ_{0} = \frac{\sqrt{b}}{s} ,  (6)

σ_{1} = \frac{\sqrt{s + b}}{s} .  (7)

Although this rule for the standard deviations is derived from the assumption of Poisson distributed data, we nevertheless consider here that the data value x follows a Gaussian distribution, as would be the case for x = μ̂ given a sufficiently large value for the mean of n.

In more complicated problems with nuisance parameters, one would nevertheless have a ML estimator μ̂ for a corresponding signal rate parameter μ, and this will have a certain standard deviation σ_{μ}. To apply the procedure here to such analyses, one must obtain the standard deviation σ_{μ} of μ̂ (i.e., what we call here x) by using, for example, the matrix of second derivatives of the log-likelihood function. The value of σ_{μ} will depend in general on the assumed values of both μ and of any nuisance parameters.

2 Definition of the test statistic q and distribution for σ_{0} = σ_{1}

We assume that the measurement can be represented by a single Gaussian distributed variable x with mean μ and standard deviation σ_{μ}, so that the likelihood function is

L(μ, σ_{μ}) = \frac{1}{\sqrt{2π}σ_{μ}} e^{−(x−μ)^{2}/2σ_{μ}^{2}} .  (8)

The two hypotheses H_{0} and H_{1} defined above thus correspond to two pairs of values for μ and σ_{μ}, namely, (μ_{0}, σ_{0}) and (μ_{1}, σ_{1}).

According to the Neyman-Pearson Lemma (see, e.g., Ref. [2]), when defining a test of a given size α of H_{0}, to obtain the highest power with respect to the alternative H_{1} the optimal test statistic is given by the likelihood ratio L(μ_{1}, σ_{1})/L(μ_{0}, σ_{0}). Equivalently we can use a monotonic function of this and therefore we define the statistic

q = −2 ln \frac{L(μ_{1}, σ_{1})}{L(μ_{0}, σ_{0})} = \frac{(x − μ_{1})^{2}}{σ_{1}^{2}} − \frac{(x − μ_{0})^{2}}{σ_{0}^{2}} + 2 ln \frac{σ_{1}}{σ_{0}} .  (9)

In the special case where the two standard deviations are the same, i.e., σ_{0} = σ_{1} ≡ σ, Eq. (9) becomes

q = \frac{μ_{1}^{2} − μ_{0}^{2} − 2(μ_{1} − μ_{0})x}{σ^{2}} .  (10)

The statistic q is thus a linear function of x and therefore also follows a Gaussian distribution with expectation value and variance given by

E[q|μ, σ_{μ}] = \frac{μ_{1}^{2} − μ_{0}^{2} − 2(μ_{1} − μ_{0})μ}{σ^{2}} ,  (11)

V [q|μ, σ_{μ}] = \frac{4(μ_{1} − μ_{0})^{2}σ_{μ}^{2}}{σ^{4}} .  (12)

These results were given in Eqs. (73) and (74) of Ref. [1] for the special cases of σ_{μ} = σ, μ_{0} = 0 and μ_{1} = μ. The purpose of the present note is to extend this result to the case where σ_{0} and σ_{1} are not equal, in which case the resulting distribution for q is no longer Gaussian.

3 Distribution of q for σ_{0} ≠ σ_{1}

In this section we treat the case of σ_{0} ≠ σ_{1}. Suppose σ_{1} > σ_{0}, as is the case if their values are assigned according to Eqs. (6) and (7). The dependence of q on x given by Eq. (9) is shown in Fig. 1. It is a parabola with a maximum (or minimum, if σ_{0} > σ_{1}) that occurs at

x_{m} = \frac{μ_{0}σ_{1}^{2} − μ_{1}σ_{0}^{2}}{σ_{1}^{2} − σ_{0}^{2}} ,  (13)

and at this point, q takes on the value

q_{m} = \frac{(μ_{1} − μ_{0})^{2}}{σ_{1}^{2} − σ_{0}^{2}} + 2 ln \frac{σ_{1}}{σ_{0}} .  (14)

To find the distribution of q, we first solve for x in terms of q. The quadratic relation (9) has two solutions, x_{+} and x_{−}, given by

x_{±} = x_{m} ± \sqrt{\frac{2(q − q_{m})}{q′′}}  (15)

where

q′′ = \frac{2}{σ_{1}^{2}} − \frac{2}{σ_{0}^{2}}  (16)

Figure 1: The test statistic q as a function of the measured variable x. The curve has been computed using μ_{0} = 0, μ_{1} = 1, σ_{0} = 0.632 and σ_{1} = 0.775, which result from Eqs. (6) and (7) using s = 5 and b = 10.

is the second derivative of q with respect to x. As x is Gaussian distributed with mean μ and standard deviation σ_{μ}, the distribution of q is given by

f(q|μ, σ_{μ}) = |\frac{dx_{+}}{dq} | \frac{1}{σ_{μ}} φ (\frac{x_{+}(q) − μ}{σ_{μ}}) + |\frac{dx_{−}}{dq} | \frac{1}{σ_{μ}} φ (\frac{x_{−}(q) − μ}{σ_{μ}}) ,  (17)

where the function φ is the standard Gaussian pdf (zero mean and unit variance), x_{+} and x_{−} are given by Eq. (15), and the derivatives of x with respect to q for the two solutions x_{+} and x_{−} are

\frac{dx_{±}}{dq} = ±\frac{1}{\sqrt{2(q − q_{m})q′′}} = ±\frac{σ_{0}σ_{1}}{2} \frac{1}{\sqrt{(σ_{0}^{2} − σ_{1}^{2})(q − q_{m})}} .  (18)

Putting together the ingredients and writing the Gaussian terms explicitly gives

f(q|μ, σ_{μ}) = \frac{1}{\sqrt{2(q − q_{m})q′′}} \frac{1}{\sqrt{2π}σ_{μ}} \{exp [−\frac{(x_{m} + \sqrt{(2/q′′)(q − q_{m})} − μ)^{2}}{2σ_{μ}^{2}} ]

+ exp [−\frac{(x_{m} − \sqrt{(2/q′′)(q − q_{m})} − μ)^{2}}{2σ_{μ}^{2}} ]\}  (19)

for q < q_{m} and zero otherwise. Distributions of q according to Eq. (19) are shown in Figs. 2 for μ_{0} = 0, μ_{1} = 1, and with values of σ_{0} and σ_{1} computed using Eqs. (6) and (7) based on the values of s and b shown.

4 Cumulative distribution of q

To compute the p-value of a hypothesized value of μ, one requires the cumulative distribution

F (q|μ, σ_{μ}) = ∫_{−∞}^{q} f(q′|μ, σ_{μ}) dq′  (20)

Because of the parabolic shape of the function q(x), the probability to find q less than a given value is equal to one minus the probability to find x between the corresponding values of x_{−} and x_{+}, i.e.,

(a)  (b)

Figure 2: Distributions of q given values of μ = 0 and μ = 1 using values of σ_{0} and σ_{1} based on the values of s and b shown. These give σ_{0} = 0.632 and σ_{1} = 0.775 for (a) and σ_{0} = 0.5 and σ_{1} = 0.548 for (b).

F (q|μ, σ_{μ}) = 1 − P (x_{−} < x < x_{+}) = 1 − Φ (\frac{x_{+} − μ}{σ_{μ}}) + Φ (\frac{x_{−} − μ}{σ_{μ}}) ,  (21)

where x_{+} and x_{−} are given by Eq. (15) and Φ is the cumulative distribution of the standard Gaussian.

If one has μ_{1} > μ_{0}, as used in the examples here, then the p-value of H_{0} (i.e., μ_{0}, σ_{0}) is the probability, assuming H_{0}, to find q less than or equal to what one observed, i.e.,

p_{0} = F (q|μ_{0}, σ_{0}) .  (22)

In a similar way, the p-value of the H_{1} hypothesis is then the probability, assuming H_{1}, to find q greater than or equal to the value observed,

p_{1} = 1 − F (q|μ_{1}, σ_{1}) .  (23)

References

[1] Glen Cowan, Kyle Cranmer, Eilam Gross and Ofer Vitells, Asymptotic formulae for likelihood-based tests of new physics, Eur. Phys. J. C 71 (2011) 1554.

[2] A. Stuart, J.K. Ord, and S. Arnold, Kendall’s Advanced Theory of Statistics, Vol. 2A: Classical Inference and the Linear Model 6th Ed., Oxford Univ. Press (1999), and earlier editions by Kendall and Stuart.