Corrigendum: The loudest event statistic: general formulation, properties and applications

During the internal collaboration review of the search for gravitational waves from low mass compact binary coalescence in LIGO's sixth science run and Virgo's science runs 2 and 3 [2], an inconsistency in results obtained from two different methods of marginalization over uncertainties in the measured value of the search's detection efficiency revealed there was an error in section 5.1, equations (23) and (24), of [1]. Here we provide a correct version of section 5.1 of [1]. At the end we give the trivial correction to section 5.2.

Our discussion begins with equation (11) of [1], which expresses the posterior distribution for the mean number of true (astrophysical) events μ with a ranking statistic above some largest observed loudness parameter $\hat{x}$ (the 'loudest event') as

$$\begin{equation} p(\mu \mid \hat{\epsilon }, \hat{\Lambda }) = \frac{ p(\mu )(1+\mu \hat{\epsilon }\hat{\Lambda })\,\mathrm{e}^{-\mu \hat{\epsilon }} }{ \int {\rm d}\mu \, p(\mu )(1+\mu \hat{\epsilon }\hat{\Lambda })\,\mathrm{e}^{-\mu \hat{\epsilon }} } \end{equation} \tag{ 1 }$$

where $\hat{\epsilon }=\epsilon (\hat{x})$ is the average detection efficiency function evaluated at the ranking statistic value of the loudest event observed, and $\hat{\Lambda }=\Lambda (\hat{x})$ is a measure of the relative probability of the loudest event arising from a foreground signal versus such an event occurring due to the experimental background; a value $\hat{\Lambda }=0$ results if the loudest event is definitely from the background, and $\hat{\Lambda }\rightarrow \infty$ in the limit that the loudest event is definitely from the foreground. The function Λ(x) is defined by equation (12) of [1]. A Poisson-distributed foreground rate with mean number of events μ is assumed, while p(μ) is the prior unconditional probability distribution for this quantity. The detection efficiency (x) is a monotonic function of the ranking statistic x that describes the probability that a single foreground event provided by nature will have a ranking statistic at least as significant as the value x; therefore, it is equivalent to describe events in terms of the detection efficiency variable , and, in particular, the loudest event will have the value $\hat{\epsilon }$, which is the smallest value of that occurs during the experiment.

Equation (1) is obtained from the likelihood function,

$$\begin{equation} p(x \mid \mu , B) = p_0(x) [1+\mu \epsilon (x)\Lambda (x)] \,\mathrm{e}^{-\mu \epsilon (x)}, \end{equation} \tag{ 2 }$$

via the application of Bayes's theorem,

$$\begin{equation} p(\mu \mid \hat{\epsilon }, \hat{\Lambda }) = \frac{ p(\mu ) p(\hat{x} \mid \mu ,B) }{ \int {\rm d}\mu \, p(\mu ) p(\hat{x} \mid \mu ,B) }, \end{equation} \tag{ 3 }$$

where B indicates dependence of the probability distribution on the background. The function p₀(x) is the unconditional probability distribution of the experimental background producing an event having a loudest event x, while the function p(x∣μ, B) is the probability distribution for x in an experiment having both a background and a foreground as described above. All dependence of the posterior distribution on the background is contained in the value $\hat{\Lambda }$.

Section 5.1 of [1] was concerned with the situation in which the value of $\hat{\epsilon }$ is not known exactly⁶. Suppose that the true value of efficiency is related to the measured value via a factor α so that

$$\begin{equation} \epsilon (x) = \alpha \epsilon _{{\rm meas}}(x). \end{equation} \tag{ 4 }$$

Given an observed loudest event, and consequently a value of $\hat{\epsilon }_{{\rm meas}}=\epsilon _{{\rm meas}}(\hat{x})$, we wish to infer the value of the efficiency and also the posterior distribution for μ. We suppose that, for a given $\hat{\epsilon }_{{\rm meas}}$, the probability distribution for the true (but unknown) efficiency $\hat{\epsilon }$ is $p(\hat{\epsilon } \mid \hat{\epsilon }_{{\rm meas}})$. For simplicity, and in keeping with [1], we choose the gamma-distribution,

$$\begin{equation} p(\hat{\epsilon }; k, \theta ) = \frac{1}{\theta ^k \Gamma (k)}\, \hat{\epsilon }^{k-1} \,\mathrm{e}^{-\hat{\epsilon }/\theta } , \end{equation} \tag{ 5 }$$

where the two parameters θ and k determine the mean, kθ, and the variance, kθ² (cf equation (22) of [1]). If we assume that the mean of this distribution is equal to the measured value of the efficiency $\hat{\epsilon }_{{\rm meas}}$, then we fix the scale parameter to be $\theta = \hat{\epsilon }_{{\rm meas}}/k$ to obtain $p(\hat{\epsilon }\mid \hat{\epsilon }_{{\rm meas}},k)=p(\hat{\epsilon }; k, \theta =\hat{\epsilon }_{{\rm meas}}/k)$; the shape parameter k then specifies the width of the distribution. When k is large, the gamma-distribution will be sharply-peaked about its mean value. Equation (5) can be re-expressed as a probability distribution for the scale factor α in equation (4) as

$$\begin{equation} p(\alpha ; k) = \frac{1}{\Gamma (k)}\,k^k \alpha ^{k-1} \,\mathrm{e}^{-k\alpha }. \end{equation} \tag{ 6 }$$

The scale factor α is a random variable describing the uncertainty in the measurement of the efficiency. We must now derive a modified version of equation (1) that accounts for this additional random variable. We may either (i) construct a marginal likelihood by marginalizing over this uncertainty in the likelihood function and then proceed to use Bayes's theorem to obtain a posterior distribution for μ, or (ii) treat α as a nuisance parameter in the posterior distribution (see, e.g., section 36.1.4.2 of [3]). These should yield the same result, but it is instructive to derive the result using both methods.

We employ method (i) first. We use equations (2) and (4) to write the likelihood function

$$\begin{equation} p(x \mid \mu , \alpha , B) = p_0(x) [1+\mu \alpha \epsilon _{{\rm meas}}(x)\Lambda (x)] \,\mathrm{e}^{-\mu \alpha \epsilon _{{\rm meas}}(x)}. \end{equation} \tag{ 7 }$$

We now construct a marginal likelihood, p(x∣μ, k, B), by marginalizing over α,

$$\begin{equation} p(x\mid \mu ,k,B) = \int _0^\infty {\rm d}\alpha \, p(\alpha ;k) p(x\mid \mu ,\alpha ,B), \end{equation} \tag{ 8 }$$

and using our gamma-distribution for α we find

$$\begin{eqnarray} &&\fl p(x\mid \mu ,k,B) = p_0(x) \frac{k^k}{\Gamma (k)} \nonumber \\ &&\quad\!\!\times \int _0^\infty {\rm d}\alpha \, \alpha ^{k-1} [1+\mu \alpha \epsilon _{{\rm meas}}(x)\Lambda (x)] \,\mathrm{e}^{-\alpha [k+\mu \epsilon _{{\rm meas}}(x)]} \nonumber \\ &&= p_0(x) \frac{ 1 + \mu \epsilon _{{\rm meas}}(x) [\Lambda (x) + 1/k] }{ (1 + \mu \epsilon _{{\rm meas}}(x) / k)^{k+1} }. \end{eqnarray} \tag{ 9 }$$

Note that in the limit k → ∞ we recover equation (2) with (x) replaced with _meas(x). We now obtain a posterior distribution for μ from Bayes's theorem,

$$\begin{equation} p(\mu \mid k, \hat{\epsilon }_{{\rm meas}}, \hat{\Lambda }) = \frac{ p(\mu ) p(\hat{x} \mid \mu , k, B) }{ p(\hat{x} \mid k, B) }, \end{equation} \tag{ 10 }$$

where the marginal probability for x is

$$\begin{equation} p(x \mid k, B) = \int _0^\infty {\rm d}\mu \, p(\mu ) p(x \mid \mu , k, B). \end{equation} \tag{ 11 }$$

For illustration we take a uniform prior, p(μ) = const, cf equation (13) of [1]. We have

$$\begin{eqnarray} p(x \mid k, B) &= \int _0^\infty {\rm d}\mu \, p_0(x) \frac{ 1 + \mu \epsilon _{{\rm meas}}(x) [\Lambda (x) + 1/k] }{ [1 + \mu \epsilon _{{\rm meas}}(x) / k]^{k+1} } \nonumber \\ &= p_0(x) \frac{k}{k - 1} \frac{1 + \Lambda (x)}{\epsilon _{{\rm meas}}(x)} \end{eqnarray} \tag{ 12 }$$

and therefore the posterior distribution is

$$\begin{equation} p(\mu \mid k, \hat{\epsilon }_{{\rm meas}}, \hat{\Lambda }) = \frac{k - 1}{k} \frac{\hat{\epsilon }_{{\rm meas}}}{1 + \hat{\Lambda }} \frac{ 1 + \mu \hat{\epsilon }_{{\rm meas}}(\hat{\Lambda } + 1/k) }{ (1 + \mu \hat{\epsilon }_{{\rm meas}}/ k)^{k+1} }. \end{equation} \tag{ 13 }$$

Equation (13) is the correct version of equation (24) of [1].

An alternative, but equivalent, method of handling the parameter α is to use method (ii) above, and treat it as a nuisance parameter in a posterior distribution $p(\mu ,\alpha \mid k, \hat{\epsilon }_{{\rm meas}}, \hat{\Lambda })$ and to integrate it out:

$$\begin{equation} p(\mu \mid k, \hat{\epsilon }_{{\rm meas}}, \hat{\Lambda }) = \int _0^\infty {\rm d}\alpha \, p(\mu , \alpha \mid k, \hat{\epsilon }_{{\rm meas}}, \hat{\Lambda }). \end{equation} \tag{ 14 }$$

From Bayes's theorem, we see that the posterior with the nuisance parameter is given by

$$\begin{equation} p(\mu , \alpha \mid k, \hat{\epsilon }_{{\rm meas}}, \hat{\Lambda }) = \frac{ p(\mu ) p(\alpha ; k) p(\hat{x} \mid \mu , \alpha , B) }{ p(\hat{x} \mid k, B) }, \end{equation} \tag{ 15 }$$

where

$$\begin{equation} p(x \mid k, B) = \int _0^\infty {\rm d}\mu \, p(\mu ) \int _0^\infty {\rm d}\alpha \, p(\alpha ; k) p(x \mid \mu , \alpha , B) \end{equation} \tag{ 16 }$$

is exactly the marginal probability for x given earlier, as can be seen by substituting equation (8) into equation (11). Finally, inserting equation (15) into equation (14) gives us

$$\begin{eqnarray} p(\mu \mid k, \hat{\epsilon }_{{\rm meas}}, \hat{\Lambda }) &= \frac{p(\mu )}{p(\hat{x} \mid k, B)} \int _0^\infty {\rm d}\alpha \, p(\alpha ; k) p(\hat{x} \mid \mu , \alpha , B) \nonumber \\ &= \frac{ p(\mu ) p(\hat{x} \mid \mu , k, B) }{ p(\hat{x} \mid k, B) } \end{eqnarray} \tag{ 17 }$$

which is the same as equation (10). As expected we obtain the same final result for the posterior distribution if we consider α to be a nuisance parameter that we integrate out or if we marginalize the likelihood over the parameter α.

We can now explain the mistake in [1]. Note that we had to modify equation (1) to account for the addition of the random variable α, while [1] used its equation (11) without modification. Thus the substitution of equation (11) into equation (23) in [1] is not equivalent to the fundamental definition of the marginalized posterior, given by equation (10) or (14) here.

The same mistake was made in section 5.2 of [1]. Its equation (25) is the incorrect starting point. Instead, if we suppose Λ(x) = βΛ_meas(x) where β is a random variable with expectation value of unity, then it is straightforward to get the correct result following the methods explained here. The easiest approach is to first compute the marginal likelihood, which in this case is obtained by substituting Λ(x) = βΛ_meas(x) in equation (2), multiplying by p(β), and integrating over β. The result has the exact same form with Λ(x) → Λ_meas(x). It then follows that the correct version of equation (27) in [1] is its equation (14) with this same substitution. Thus, the posterior marginalized over uncertainties in $\hat{\Lambda }_{{\rm meas}}$ does not depend on the variance of this uncertainty at all (provided that the expectation value of β is unity), and equations (26) and (28) in [1] are superfluous. This also means marginalizing over uncertainties in $\hat{\Lambda }_{{\rm meas}}$ will not change upper limits on the rate at all. The conclusion of section 5.2 of [1] is still correct: marginalization over uncertainties in $\hat{\Lambda }_{{\rm meas}}$ can be neglected.

We would like to acknowledge many useful discussions with members of the LIGO Scientific Collaboration that took place during the internal collaboration review of [2] which led to the discovery of the error in [1]. This work has been supported in part by NSF grant PHY-0970074. SF would like to thank the Royal Society for support.