# Data analysis strategies for the detection of gravitational waves in non-Gaussian noise

## Abstract

In order to analyze data produced by the kilometer-scale gravitational wave detectors that will begin operation early next century, one needs to develop robust statistical tools capable of extracting weak signals from the detector noise. This noise will likely have non-stationary and non-Gaussian components. To facilitate the construction of robust detection techniques, I present a simple two-component noise model that consists of a background of Gaussian noise as well as stochastic noise bursts. The optimal detection statistic obtained for such a noise model incorporates a natural veto which suppresses spurious events that would be caused by the noise bursts. When two detectors are present, I show that the optimal statistic for the non-Gaussian noise model can be approximated by a simple coincidence detection strategy. For simulated detector noise containing noise bursts, I compare the operating characteristics of (i) a locally optimal detection statistic (which has nearly-optimal behavior for small signal amplitudes) for the non-Gaussian noise model, (ii) a standard coincidence-style detection strategy, and (iii) the optimal statistic for Gaussian noise.

The reliable detection of weak gravitational wave signals from broad-band detector noise from kilometer-scale interferometers such as LIGO  and VIRGO  is the primary concern in developing gravitational wave data analysis strategies. Because of the weakness of the expected gravitational wave signals, it is critical that the detection strategy should be nearly optimal. However, the detector noise may not be purely stationary and Gaussian, so it is important that the detection strategy also be robust so that detections will be reliable.

Until now, most work on the development of data analysis strategies has been limited to the case of stationary Gaussian noise, though there has been some work  on the creation of vetoes that will discriminate between expected gravitational wave signals and non-Gaussian noise bursts. Because the properties of the noise in the LIGO and VIRGO detectors will not be known in advance, it is difficult to assess how well these strategies and vetoes will perform. When real interferometer data from the 40-meter Caltech prototype is used, it is found that additional vetoes are needed to deal with the abundance of false alarms arising from the non-Gaussian noise .

In this paper, I present a simple non-Gaussian noise model, consisting of Poisson-distributed noise bursts, that represents a number of potential non-Gaussian noise sources that may be present in future interferometers. This noise model is less naïve than the usual assumptions of Gaussian noise alone, but is simple enough that many analytical results can be obtained. By using this noise model, the robustness of various detection strategies can be assessed. A general introduction to signal detection in non-Gaussian noise can be found in Ref. . The new result in this paper is the use of the non-Gaussian noise model in examining the performance of simple multi-detector search strategies, which will be important for gravitational wave searches.

I will adopt the following notation in this paper. The detector output is a set of $`N`$ values that are written collectively as a vector $`𝐡`$ in a $`N`$-dimensional vector space $`V`$. This vector can be thought of as having components $`h_j=(𝐡,𝐞_j)`$, $`j[0,N1]`$, which represent a time series of sample measurements made by the detector. Here, $`𝐞_j`$ is the appropriate Cartesian basis on $`V`$. Alternatively, the vector can be expressed as the set of components $`\{\stackrel{~}{h}_{2k},\stackrel{~}{h}_{2k+1}\}`$, $`k[0,N/21]`$—the real and imaginary Fourier transform components of $`\{h_j\}`$. These components are given by $`\stackrel{~}{h}_j=(𝐡,\stackrel{~}{𝐞}_j)`$, where $`\stackrel{~}{𝐞}_{2k}=(2/N)^{1/2}_{j=0}^{N1}𝐞_j\mathrm{cos}(2\pi jk/N)`$ and $`\stackrel{~}{𝐞}_{2k+1}=(2/N)^{1/2}_{j=0}^{N1}𝐞_j\mathrm{sin}(2\pi jk/N)`$, $`k[0,N/21]`$. Thus, the vectors can be treated in either a time representation or a frequency representation.

A natural inner product of two vectors, $`𝐚`$ and $`𝐛`$, in this vector space is defined as $`(𝐚,\mathrm{𝐐𝐛})`$. The kernel $`𝐐`$ is the inverse of the auto-correlation matrix, $`𝐑`$, of the Gaussian component, $`𝐧_\text{G}`$, of the detector noise. Thus, $`𝐑=𝐧_\text{G}𝐧_\text{G}`$, where the angle brackets denote an average over an ensemble of realizations of detector noise. Vector norms are defined in terms of this inner product, i.e., $`𝐚^2=(𝐚,\mathrm{𝐐𝐚})`$, as are unit vectors: $`\widehat{𝐚}=𝐚/𝐚`$.

The detector noise $`𝐧`$ consists of two components: (i) the Gaussian component $`𝐧_\text{G}`$ (that is always present) and (ii) a possible noise burst component $`𝐧_\text{B}`$ that is present with probability $`P_\text{B}`$. The Gaussian component has the probability distribution $`p[𝐧_\text{G}]\mathrm{exp}(𝐧_\text{G}^2/2)`$. If the burst component is randomly distributed in the vector space $`V`$ with normalized measure $`\widehat{D}[𝐧_\text{B}]`$, then the probability distribution for the noise is

$$p[𝐧]e^{𝐧^2/2}+\frac{P_\text{B}}{1P_\text{B}}\widehat{D}[𝐧_\text{B}]e^{𝐧𝐧_\text{B}^2/2}.$$
(1)   
Suppose that the noise bursts are uniformly-distributed in the vector space $`V`$ out to a large radius $`R`$. Such bursts will typically last the entire duration of interest ($`N`$ samples) and fill the entire frequency band. The noise distribution is approximately

$$p[𝐧]e^{𝐧^2/2}(1+ϵe^{𝐧^2/2})$$
(2)   
for $`𝐧<R`$, where $`ϵ2^{N/2}\mathrm{\Gamma }(N/2+1)R^NP_\text{B}/(1P_\text{B})`$ for large $`R`$.

An alternative derivation of Eq. (2) is the following. Suppose the noise burst is simply an additional transient source of Gaussian noise. Then Eq. (2) can be written as

$$p[𝐧]e^{(𝐧,\mathrm{𝐐𝐧})/2}(1+\epsilon e^{(𝐧,[𝐐𝐐^{}]𝐧)/2}),$$
(3)   
where $`\epsilon =(det|𝐑|/det|𝐑^{}|)^{N/2}P_\text{B}/(1P_\text{B})`$, $`𝐑^{}`$ is the auto-correlation matrix for both ambient and the burst components of Gaussian noise together, and the matrix $`𝐐^{}`$ is the inverse of $`𝐑^{}`$. If the typical noise burst is much louder than the ambient Gaussian noise, then $`𝐐𝐐^{}𝐐`$. Therefore, in the case of loud Gaussian noise bursts, Eq. (3) has the same form as Eq. (2).

The likelihood ratio can now be computed for the two alternative hypotheses: that the output $`𝐡=𝐧`$ is noise alone ($`H_0`$), or that the output $`𝐡=A\widehat{𝐮}+𝐧`$ is a signal $`A\widehat{𝐮}`$ of amplitude $`A`$ embedded in noise ($`H_1`$).<sup>*</sup><sup>*</sup>*In this paper, I consider only signals which are completely known (up to their amplitude). The generalization of this case to one in which the signals have an unknown initial phase is straightforward: see, e.g., Ref. . The likelihood ratio $`\mathrm{\Lambda }(A)=p[𝐡H_1]/p[𝐡H_0]`$ is the ratio of the posterior probability of obtaining the observed output $`𝐡`$ given hypothesis $`H_1`$ to the posterior probability of obtaining $`𝐡`$ given hypothesis $`H_0`$. One then finds that the likelihood ratio is

$$\mathrm{\Lambda }(A)=\frac{\mathrm{\Lambda }_\text{G}(A)+\alpha }{1+\alpha }$$
(4)   
where

$$\mathrm{\Lambda }_\text{G}(A)=e^{A(\widehat{𝐮},\mathrm{𝐐𝐡})A^2/2}$$
(5)   
and

$$\alpha =ϵe^{𝐡^2/2}.$$
(6)   
The quantity $`\mathrm{\Lambda }_\text{G}(A)`$ is the likelihood ratio one would obtain if only the Gaussian noise component were present; it is a monotonically increasing function of the matched filter $`(\widehat{𝐮},\mathrm{𝐐𝐡})`$. When a non-Gaussian noise burst component can occur, the likelihood ratio depends on an additional measured quantity: the magnitude of the output vector, $`𝐡`$. In addition to these functions of the detector output, the likelihood ratio also depends on the expected amplitude $`A`$ of the signal and the proportionality factor $`ϵ`$, which encapsulates the probability of a noise burst and the maximum possible amplitude $`R`$ of the burst.

The quantity $`\alpha `$ plays a central role in the modified likelihood ratio of Eq. (4): it acts as a detector of noise bursts, and vetoes events that are more likely due to the burst. The logarithm of $`\alpha `$ is proportional to the amount of power in the detector output. For Gaussian noise alone, the expected value of the power is $`𝐧_\text{G}^2=N`$, and thus $`\mathrm{ln}\alpha \mathrm{ln}P_\text{B}\frac{1}{2}N\mathrm{ln}(R^2/N)`$ for large $`N`$. Thus, for $`R^2>N`$ (as was assumed above), the value of $`\alpha `$ will be small. However, when a noise burst is present, the value of $`\mathrm{ln}\alpha `$ is increased by a typical amount $`\frac{1}{2}R^2`$, so $`\alpha `$ will typically become large. In fact, $`\alpha `$ is the excess-power statistic for detection of arbitrary noise bursts . For small values of $`\alpha `$, the likelihood ratio approaches the usual Gaussian noise likelihood ratio. However, for large values of $`\alpha `$, the likelihood ratio approaches unity.

The likelihood ratio is a function of the detector output via the two quantities $`x=(\widehat{𝐮},\mathrm{𝐐𝐡})=𝐡\mathrm{cos}\theta `$ and $`𝐡`$. The likelihood ratio also depends on the expected signal amplitude $`A`$ and the factor $`ϵ`$. Thus, the likelihood ratio is

$$\mathrm{\Lambda }(A)=\frac{e^{AxA^2/2}+ϵe^{(x\mathrm{sec}\theta )^2/2}}{1+ϵe^{(x\mathrm{sec}\theta )^2/2}}.$$
(7)   
Figure 1 shows the likelihood ratios as functions of the matched filter statistic $`x`$ and the angle $`\theta `$ for a given value of $`ϵ`$ and $`A`$. Notice that the likelihood ratio is attenuated when the magnitude of the output, $`𝐡`$, is much larger than the largest expected signal; this attenuation occurs at smaller signal-to-noise ratios for smaller absolute values of the direction cosine $`\mathrm{cos}\theta `$.

Because of the non-trivial dependence of the likelihood ratio on the expected signal amplitude, the optimal statistic for signals of some amplitude $`A`$ will not be the same statistic as the optimal statistic for a different amplitude $`A^{}`$. This situation is unlike the case of purely Gaussian noise, for which the likelihood ratio grows monotonically with the matched filter output, regardless of the expected signal amplitude. Because the amplitude of a signal will not typically be known in advance, it is useful to consider a *locally optimal* statistic , which provides the optimal performance in the limit of small amplitude signals. The rationale behind such a choice is that large amplitude signals pose no challenge for detection, so any reasonable statistic will suffice for them. The locally optimal statistic, which is defined as $`\lambda =d\mathrm{ln}\mathrm{\Lambda }/dA|_{A=0}=(\widehat{𝐮},\mathbf{}\mathrm{ln}p[𝐡])`$, is given by

$$\lambda =\frac{x}{1+ϵe^{(x\mathrm{sec}\theta )^2/2}},$$
(8)   
where $`x=(\widehat{𝐮},\mathrm{𝐐𝐡})`$ is the matched filter. Notice that the locally optimal statistic also incorporates a veto based on the value of $`\alpha =ϵ\mathrm{exp}[(x\mathrm{sec}\theta )^2/2]`$. The locally optimal statistic grows approximately linearly for $`x<(2\mathrm{ln}ϵ)^{1/2}\mathrm{cos}\theta `$, at which point it is effectively cut off. This is the same general behavior as was seen before for the likelihood ratio $`\mathrm{\Lambda }(A)`$, but there is no longer any scaling with a prior amplitude.

Another approach is to use a prior distribution $`f(A)`$ for the amplitude, if such a distribution is known, to construct an integrated likelihood ratio $`\mathrm{\Lambda }=\mathrm{\Lambda }(A)f(A)𝑑A`$; this integrated likelihood ratio would then be the optimal statistic for detection of all potential signals. Even if the prior distribution is not known, it is possible to obtain an approximate likelihood ratio if $`f(A)`$ is a slowly-varying function. One can show that $`\mathrm{\Lambda }_\text{G}=\mathrm{\Lambda }_\text{G}(A)f(A)𝑑Af(x)\mathrm{exp}(x^2/2)`$. A further approximation neglects the slowly-varying function entirely; one then obtains $`\mathrm{\Lambda }_\text{G}\mathrm{exp}(x^2/2)`$. Thus,

$$\mathrm{\Lambda }\frac{e^{x^2/2}+ϵe^{(x\mathrm{sec}\theta )^2/2}}{1+ϵe^{(x\mathrm{sec}\theta )^2/2}}.$$
(9)   
It is important to generalize the above analysis to the case in which there are multiple detectors. Such a situation is deemed to be essential for gravitational wave searches since a coincident detection is an important corroboration that a signal is of an astronomical origin rather than internal to the detector. As a simple model,A more general model is being considered by Finn . consider two gravitational wave detectors with identical but independent noise (including possible noise bursts) that also have the same response to gravitational wave signals. If the noise was purely Gaussian, the likelihood ratio for the two detectors together would be $`\mathrm{\Lambda }_\text{G}^{\text{AB}}=\mathrm{\Lambda }_\text{G}^\text{A}\times \mathrm{\Lambda }_\text{G}^\text{B}`$, i.e., the product of the likelihood ratios for Gaussian noise in each separate detector. The resulting likelihood ratio is thus an monotonically increasing function of the *sum* of the matched filter values in each detector: $`\mathrm{ln}\mathrm{\Lambda }_\text{G}^{\text{AB}}x^\text{A}+x^\text{B}`$ where $`x^\text{A}=(\widehat{𝐮},𝐐^\text{A}𝐡^\text{A})`$ is the signal-to-noise ratio in the first detector (A) and $`x^\text{B}=(\widehat{𝐮},𝐐^\text{B}𝐡^\text{B})`$ is the signal-to-noise ratio in the second detector (B). If the likelihood ratios are integrated over a slowly-varying prior distribution of possible signal amplitudes, then one finds $`\mathrm{ln}\mathrm{\Lambda }_\text{G}^{\text{AB}}(x^\text{A}+x^\text{B})^2`$, which also depends on the sum of the two detectors’ matched filter values. (In this case, it is important to integrate over the distribution of amplitudes *after* multiplying the two detector’s likelihood ratios together since the same signal amplitude should be present in each detector.) Thus, one would decide that a signal was present if the sum of the matched filter values in the two detectors exceeded some threshold. An alternative strategy would decide that a signal was present only if *each* of the matched filter values exceeded some threshold, i.e., one thresholds on the quantity $`\mathrm{min}(x^\text{A},x^\text{B})`$ rather than on $`x^\text{A}+x^\text{B}`$.

Suppose the likelihood ratios for both detectors, $`\mathrm{\Lambda }^\text{A}(x^\text{A},\theta ^\text{A})`$ and $`\mathrm{\Lambda }^\text{B}(x^\text{B},\theta ^\text{B})`$, are given by Eq. (7). Then, the likelihood ratio for the combined detector is given by the product of $`\mathrm{\Lambda }^\text{A}`$ and $`\mathrm{\Lambda }^\text{B}`$. Figure 2 shows a plot of contours of constant likelihood ratio as a function of the two measured signal-to-noise ratios for $`\theta ^\text{A}=\theta ^\text{B}=0`$ and $`\theta ^\text{A}=0`$, $`\theta ^\text{B}=\pi /6`$. At low signal-to-noise ratios, the contours are approximately lines of constant $`x^\text{A}+x^\text{B}`$—in this regime, the likelihood ratio behaves like the Gaussian likelihood ratio. However, for large values of either $`x^\text{A}`$ or $`x^\text{B}`$, the contours are approximately lines of constant $`\mathrm{min}(x^\text{A},x^\text{B})`$.

The locally optimal two-detector statistic is $`\lambda ^{\text{AB}}=\lambda ^\text{A}+\lambda ^\text{B}`$, where $`\lambda ^\text{A}`$ and $`\lambda ^\text{B}`$ both have the form of Eq. (8). This statistic also depends on $`x^\text{A}`$, $`x^\text{B}`$, $`\theta ^\text{A}`$, and $`\theta ^\text{B}`$. As before, the contours of constant $`\lambda ^{\text{AB}}`$ for fixed $`\theta ^\text{A}`$ and $`\theta ^\text{B}`$ are lines of constant $`x^\text{A}+x^\text{B}`$ for low signal-to-noise ratios but are lines of constant $`\mathrm{min}(x^\text{A},x^\text{B})`$ if the signal-to-noise ratio in one of the detectors is large.

These results somewhat justify the use of the minimum statistic $`\mathrm{min}(x^\text{A},x^\text{B})`$ when there are two detectors with independent noise operating. When one of the signal-to-noise ratios is large but the other is moderate, the minimum statistic will give approximately the same value as the locally optimal statistic. In this case, the effect of the $`\alpha `$-based veto in the locally optimal statistic is mimicked in the minimum statistic. The locally optimal statistic is better able to deal with the case when there are noise bursts present in both detectors, but this case occurs with probability $`P_\text{B}^2`$, so it should be extremely rare. Moreover, the minimum statistic is less powerful than the locally optimal statistic for small signal-to-noise ratios in both detectors. In this case, the locally optimal statistic is approximately $`x^\text{A}+x^\text{B}`$, which is the same as the optimal statistic for Gaussian noise. However, the optimal statistic for Gaussian noise may be unsuitable for the noise model I have considered because the false alarm probability cannot be reduced below $`P_\text{B}`$ for reasonable thresholds. (The minimum statistic suffers the same problem but around the much lower probability $`P_\text{B}^2`$.)

To see the relative performance of these statistics, it is useful to consider the operational characteristics of the tests. For some fixed false alarm probability $`Q_0`$, the operational characteristic is the detection probability $`Q_1`$ as a function of signal amplitude. I have computed these curves using Monte Carlo techniques. The noise model I used corresponds to that of Eq. (2) with burst probability $`P_\text{B}=1\%`$ and maximum burst amplitude $`R=25`$. I fixed the false alarm probability to $`Q_0=10^3`$ and I examined a vector space dimension of $`N=4`$. Figure 3 shows the relative performances of the three statistics in terms of the detection probability as a function of signal amplitude. One sees that the locally optimal statistic performs the best for small amplitude signals (as expected), but that it has poor performance for large amplitude signals. This is because the large amplitude deviations are interpreted as noise bursts and are suppressed. (If the locally optimal statistic were used in a real search, these large events would not be rejected outright but would rather be subjected to further scrutiny; thus the attenuation of the locally optimal statistic for large signal amplitudes is somewhat misleading.) The sum of the signal-to-noise ratio, which is optimal statistic for Gaussian noise, is seen to have poor performance. The reason is that because the false alarm probability is much smaller than the burst probability, the threshold required for this statistic becomes unreasonably large. However, since the false alarm probability is much higher than the burst probability *squared*, the minimum statistic performs reasonably well for both small and large amplitude signals.

The conclusion that can be drawn from this analysis is that search strategies based on a Gaussian noise model can be significantly improved by considering a more realistic non-Gaussian noise model, which may contain noise bursts that occur relatively frequently. If a method can be devised to detect these noise bursts and veto data that contains a noise burst, then the use of such a veto effectively reproduces the locally optimal strategy for the noise model containing bursts. However, the necessary veto may not be easily found since it may not be possible to determine the non-Gaussian properties of the detector noise to sufficient accuracy. A viable alternative when two detectors are operating is to use a coincidence strategy of detection. Since the minimum statistic used in the coincidence strategy is a robust statistic, one would expect that it should have good performance for any realistic noise model.

I would like to thank Bruce Allen, Patrick Brady, Éanna Flanagan, and Kip Thorne for their many useful comments on this paper. This work was supported by the National Science Foundation grant PHY-9424337.