A Bayesian approach to the follow-up of candidate gravitational wave signals
aa r X i v : . [ g r- q c ] J u l A Bayesian approach to the follow-up of candidate gravitational wave signals
J. Veitch , A. Vecchio , School of Physics and Astronomy, University of Birmingham, Edgbaston, Birmingham B15 2TT, UK Northwestern University, Department of Physics and Astronomy,2145 Sheridan Road, Evanston, IL 60208, USA
Ground-based gravitational wave laser interferometers (LIGO, GEO-600, Virgo and Tama-300)have now reached high sensitivity and duty cycle. We present a Bayesian evidence-based approach tothe search for gravitational waves, in particular aimed at the follow-up of candidate events generatedby the analysis pipeline. We introduce and demonstrate an efficient method to compute the evidenceand odds ratio between different models, and illustrate this approach using the specific case of thegravitational wave signal generated during the inspiral phase of binary systems, modelled at theleading quadrupole Newtonian order, in synthetic noise. We show that the method is effective indetecting signals at the detection threshold and it is robust against (some types of) instrumentalartefacts. The computational efficiency of this method makes it scalable to the analysis of all thetriggers generated by the analysis pipelines to search for coalescing binaries in surveys with ground-based interferometers, and to a whole variety of signal waveforms, characterised by a larger numberof parameters.
PACS numbers: 04.80.Nn, 02.70.Uu, 02.70.Rr
I. INTRODUCTION
Ground-based gravitational wave (GW) laser interfer-ometers – LIGO [1], Virgo [2], GEO-600 [3] and TAMA-300 [4] – have been in operation for a few years, alter-nating times of data taking for science analysis at pro-gressively greater sensitivity and duty cycle, with com-missioning periods to improve the instruments’ perfor-mance. Recently LIGO has completed its fifth sciencerun (S5) which lasted for ≈ ≈ ≈
10 across thewhole frequency band and shift the low-frequency cut-offto about 10 Hz.The search for GWs has therefore reached a stage inwhich the first direct detection is plausible and there areexpectations that by the time the instruments operatein advanced configuration it will be possible to routinelyobserve a wide variety of sources in this new observationalwindow and study them in detail, see [9] and referencestherein.One of the most promising classes of sources are binarysystems of compact objects, which are reviewed in [9].GW laser interferometers are omni-directional detectorsthat continuously monitor the whole sky looking for rareand/or weak signals. The general approach to the anal-ysis of the data is to employ efficient algorithms able tokeep up with the data flow and to identify events above a given signal-to-noise ratio threshold. The candidateevents (or “triggers”) are then followed up with a rangeof techniques to decide whether or not a GW signal ispresent.Most of the effort on the data analysis side has so farbeen devoted to the development and implementation oftechniques for the mass analysis of the data. These algo-rithms have been reaching maturity and are being appliedto an increasingly large volume of data [10, 11, 12, 13, 14];however comparatively less experience has been gainedon follow-up methods, see e.g.
Section VIII of Ref [10](and references therein) and Ref [15]. So far the approachto candidate follow-ups is based [10, 15, 16, 17, 18] on(i) addressing the probability of false alarm of an eventagainst the background, usually quantified by repeatingthe analysis on data from multiple interferometers shiftedin time by an unphysical amount(so-called ”time-slides”),(ii) looking for possible correlations between monitor-ing channels and the GW channel that could reveal ananomalous behaviour of the instrument (or sub-system)or anthropogenic causes around the time of the candidateand (iii) checking that the signal at different interferome-ters shows consistent behaviour. In this paper we proposea conceptually and practically very different approach tothe follow-up of candidates: Bayesian model selection (orhypothesis testing) that is based on the evaluation of the marginal likelihood or evidence of a specific model andon the odds ratios between competing hypotheses. Notethat in principle the method discussed here could alsobe used to search for signals in the whole data set; how-ever the computational costs involved in adopting suchan approach and the performance of the existing analy-sis algorithms do not strongly support an effort in thisdirection at present.The key issue that we address in this paper is, givena data set and a set of prior information, how we cal-culate the probability of a GW signal being present.In the formalism of Bayesian inference this is trans-lated into considering two models – Model 1: “there isa GW signal and noise”, and Model 2: “only noise ispresent” – and computing the probability of each of them.Bayesian inference provides a conceptually straightfor-ward prescription to evaluate the probabilities, or theodds ratio of the two models. In the context of GWobservations, Bayesian inference is starting to be con-sidered [19, 20, 21, 22, 23] as a powerful approach forground-based observations. The heavy computationalload involved in this method when using exhaustive inte-gration has however limited its use on real data. Recentlyapplications of reversible jumps Markov-Chain Monte-Carlo methods (MCMC) [24] have been considered totackle this issue in the context of searches for binary in-spirals [22]: in this approach Markov-chains are free tomove between models, and therefore one can estimate theBayes factor from the relative time spent by the chainsin each one. This technique is promising, still the com-putational burden is quite significant. In this paper weconsider a different and efficient numerical implementa-tion of the direct calculation of odds ratios that makesthis approach realistically applicable to extensive follow-up studies of triggers and that we show is robust againsta selection of noise artefacts.The paper is organised in the following way: in Sec-tion II we introduce the key concepts about model selec-tion and hypothesis testing in the Bayesian framework;we also introduce the simple case – possible observationsof inspiralling binary systems in a single interferometer– that we will consider in the paper; Section III con-tains the key results of this paper, where we show thepower and robustness of Bayesian hypothesis testing forfollow-up studies; Section IV contains our conclusionsand pointers to future work. II. MODEL SELECTIONA. Overview
Let us consider a set of hypotheses (or models) { H } and denote with I all relevant prior information we hold.The predictions made by a model depend on a set of0 or more parameters, and the possible combinations ofparameters define the parameter space of that model. Aparameter vector ~θ represents a point in parameter space,and although each model has its own set of parameters,for ease of notation we shall use ~θ for them all, and werepresent the data under consideration by d .A straightforward application of Bayes theorem yieldsthe probability of a given model: P ( H j | d, I ) = P ( H j | I ) P ( d | H j , I ) P ( d | I ) . (1)If the models under consideration are exhaustive and mu-tually exclusive, then the probabilities above are clearly normalised and satisfy: P Nj =1 P ( H j | d, I ) = 1. In Equa-tion 1, P ( H j | d, I ) is the posterior probability of model H j given the data, P ( H j | I ) is the prior probability of hy-pothesis H j and P ( d | H j , I ) is the marginal likelihood or evidence for H j that can be written as: P ( d | H j , I ) = L ( H j )= Z d~θ p ( ~θ | H j , I ) p ( d | ~θ, H j , I ) , (2)where p ( ~θ | H j , I ) is the prior probability density distri-bution for the parameter vector ~θ and is normalised to1. p ( d | ~θ, H j , I ) is called the likelihood , and is an un-normalised measure of the fit of the data to the model.If we had an exhaustive set of models, we could sim-ply calculate the probability of each model and comparethem to find the most likely. Unfortunately we do nothave a complete set of models of the data from a GWdetector, but we still want to compare the models we do have. This is a normal procedure in Bayesian inference,which gives what is termed the odds ratio between hy-potheses, which is a quantification of their relative prob-ability.For example, the odds of model H j against model H k is O jk = P ( H j | d, I ) P ( H k | d, I ) = P ( H j | I ) P ( H k | I ) P ( d | H j , I ) P ( d | H k , I )= P ( H j | I ) P ( H k | I ) B jk (3)where P ( H j | I ) /P ( H k | I ) encodes the ratio of theprior state of belief of the two models and B jk ≡ P ( d | H j , I ) /P ( d | H k , I ) is known as the Bayes factor ; themarginal likelihoods P ( d | H j,k , I ) are given by Eq. (2).The practical advantage of considering the odds ratio (3)is the fact that one does not need to evaluate P ( d | I ). Ofcourse, if one wanted to evaluate the odds of model H j with respect to all the other independent alternatives,the appropriate quantity is:O j, other = P ( H j | d, I ) P k = j P ( H k | d, I ) = 1 P j = k O kj . (4)An interesting consequence of using Bayesian model se-lection is that it automatically includes a quantitativeversion of “Occam’s Razor”, the principle that the sim-pler model should be prefered. Because the evidence isfound by integrating over the entire prior domain of amodel, and this prior is normalised to unity, the largerthe volume of parameter space which the model spreadsits prior over the lower the resulting evidence will be,all else being equal. If two models can assign the samemaximum likelihood to some data by fitting it equallywell, the one which makes the more precise prediction ofthe parameters which produce the maximum likelihoodwill benefit the more than one which makes the broaderprediction. [30, 37] B. A concrete example
The specific problem that we are addressing here ishow to answer the question: what are the odds that thereis a signal present in the given observations ? We actuallyneed to be more specific and spell out the whole set ofbackground information I for the problem at hand. Herewe concentrate on searches for inspiralling compact ob-jects, though the method outlined here has the potentialof much wider applications. The background information I is therefore as follows: (i) the data set d consist of thesuperposition of noise n and (possibly) an inspiral GWsignal h ; (ii) if present, only one GW signal is present atany one time; (iii) the waveform model is exactly known,though the actual parameters characterising the sourceare unknown, within some prior range; (iv) gravitationalradiation and noise are statistically independent; (v) thenoise is a Gaussian and stationary random process withzero mean and variance (at any frequency) described by a known spectral density ; (vi) observations are carried outwith a single instrument.We are therefore considering a situation in which(schematically) there are only two models: • The first hypothesis – that we will label as H S , fornoise and GW signal – corresponds to “there is asignal and noise present”; the data set in the timedomain is therefore described by d ( t ) = n ( t ) + h ( t ) . (5) • The second hypothesis – that we will indicate as H N , for noise only – corresponds to ”there is onlynoise present”, and the data set is therefore: d ( t ) = n ( t ) . (6)The goal of the analysis is therefore to compute the oddsratio, Equation (3), between H S and H N that we willindicated with O SN .In this paper, for simplicity we consider the obser-vations of gravitational radiation produced during theinspiral phase of a binary system of non-spinning com-pact objects. The two objects have individual mass m , , and the binary’s chirp and total mass are there-fore M = ( m m ) / ( m + m ) − / and m = m + m ,respectively. We have also assumed m = m . We modelradiation at the leading Newtonian quadrupole order. Asthe analysis is more conveniently carried out in the fre-quency domain, the model that we adopt is the station-ary phase approximation to the radiation in the Fourierdomain, see e.g. [31]. Describing with ˜ g ( f ) the Fouriertransform at frequency f of the time-domain function g ( t ) we can express the GW signal as:˜ h ( f ; ~θ ) = (cid:26) A M / f − / e iψ ( f ; ~θ ) ( f ≤ f LSO )0 ( f > f
LSO ) , (7) where ψ ( f ; ~θ ) = 2 πf t c − φ c − π π M f ) − / (8)is the signal phase in the frequency domain and f LSO isthe frequency of the last stable orbit; in this paper weset f LSO equal to the last stable circular orbit for theSchwarzschild space-time and equal mass non-spinningobjects, so that f LSO = (6 / / π M ) − . (9)In Eqs (7) and (8), ~θ is the 4-dimensional parameter vec-tor describing the GW signal; as parameters we choose ~θ = { A, M , t c , φ c } , (10)where A is an overall amplitude that depends on thesource position in the sky and the inclination and po-larisation angle of the source, and t c and φ c are thetime and phase at coalescence. Note that we are us-ing geometrical units in which c = G = 1 and therefore M ⊙ = 4 . × − s.We model the noise as a Gaussian and stationary ran-dom process with zero mean, and variance at frequency f k given by σ k = (cid:16) f k f (cid:17) − + 1 + (cid:16) f k f (cid:17) , which roughlyrepresents the shape of a first generation interferometernoise curve up to a multiplicative constant of the orderof 10 − , which is irrelevant as it cancels in the likelihoodfunction; f = 150 Hz is chosen to pick the frequency ofmaximum sensitivity.In a real application, the variance of each point canbe estimated from the one-sided power spectral density,calculated with Welch’s method (for example) from thedata around the segment of interest [38]. Notice that inthe next section we will consider deviations of the actualnoise affecting the measurements (but the noise modelused to construct the likelihood will remain unchanged)from the Gaussian and stationary noise shown above. Wewill discuss this in more detail in Section III. C. Evaluation of the odds ratio
We can now spell out the odds ratio that needs to becalculated; from Equation (3), this is given byO SN = P ( H S | I ) P ( H N | I ) P ( d | H S , I ) P ( d | H N , I ) . (11)For the noise model in this case, there are no free param-eters since the noise profile is known, and the evidence isgiven for Gaussian noise by, P ( d | H N , I ) = N Y k =1 (cid:0) πσ k (cid:1) − exp − (cid:12)(cid:12)(cid:12) ˜ d k (cid:12)(cid:12)(cid:12) σ k , (12)where the index k is a short-hand representation for thosequantities dependent on frequency f k , and N is the totalnumber of data points. Each point k has a real andimaginary observation, the variance of each part will beequal in any realistic data, σ k = σ k R = σ k I , although inprinciple they could differ. Evaluating the evidence forthis case requires no integration.For the signal model however, the evaluation requiresintegration over the prior domain (denoted Θ ) of all pa-rameters, codified in a vector ~θ and is given by P ( d | H S , I ) = Z Θ d~θ p ( ~θ | H S , I ) N Y k =1 (cid:0) πσ k (cid:1) − × exp − (cid:12)(cid:12)(cid:12) ˜ h k ( ~θ ) − ˜ d k (cid:12)(cid:12)(cid:12) σ k . (13)The prior probability density function p ( ~θ | H S , I ) is forthis proof of concept case uniform on all four parameters,within the domain defined below. It should be noted thatthe prior must be normalised to one when computing theBayes factor, and that increasing the range of the priorwill decrease the Bayes factor, as the model becomes lesspredictive and is penalised in automatic accordance withOccam’s razor.In the case of the Newtonian inspiral model, the com-putation of the marginal likelihood (13) requires the eval-uation of a 4-dimensional integral which could be calcu-lated using a grid-based approach. However we are in-terested in developing a general and flexible approach,and for inspiral binaries we will need to expand the evi-dence calculation to include more realistic waveforms forthe analysis of real data with post-Newtonian effects andspins, which have many more parameters: the integral inEq. (13) rapidly becomes unfeasible to integrate exhaus-tively. To avoid this problem, we have used a probabilis-tic algorithm called Nested Sampling [25, 26, 27], whichhas been used before in the context of Bayesian hypoth-esis testing in cosmology [28], but not considered so farin gravitational wave data analysis. The application ofthis technique has required significant development andtuning of the algorithm; the details of this work are be-yond the scope of this paper and will be documented ina separate publication [29]. Here we shall focus on theadvantages of using the evidence in the context of inspi-ral analyses, which in the next section will be presentedwith discussion of results from some fiducial data analysisproblems. III. RESULTS
In this Section, we will document the results of fourparticular problems of interest on which the algorithmhas been tested: (i) A test on pure Gaussian noise, (ii)the same noise with added signals of different strengths introduced in order to explore the sensitivity of the al-gorithm to different signal-to-noise ratios. In addition tothese cases of obvious relevance, we have also consideredthe situation where some noise artefacts – not includedinto the model for the computation of the marginal likeli-hood and the Bayes factor – were added to the data (thatis to the simple Gaussian and stationary noise contribu-tion). In particular, (iii) a decaying sinusoid waveformwhich might approximate an instrumental ringdown, and(iv) a Poissonian noise component. The two latter casesin particular are intended to look at the situation wherethe models to be tested do not fit the data well at all, inother words they are not exhaustive. This is of interestbecause in reality there will be more types of interferencewith the data than we could possibly hope to model, sorobustness against such artefacts is an essential charac-teristic of a search algorithm. This is directly related tothe false-alarm rate of existing searches.In this paper, we have used synthetic data sets in thefrequency band 40 - 500 Hz, and used 30 s of data. Thepriors defining the range of the model were such that A ∈ [0 , × ], M ∈ [7 . , . M ⊙ , t c ∈ [19 . , . s , φ c ∈ [0 , π ). The large amplitude reflects the conversionof mass into seconds in Equation (7). The injected value M = 8 . M ⊙ was chosen purely to speed up the analysis,as its highest frequency f LSO ( M = 8 . M ⊙ ) = 478 .
45 Hzallowed frequency components above this value to beeliminated from the innermost loop of the calculation.These prior ranges were chosen to approximate the size ofthe parameter space that would have to be searched in areal application, which we suggest would be run on candi-dates triggered by a higher stage in a pipeline, providingsome information about the possible characteristics of thesignal, namely M and t c . On these data sets and withsuch prior ranges, the calculation of the evidence tookbetween several minutes and several hours on a single 2GHz CPU; the efficiency of the algorithm varies howeverwith the specific tuning in a complicated fashion whichwill be described in a future paper [29]. The speed of thecomputation makes this approach amenable to extensionto waveforms characterised by a much larger number ofparameters, and to full lists of triggers generated by aninspiral search pipelines, see e.g. Refs. [10, 11, 12, 13].
A. Sources of Uncertainties
Before presenting the results of the example trials, weshall emphasise that there are two contributions to thedistribution of results which are obtained when using thismethod on any data set. These are the inherent uncer-tainties from the probabilistic nature of the algorithmitself, and the variations in results due to different noiserealisations, produced by the random nature of the noise.The contribution from the noise is an inherent partof any search and the results will naturally vary fromdataset to dataset depending on the particular observa-tions made. This cannot be averted, nor should it be.The uncertainty inherent in the algorithm itself however,is a true “error” in the result, since there is some par-ticular number which represents the odds ratio, and anydeviation from this is undesired. Fortunately, by increas-ing the running time of the algorithm, these errors can beminimised to a level similar to, or below the variation dueto different noise realisations. The inherent uncertaintycan in principle be made arbitrarily small, but in the in-terest of keeping a low run-time, in this paper we havechosen a level roughly similar to the changes caused bythe random fluctuations of the noise; this is appropriatefor the issues discussed here, and none of our conclusionsare affected by the approximation in the evaluation ofthe integral in Eq. (13).We have run two simple test cases to illustrate the dif-ferent sources of uncertainties that affect the results; inboth cases we have run the algorithm a number of timeson signal-free data sets; in the first instance, data setswere generated to contain different realisations of Gaus-sian and stationary noise; in the second, the evidenceevaluation was performed multiple times on the same noise realisation; Figures 1 and 2 summarise the results.Figure 1 shows the distribution of log Bayes factorsrecovered from the analysis of 100 different noise realisa-tions, and represents the background level of changes inthe odds which is caused by random fluctuations in thenoise. Figure 2 shows the inherent uncertainty causedby the algorithm when repeatedly running on a singlenoise realisation of the same characteristics. The distri-bution is shown for 200 runs, and the accuracy level hasbeen chosen such that the errors from this effect are notgreater then those in Figure 1.
B. Stationary Gaussian noise with no signal
The most basic test of whether the evidence-based ap-proach can be helpful in a follow-up analysis is to findhow it responds in the absence of a signal. At the presenttime, this is still the routine situation in gravitationalwave data analysis. The results shown in Figure 2 arethe values of the Bayes factor, cf. Eqs (3) and (11), be-tween the signal hypothesis and the noise hypothesis, ascalculated by performing 200 runs on identical noise re-alisations.The recovered Bayes factors for the hypotheses favour,as expected, the noise only hypothesis, but are (rela-tively) close to unity. This number is then multiplied bya prior odds ratio, see Eq. (11), to get a total posteriorodds ratio. A reasonable prior would give much largercredence to the hypothesis that no signal is present (i.e.have a value ≪ −1.9 −1.8 −1.7 −1.6 −1.5 −1.40510152025 log Bayes factor C oun t s FIG. 1: The distribution of log Bayes factors recovered fromrunning the algorithm on 100 different Gaussian and station-ary noise realisations, each of length 13 800 complex samplesas described in the text. The distribution in values arisesfrom random fluctuations of the noise, which is contrastedwith Figure 2, where the inherent uncertainty in the algo-rithm is shown. In both cases, the variation in Bayes factorsfound in datasets with no signal is extremely small comparedto the change in Bayes factor that the presence of a signalproduces, which is shown in Figure 3 for a range of signal-to-noise ratios. point of view, so as to obtain a desired “false alarm rate”,by performing simulations where known GW signals areadded to the data either in hardware or in software.In the next Section we test the response of the al-gorithm to a GW inspiral signal added to noise, usingsignal-to-noise ratios ranging from 1 to 7. From this weshall see that the evidence is an extremely sensitive in-dicator of the presence of a signal, indeed it is provablythe optimal inference for model comparison [30].
C. Signal injected into stationary Gaussian noise
In this section we discuss the results obtained byadding GW inspiral signals, with signal-to-noise ratio(SNR) varying between 1 and 7, to stationary and Gaus-sian noise, of the same statistical properties as in sectionIII B. The three non-amplitude parameters were keptthe same for all of the injections, and had the values M = 8 . ⊙ , t c = 20 . φ c = 0 .
0. The algorithm wasrun on each of these datasets 20 times, and the resultsare summarised in Figure 3. The optimal signal-to-noiseratio [31] is given by −1.75 −1.7 −1.65 −1.6 −1.5505101520253035 log Bayes factor C oun t s FIG. 2: The distribution of log Bayes factors recovered fromrunning the algorithm 200 times on the same dataset whichcontained no signal. The noise is modelled as a Gaussianand stationary random process. The distribution of resultsis caused by the probabilistic nature of the algorithm, andthe range of this distribution can be reduced to have an arbi-trarily small range at the expense of increased computationtime. In this example and throughout the rest of the paper,the distribution width is chosen to be the same order as theuncertainty caused by different noise realisations, shown inFigure 1.
SNR = vuut N X k =1 | ˜ h k | σ k . (14)It is very clear from Figure 3, that the Bayes factor B S N rapidly increases with SNR, hence the logarithmicscale on the vertical axis. This shows that the odds ratiobetween the two candidate models in this case climbsvery rapidly to favour the signal model H S once sufficientevidence is present.We can now see the effect that the prior odds ratiohas on the posterior odds ratio: it acts as a thresholdvalue, above which the overall odds are in favour of thesignal model. The value of the prior odds effectively setsa limit on SNR above which the odds favour the signalmodel, but once this threshold has been reached, the in-creasing Bayes factor will rapidly climb by many ordersof magnitude.In the cases discussed above, the SNR is chosen bychanging the overall amplitude of the injection, whileholding the injected chirp mass constant between eachcase. It should be noted that there is an additionalchange in the odds ratio which comes about when theinjected mass changes. This is due to the fact that thenoise evidence P ( d | H N , I ) is dependent on the specific l og B a y e s f ac t o r FIG. 3: The mean recovered log Bayes factors from addingan inspiral signal of increasing signal-to-noise ratio, shownon the x -axis, to Gaussian and stationary noise. The algo-rithm was run on each dataset (identical noise realisation) 20times: each point on the plot corresponds to the mean of theBayes factor. Note that the Bayes factor grows exponentiallywith increasing SNR, such that it is an extremely sensitiveindication of the presence of a signal. The signal and dataparameters were chosen as described in the text. shape of the waveform that is present in the data andits inner product with the individual noise realisation.Writing the evidence for the noise model when the datais composed of noise and signal, ˜ d k = ˜ n k + ˜ h k , shows theorigin of the effect: P ( d | H N , I ) ∝ Y k exp " − | ˜ n k | σ k + | ˜ h k | σ k + ˜ n ∗ k ˜ h k + ˜ n k ˜ h ∗ k σ k ! (15)The first term in bracket is a measure of the fit of thenoise to the estimated background noise spectrum { σ k } independent of the signal. The second term is the signal-to-noise ratio squared, which is independent of the noiserealisation. The third term is the important one in thiseffect: it measures the contrast between the signal wave-form and the noise realisation. It is possible to have aconstant SNR while the Bayes factor changes if this termvaries, although because the noise and signal are uncorre-lated ( h ˜ n ∗ k ˜ h k i = h ˜ n k ˜ h ∗ k i = 0), it should be small in com-parison to the other terms in this equation. This termcan change through either differing noise realisations, ora change in the waveform’s shape , the phase, or the over-all amplitude of the signal. Since we have used the samevalue of M in these simulations, this effect should befurther reduced.The choice of the prior odds is an interesting questionwhich we do not attempt to answer conclusively here,but at least two possibilities seem reasonable. Since theprior hypothesis probability P ( H S | I ) on the signal modelreflects our knowledge before we examine the dataset, wecould consider our knowledge of the rate of inspiral eventswhich would be visible to the interferometer in question(see e.g. [32] and references therein), and the length ofthe prior window on t c . Multiplying these would give usan estimate of the number of inspirals we would expectto see in that time period and mass range, which couldbe used as the prior odds on seeing an inspiral event inthat time.Alternatively, one could perform a large number of tri-als on different data realisations and find the distributionof Bayes factors for “false alarm” signals with their fre-quency. A desired false alarm rate could then be setby choosing the prior odds ratio to be the inverse of theBayes factor for that false alarm rate, see e.g. [39]. In realdata where it is not known a-priori if a signal is present,were the search algorithm extended to coherent multipleinterferometer models, this could easily be accomplishedby offsetting the data from two or more interferometers,so that any real signal would not appear with phase co-herence. This is similar to the approach taken in theexisting inspiral analysis pipeline to perform backgroundrate estimation, see e.g. Refs. [10, 11, 12, 13, 14]. Thischoice also has appeal as it would automatically includethe effect of noise artefacts present in real data whichmight cause the odds ratio to increase even in the casewhere there is no astrophysical signal, partially matchingthe inspiral model.A more Bayesian treatment of spurious artefacts mightattempt to model them too, and then compare theirmodel evidences to those for the GW signal and noisemodels to test each hypothesis against the others. How-ever, it is unlikely to be possible to model every singlenoise artefact that appears in a real interferometer, so itis desirable that the basic algorithm functions well evenin the presence of such disturbances. In the followingsections we report results of tests of the performance ofthe algorithm in the presence of some typical (thoughsimulated here) deviations of the noise from the simpleGaussian and stationary behaviour. As we will see, thealgorithm is less sensitive to these noise artefacts thanthe true signal, though the Bayes factor decreases as itshould. This result indeed supports the usefulness ofevidence-based algorithms as part of the tools for follow-up analyses in searches for GWs in the data of currentinterferometers.
D. Stationary Gaussian noise with ringdowninjected
One type of artefact which is often encountered in realdata is the instrumental ringdown, where some compo-nent of the detector is excited and produces dampedoscillations in the gravitational wave channel, gradually decaying. We have modelled the resulting signal for thepurposes of adding it to the synthetic data set as a genericdecaying sinusoid waveform in the frequency domain,˜ R ( f ; A R , f R , t , τ ) = A R e − πi ( f R − f ) t τ − + 2 πi ( f − f R ) . (16)The parameters of this signal are amplitude A R , fre-quency f R , starting time t , and decay time τ after whichthe envelope amplitude falls to e − of its original value.Using this waveform is equivalent to setting the initialphase of the wave to be 0, but this does not affect thegenerality of the results. In the simulations presentedhere, a central frequency f R = 150 Hz was chosen, sothat the peak of the ringdown would lie in the sensitiveregion of the frequency domain range. Ringdowns with τ = 1 s and 10 s were both added to the dataset, and thealgorithm was run as before to search for the presence ofa Newtonian inspiral. We stress again that the instru-mental ringdown was not part of the models consideredfor the computation of the Bayes factor.The results of this test are reported in Figure 4; thesignal-to-noise ratio of the instrumental ringdown shownon the x -axis is calculated in the same way as for the GWinspiral signal, using Equation 16, but replacing ˜ h with˜ R . The figure shows the response of the Bayes factor B SN of the inspiral model vs. noise model (for both of whichthe noise component is modelled as exactly Gaussian andstationary) to a dataset containing Gaussian noise and aninstrumental ringdown signal, which is much more sub-dued than that to an injected inspiral signal. The verticalaxis of the plot of Figure 4 is therefore the same quan-tity represented on the corresponding axis of Figure 3. Aringdown SNR of 1000 is necessary to produce a Bayesfactor comparable to that for an inspiral SNR of 4.5 (c.f.figure 3), when a ringdown with decay constant τ =10 sis used (for τ =1 s the effect is even smaller, and hardlynoticeable on the scale of the plot). E. Stationary and Gaussian noise with additionalPoissonian component
We now explore a different case in which the model’sassumption that the noise follows a Gaussian probabil-ity distribution is not fully accurate, by injecting an ad-ditional Poissonian noise component into the simulateddata. To accomplish this, a time-series stretch of noisewith amplitude drawn from a Poisson distribution in thetime domain, scaled by a factor 100 to increase its ampli-tude, and uniform random phase in the interval [0 , π ) isgenerated and Fourier transformed, before being addedto the standard Gaussian stationary noise generated inthe frequency domain. The root-mean-square of theGaussian and Poisson components are 4.6 and 0.3 re-spectively. The spectrum of the Poisson noise showed anapproximately uniform power across frequency in con-trast with the shaped power spectrum of the Gaussian −50050100150200 signal to noise ratio l og B a y e s f ac t o r FIG. 4: The Bayes factors B SN for inspiral model vs. noisemodel in the presence of an injected ringdown to simulate suchevents in instrumental noise. Shown here are the results fromtwo ringdowns, with decay times τ = 1 s (dashed line) and10 s (solid line), and varying signal to noise ratios plotted onthe x -axis. From these results it is clear that the algorithmis robust against interference from this source, as only the τ = 10 s caused an increase in the Bayes factor, and thenonly at very high SNR. −2 −1 Frequency (Hz) S t r a i n FIG. 5: The amplitude spectra of the Gaussian (upper black)and Poisson (lower grey) contributions to the data used insection III E. noise. The amplitude spectrum of the two noise contri-butions are shows in figure 5.The intention of this procedure is to simulate the ef-fects of outliers from the Gaussian noise distribution,which we know occur in the real interferometer data, andwhich we will be very unlikely to be able to model, mak-ing them in effect random. In this test we have chosen aPoisson distribution P ( λ ) with a mean λ of one point forevery ten time stamps, i.e. λ = 0 .
1. We have exploredboth the signal-free case and the situation in which aGW signal is added to the (Gaussian + Poissonian) noise,analogous of the studies presented in Secs. III B and III C, −2.0 −1.95 −1.9 −1.85 −1.8 −1.75024681012 log Bayes factor C oun t s FIG. 6: The results from injecting a Poissonian componentwith distribution P (0 .
1) into a segment of Gaussian noise (andno GW signal) and running the algorithm 100 times. Thepresence of Poisson noise does not increase the Bayes factorfor an inspiral signal, as the Poissonian noise does not matchthe template for this waveform. The results of plot should becompared with those obtained on pure Gaussian and station-ary noise reported in Figure 2. respectively.Figure 6 shows the results of analysing simply the com-bination of Gaussian and Poisson noise, without a GWsignal present. We can see that the change in Bayes fac-tor due to this additional noise component is minimal(cf. Figure 2), and does not trigger the signal modelat all, conversely it depresses the odds. The estimatedprobability of a signal being present remains low.Figure 7 shows the results of tests similar to those inSection III C (GW signal and noise). In comparison toFigure 3, the recovered Bayes factors are reduced for thesame signal-to-noise ratio – as an example, log B SN ≈
10 is reached for SNR ≈ . ≈ . l og B a y e s f ac t o r FIG. 7: The results from adding GW signals of varying signal-to-noise ratio onto the combined noise data from a Poissonianand Gaussian distribution, as described in the text. The re-sulting Bayes factors are lowered in comparison to Figure 3. rigourously quantified by performing Monte-Carlo simu-lations on real data.Since we do not expect the algorithm to perform better when assumptions of the models are violated, it is desir-able that it would degrade in the least harmful manner.The results of the tests presented here are therefore reas-suring: the Bayes factor will not spuriously generate falsealarms in the presence of this type of noise, but insteadwill simply not perform as well at detecting actual sig-nals. It is worth mentioning again that ideally one wouldalso compute the evidence for this noise model, and thatif this was done some of the sensitivity would be restored.
IV. CONCLUSIONS
Bayesian inference using evidence and odds ratios eval-uations provides a clear and justifiable means of deter-mining the probability of a signal being present in thedata: it is the optimal inference and, as we have shown(although only in simple cases), is also robust againstinterference from non-gravitational wave signals presentin the data, and against deviations in the noise profile.Through the method used in this paper the calculation ofodds ratios between hypotheses is made feasible withinuseful time-scales. The run-time of the code that we have developed can vary depending on factors such as the de-sired accuracy, whether or not a signal is present, theactual GW waveform and relevant number of parametersthat describe the model; however, to perform a single runof one of the analyses typical in this paper takes (much)less than one day on a single 2 GHz CPU. This is signif-icantly more efficient than other approaches consideredso far for Bayesian model selection, e.g. [22, 40]. Thisspeed makes it possible to perform thousands of odds ra-tios calculations per day on Beowulf-type clusters: thisgives the method a very good combination of sensitivityand speed which we hope will allow the method to beused as one of the techniques for follow-up studies in theanalysis of real data.It is our intention now to further develop this methodtowards the evaluation of odds ratios in real interferom-eter data, and integrate it with the existing analyses toprovide a robust Bayesian follow-up capable of determin-ing the odds of a signal being present. In order to achievethis, we need primarily to include additional models ofGW signals (such as post-Newtonian waveforms) andpossibly of instrumental artefacts, which will allow theanalysis to distinguish between different types of sourcesand to eliminate or detect contamination of the noisefrom unwanted sources.Another important feature of the method introducedin this paper which has not been discussed here but maybe useful in a combined Bayesian analysis is its abilityto find the maximum likelihood values of the parame-ters, which would integrate well with a Markov-ChainMonte-Carlo approach [24] to full parameter estimation,see e.g. [33, 34, 35, 36] and references therein. We intendto include such an interface with MCMC estimation toprovide a combined Bayesian follow-up package for inspi-ral analysis.
Acknowledgments
This work has benefited from many discussions withmembers of the LIGO Scientific Collaboration. We wouldlike in particular to thank Nelson Christensen and JamesClark for fruitful and stimulating discussions. Many ofthe simulations presented in this paper were performedon the Blue BEAR and Tsunami Beowulf clusters ofthe University of Birmingham. This work has beensupported by the UK Science and Technology FacilitiesCouncil. While at Northwestern University, AV was sup-ported by the David and Lucile Packard Foundation andby NASA grant NNG06GH87G. [1] B. C. Barish and R. Weiss, Physics Today , 44 (1999).B. Abbott et al. [The LIGO Scientific Collaboration],arXiv:0711.3041.[2] F. Acernese, et al., Class. Quant. Grav. , 385 (2004)[3] B. Willke, et al. Class. Quant. Grav. , 417 (2004) [4] R. Takahashi, et al. Class. Quant. Grav. , 403 (2004)[5] D. Sigg, Class. Quant. Grav. , S51 (2006); S. Waldman,Class. Quant. Grav. , S653 (2006).[6] R. Adhikari, P. Fritschel and S. Waldman, LIGO Tech-nical Document T060156-01 (2006). [7] [8] [9] C. Cutler and K. S. Thorne, in General Relativity andGravitation
Ed N. T. Bishop and D. M. Sunil (Singapore:World Scientific) p 72 (2002). arXiv:gr-qc/0204090.[10] B. Abbot et al. [The LIGO Scientific Collaboration],Phys. Rev. D , 082001 (2005).[11] B. Abbott et al. [The LIGO Scientific Collaboration],Phys. Rev. D , 062001 (2006).[12] B. Abbott et al. [The LIGO Scientific Collaboration],Phys. Rev. D , 082002 (2005)[13] B. Abbott et al. [The LIGO Scientific Collaboration],submitted to Phys. Rev. D , arXiv:0704.3368[14] B. Abbott et al. [The LIGO Scientific Collaboration],arXiv:0712.2050 .[15] R. Gouaty [for the LIGO Scientific Collaboration],”Detection confidence tests for burst and inspiralcandidate events”, submitted to Class. Quant. Grav.,LIGO Document LIGO-P080042-00-Z, available at [16] B. Allen, Phys. Rev. D , 062001 (2005).[17] A. Pai, S. Dhurandhar and S. Bose, Phys. Rev. D ,042004 (2001) [arXiv:gr-qc/0009078].[18] S. Chatterji, et al, Phys. Rev. D , 043003 (2007).[20] J. Clark, I. S. Heng, M. Pitkin and G. Woan, to appearin the proceedings of the , arXiv:0711.4039 .[21] A. C. Searle, P. J. Sutton, M. Tinto and G. Woan, to ap-pear in the proceedings of the , arXiv:0712.0196 .[22] R. Umst¨atter and M. Tinto, submitted to Phys. Rev. D, arXiv:0712.1030[23] K. C. Cannon, submitted to Phys. Rev. D[24] W. R. Gilks, S. Richardson, and D. J. Spiegelhalter, Markov chain Monte Carlo in practice , (Chapman &Hall / CRC, Boca Raton, 1996). A. Gelman, J B. Car- lin, H. Stern, and D. B. Rubin,
Bayesian data analysis ,(Chapman & Hall, 1997).[25] J. Skilling, in
AIP Conference Proceedings: 24th Inter-national Workshop on Bayesian Inference and MaximumEntropy Methods in Science and Engineering , Volume735 pp. 395-405 (2004).[26] [27] D. S. Sivia with J. Skilling,
Data Analysis: A BayesianTutorial (Second edition, Oxford University Press 2006)[28] P. Mukherjee, D. Parkinson, A. R. Liddle, Ap. J.
Probability theory: The logic of science (Cambridge University Press, 2003)[31] C. Cutler and E. E. Flanagan, Phys. Rev. D , 2658(1994)[32] V. Kalogera, K. Belczynski, C. Kim, R. O’Shaughnessyand B. Willems, Phys. Rept. , 75 (2007). C. Kim,V. Kalogera and D. R. Lorimer, arXiv:astro-ph/0608280[33] C. R¨over, R. Meyer, and N. Christensen, Class. Quant.Grav. , 4895 (2006).[34] C. R¨over, R. Meyer, and N. Christensen, Phys. Rev. D , 062004 (2007)[35] C. R¨over, R. Meyer, G. M. Guidi, A. Vicer´e and N. Chris-tensen, Class. Quant. Grav. , S607 (2007)[36] M. van der Sluys, et al. arXiv:0710.1897,arXiv:0805.1689.[37] R. Umst¨atter et al, Class.Quant.Grav. S901-S912(2005)[38] P. D. Welch, IEEE Trans. Audio and Electroacoustics,
AU-15 [40] N. J. Cornish and T. B. Littenberg, Phys. Rev. D76