{'fraction_non_alphanumeric': 0.02634801142805315, 'fraction_numerical': 0.008117545689537889, 'mean_word_length': 4.680893087161872, 'pattern_counts': {'":': 1, '<': 3, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
BayesWave: Bayesian Inference for Gravitational Wave Bursts and Instrument Glitches

A central challenge in Gravitational Wave Astronomy is identifying weak signals in the presence of non-stationary and non-Gaussian noise. The separation of gravitational wave signals from noise requires good models for both. When accurate signal models are available, such as for binary Neutron star systems, it is possible to make robust detection statements even when the noise is poorly understood. In contrast, searches for"un-modeled"transient signals are strongly impacted by the methods used to characterize the noise. Here we take a Bayesian approach and introduce a multi-component, variable dimension, parameterized noise model that explicitly accounts for non-stationarity and non-Gaussianity in data from interferometric gravitational wave detectors. Instrumental transients (glitches) and burst sources of gravitational waves are modeled using a Morlet-Gabor continuous wavelet frame. The number and placement of the wavelets is determined by a trans-dimensional Reversible Jump Markov Chain Monte Carlo algorithm. The Gaussian component of the noise and sharp line features in the noise spectrum are modeled using the BayesLine algorithm, which operates in concert with the wavelet model.


Introduction
The LIGO [1] and Virgo [2] gravitational wave detectors are currently undergoing major upgrades, with the goal of improving the broadband sensitivity by an order of magnitude [3,4]. But it takes more than lasers, mirrors, photodiodes and suspension systems to detect gravitational waves. Another key component is the signal processing used to tease the faint gravitational wave signals out of the instrument noise.
Searches for gravitational waves using data from the initial LIGO-Virgo science runs have shown that non-Gaussian features in the data, such as noise transients or "glitches", impact the ability to detect weak signals as the various detection statistics used in the searches end up with heavier tails than they would for Gaussian noise [5,6,7]. The standard procedure for accounting for these effects is to perform Monte Carlo trials with scrambled data produced by introducing time lags that break signal coherence between the detectors in the network. These trials are used to produce background probability distributions for search statistics -such as the cross correlation of signal templates with the data -which can then be used to assess the significance of events detected in the zero-lag data. Extensive testing with simulated signal injections into the data has shown that this standard approach is reliable and robust, though with a sensitivity that is less than what would be possible if the noise were stationary and Gaussian. In addition to degrading detection efficiency, the non-Gaussianity of the noise can also impact parameter estimation of quantities like the sky location and masses of a binary system [8].
The detection statistics used in the searches are often motivated by analytical studies assuming Gaussian noise, and then the effects of non-Gaussianity are taken into account by using the Monte Carlo derived distributions estimated from the actual detector data. In effect the data are used to infer the noise properties, or at least, the properties of the noise projected onto the signal manifold defined by the detection statistic. Here the noise modeling is taken a step further in a "holistic" approach to gravitational wave detection that employs a sophisticated noise model in addition to a signal model. The data are used to jointly infer the parameters of both the noise and signal models, following the motto "model everything and let the data sort it out." In an earlier study we introduced a discrete orthogonal wavelet model to account for short-duration, non-stationary features in the instrument noise [9]. This approach was developed further in collaboration with P. Baker, and a description of the study can be found in Baker's Ph.D. thesis [10]. From these foundations we have developed the BayesWave algorithm for the joint detection and characterization of instrument glitches and gravitational wave bursts. The name "BayesWave" is derived from a concatenation of "Bayesian" and "Wavelets." The key to the success of the new algorithm is the adoption of a non-orthogonal, continuous Morlet-Gabor wavelet basis, which provides an extremely compact representation of bursts and glitches, and has the additional advantage that the entire analysis can be implemented in the Fourier domain. Bayesian model selection and parameter inference is implemented using a trans-dimensional Reversible Jump Markov Chain Monte Carlo (RJMCMC) [11] algorithm that varies the number and parameters of the wavelets used to model the glitches and gravitational wave signals. Using a variable number of wavelets builds in a natural parsimony that balances the complexity of the model against the quality of the fit. Much of the BayesWave algorithm has been re-purposed from techniques we developed to detect and characterize thousands of overlapping signals from compact white dwarf binaries in simulated data from a future space based gravitational wave detector [12,13]. The BayesWave algorithm is run in concert with the BayesLine spectral whitening algorithm, which we describe in a companion paper. The BayesLine algorithm models the noise spectrum using a combination of a smooth spline curve and a collection of Lorentzians to fit sharp line features. BayesLine employs a RJMCMC algorithm to vary the number and placement of the spline control points and the number of Lorentzians and their parameters.
We begin by reviewing the application of Bayesian inference in gravitational wave astronomy, with an emphasis on the importance of choosing an appropriate likelihood function. Next we describe the glitch and gravitational wave burst models and lay out our reasons for choosing the Morlet-Gabor wavelet basis. This is followed by an investigation of various time-frequency clustering priors using simulated Gaussian noise and simulated gravitational wave signals. The effectiveness of the BayesWave algorithm is then demonstrated using real-world data from the sixth LIGO science run, including glitch detection and reconstruction, and the recovery of hardware and software signal injections.

Bayesian Inference for Gravitational Wave Astronomy
The posterior distribution function p(h|s, M ) summarizes the properties of a gravitational waveform h that can be inferred from data s under model M : where p(h|M ) is the prior distribution for h in model M , p(s|h, M ) is the likelihood of the data for the given waveform, and p(s|M ) = p(h|M )p(s|h, M )dh is the marginal likelihood, or evidence, for the model M under consideration. Here s is a vector of data collected from a network of gravitational wave detectors, which is a combination of instrument noise n and a gravitational wave signal h such that s = R h + n. Here R is a time-delay operator that describes the network response to the two gravitational wave polarizations h.
The Bayesian approach to data analysis is entirely mechanical. Once you have specified your signal model and the likelihood function, the posterior distribution can be computed using numerical techniques such as Markov Chain Monte Carlo [14] or Nested Sampling [15]. From the posterior distribution function we can extract point estimates and "error bars" for the parameters that characterize the waveform h, and by comparing the evidence for alternative models, we can quote odds ratios for possibilities such as "the data contains the signal from a black hole merger" and "the data is purely instrumental noise." On the other hand, the quality of what you get out is only as good as what you put in. If your signal model is flawed, or your likelihood function is too simple, the posterior distribution will produce misleading results. The frequentist approach is more forgiving, as the likelihood function is only used to motivate the form of a search statistic, which is then tuned using signal injections and time slides of the data.
The likelihood function is defined by demanding that the residual formed from subtracting the gravitational wave model from the data, r = s − R h, is consistent with the model for the noise n. For stationary, Gaussian noise this leads to the standard likelihood function [16] where C is the noise correlation matrix and N is the length of the data s. For stationary noise, transforming the data to the Fourier domain yields a diagonal correlation matrix, C f f = S n (f ) δ f f , which allows for faster likelihood evaluations. The existing Bayesian parameter estimate algorithms that have been applied to LIGO/Virgo data [8] assume that the noise is stationary and Gaussian and use off-source estimates for the noise spectrum S n (f ). However, non-stationary, non-Gaussian features such as glitches or wandering lines can bias this type of analysis [9,17]. The need for more realistic noise modeling has been emphasized in several studies. Allen et al. [18,19] have argued that parameterized noise models are needed to derive robust search statistics. They found that distributions with heavier tails than a Gaussian were more robust against noise transients, a result that has been confirmed by Röver using a Student-t likelihood function [20]. Clark et al. [21] included a model for Sine-Gaussian instrument glitches as an alternative hypothesis when computing Bayesian odds ratios for gravitational wave signals from pulsar glitches. Clark et al. also suggested that it would be valuable to have a classification of instrument glitches that could be used to construct better models for the instrument noise. Similar "glitch hypotheses" were considered by Veitch and Vecchio [22] in a study of Bayesian model selection applied to the search for black hole inspiral signals. The possibility of subtracting instrument glitches from the data using a physical model of the detector has been investigated by Ajith et al. [23,24]. Principe and Pinto [25] have introduced a physically motivated model for the glitch contribution to the instrument noise, and have used this model to derive what they refer to as a "locally optimum network detection statistic" [26]. In their approach the glitches are treated in a statistical sense, while our approach directly models the glitches present in the data. The possibility of directly deriving likelihood functions from the data in a Bayesian setting has previously been considered by Cannon [27]. In Cannon's approach the data is first reduced to an n-tuple of quantities such as the parameters produced by a matched filter search. Then, using time slides and signal injections, the likelihood distributions for the signal and noise hypotheses are directly estimated from the data. These are then used to estimate the posterior probability that a measured set of parameters corresponds to a gravitational wave event.
The BayesWave algorithm models the data s as linear combination of gravitational wave signals, h; stationary, colored Gaussian noise, n G ; and glitches, g: The spectrum of the Gaussian component n G is modeled by the BayesLine algorithm, running on relatively short stretches of data (tens to hundreds of seconds) to allow for broad spectrum, long term non-stationary behavior. We have experimented with using more complicated noise models for n, such as the sum of Gaussians proposed by Allen et al. [18,19], but found little benefit from the added complexity [9]. Studies of the LIGO/Virgo data show that glitches are structured features that occupy a compact region in the time-frequency plane. It is natural therefore to use some kind of wavelet representation for the glitches. This could be a continuous wavelet representation such as the Morlet-Gabor "Sine-Gaussian" wavelets used in the LIGO Omega pipeline [28] and for performing time-frequency decomposition of the data via Q-scans [29], or a discrete wavelet representation such as those used in the Kleine-welle studies of environmental channels or the coherent WaveBurst burst search pipeline [30,31,32]. Principe and Pinto [25] have shown that typical glitches can be represented as the sum of a small number of Morlet-Gabor "glitch atoms". The BayesWave algorithm uses a transdimensional RJMCMC to model glitches and signals with a variable number of Morlet-Gabor wavelets, where the number and parameters of the wavelets are not pre-specified, but rather, determined from the data. In principle the RJMCMC algorithm can produce posterior distributions for the components of the nested model (3) which can be used for model selection. For example, we can define a Gaussian noise model to be the state that employs zero glitch or signal wavelets. Similarly, we can define a pure glitch model to be one that uses zero signal wavelets, and likewise, a pure signal model to be any state with zero glitch wavelets. In practice this is only useful for cases where the relative evidence for each model (the pairwise Bayes Factors) are within a factor of a thousand or so. Beyond that point, the less favored models are visited so infrequently that the RJMCMC derived Bayes Factors become suspect. It is possible to extend the RJMCMC analysis to higher Bayes Factors by imposing a false prior that weights the less favored model more heavily, then correct for the false prior when computing the Bayes Factor [33], but it gets increasingly difficult to achieve good mixing. To handle situations with large Bayes Factors we use thermodynamic integration to directly calculate the evidence for each model. The various models -Gaussian noise, glitches plus Gaussian noise and GW signal plus Gaussian noise -are themselves composite models formed by the union of all possible spectral estimates from BayesLine, and all possible combinations of wavelets in the glitch or signal models. Each composite model is explored by a RJMCMC, and thermodynamic integration is used to compute the evidence for the model. We are unaware of any other applications where RJMCMC and thermodynamic integration have been combined in such a fashion. The evidence is used to select the model that best describes the data.

Glitch and Signal Modeling
Earlier versions of the BayesWave algorithm [9,10] used orthogonal bases of discrete Meyer wavelets. Initially we used a dyadic basis uniformly spaced in log frequency [9], and later we experimented with a binary basis uniformly spaced in the frequency [10]. The parameters in these models were the number, amplitude and time-frequency location of the "active" wavelets. There are several advantages to using a discrete wavelet grid: the parameters of the wavelet model are completely uncorrelated, resulting in good mixing of the Markov Chains; and it is easy to identify clusters of wavelets. There are also several disadvantages: the shape of the pixels in time-frequency are fixed in advance by the decomposition level of the wavelet transform, leading to suboptimal representation of signals/glitches [10]; and projecting the signal wavelets onto each detector in the network required the application of a computationally expensive time translation operator [10]. We eventually concluded that the negatives outweighed the positives and abandoned discrete wavelets in favor of non-orthogonal Morlet-Gabor wavelets which have a number of nice properties: their shape in time-frequency is variable and can adapt to fit the signal; they have an analytic Fourier representation so the entire analysis can be performed in the Fourier domain; and the time shifts needed to project the signal wavelets onto the various detectors are described by a simple phase shift. The lack of orthogonality is not a major impediment to achieving good mixing in the Markov Chains, and clusters can be defined using a metric distance on the wavelet space.
Morlet-Gabor wavelets can be expressed in the time-domain as with τ = Q/(2πf 0 ). Here A is the amplitude, Q is the quality factor, t 0 and f 0 are the central time and frequency, and φ 0 is the phase offset. Their Fourier transform is given by Since we typically restrict Q > 1, the wavelets are highly peaked around f = f 0 , and we may neglect the second term in the Fourier transform. Morlet-Gabor wavelets have the minimum time-frequency area, πσ t σ f = 1/2, allowed by the Heisenberg-Gabor uncertainty principle (their form is closely related to coherent states in quantum mechanics). The glitch model employs a variable number of independent wavelets, N g,i in each detector. The signal model can either use two sets of independent wavelets to model h + and h × at the Geocenter, or for signals with a definite polarization, a single set of wavelets to model h + , and an ellipticity parameter . In the frequency domain we define h × = h + e iπ/2 , with the polarization going from linear to circular as runs from 0 to 1. Our motivation for assuming elliptical polarization is two-fold: Firstly, the the small number of detectors available in the early advanced detector era, and the near alignment of the two LIGO detectors, results in poor polarization sensitivity, making it difficult to reconstruct h + and h × separately. Secondly, we expect most sources of gravitational waves to show some degree of polarization. For example, the signals from binary inspirals are elliptically polarized with ellipticities that are approximately constant across multiple wave cycles, while signals from core collapse supernovae are thought to be dominantly linearly polarized. In what follows we will assume that the signals are polarized, and that is time independent. The Geocenter signal wavelets are projected onto each detector in the network using sky-location-dependent time-delays ∆t(θ, φ), and antenna patterns F + (θ, φ, ψ), F × (θ, φ, ψ): where All N s wavelets in the signal model share the same set of extrinsic parameters describing the sky location (θ, φ), ellipticity , and polarization angle ψ. A gravitational wave signal can be fit equally well (in terms of the likelihood) by the glitch or signal models, but the evidence for the signal model will be larger as it provides a more parsimonious description of the data. For a signal described by N s wavelets incident on a N d detector network, the signal model uses D s = 5N s + 4 parameters (the 5 extrinsic parameters per wavelet and the 4 common extrinsic parameters), while the glitch model uses D g = 5N d N s parameters. When N s 1 or N d 1, the glitch model requires vastly more parameters than the signal model. When N s = 1 and N d = 1 the extrinsic parameters are irrelevant and D s = D g , which makes it impossible to tell the difference between an un-modeled gravitational wave signal and a glitch using a single detector. The most marginal case for detection is a two detector network and a single sine-Gaussian signal, as then D s = 9 and D g = 10. However, despite the near equality of the model dimensions in this case, the signal model will still be preferred as not all parameters are created equal -models pay a greater penalty for parameters that are well constrained (measured by the ratio of prior to posterior volume). Intrinsic parameters such as t 0 and f 0 fall into this category.

Priors
Having defined the likelihood and the building blocks of the signal and glitch models, all that remains is to specify prior distributions on the model parameters. For the extrinsic parameters of the signal model we adopt the prior that sources are uniformly distributed on the sky, with ellipticities uniform in the range ∈ [0, 1], and polarization angles uniform in the range ψ ∈ [0, π]. Similarly, we take the prior on the intrinsic parameters t 0 and f 0 to be uniformly distributed across the observation time and frequency band being analyzed. These priors should be modified if we have an electromagnetic counterpart such as a core collapse supernovae, which would pin down the sky location and perhaps the ellipticity, and narrow the range on t 0 . We take the quality factors to be uniform in the range Q ∈ [Q min , Q max ], where we typically set Q min = 2 and Q max = 100. The phase parameters are taken to be uniform in the range φ 0 ∈ [0, 2π]. We take the number of parameters in the glitch and signal models to be uniform across a wide range, typically N s ∈ [0, 100] and N g ∈ [0, 100 × N d ]. When running on the signal model the ranges are set such that N s ∈ [1, 100], N g = 0, while the glitch model has N s = 0 and N g ∈ [1, 100 × N d ] (in other words, a single glitch wavelet in any detector is acceptable). The Gaussian noise model has N g = 0 and N s = 0. All that remains is to specify an amplitude prior for the wavelets. In this instance, a uniform prior has little to recommend it, and we need to specify physically well motivated priors for the signal and glitch amplitudes. In each case we specify a prior on the amplitude of individual wavelets, but these may be supplemented by priors on the total power in the signal and glitch models.

Glitch amplitude prior
From previous studies we know that loud glitches are less common than quiet glitches, and that glitches with amplitudes that are comparable to the Gaussian noise level tend to blend into the background. These considerations motivate an amplitude prior expressed in terms of the signal-to-noise ratio (SNR) of the wavelet. For glitches we adopt the form which peaks at SNR = SNR * . The SNR of a Morlet-Gabor wavelet can be estimated as which means that the SNR prior is actually a joint prior on A, f 0 , Q. A non-flat prior on the SNR will result in a bias in the recovered amplitude which is tolerable so long as it is not too large. To estimate the bias, consider a single sine-Gaussian injection with signal-to-noise SNR, and further assume that the template exactly matches the signal in all parameters other than the amplitude. The posterior distribution for SNR then has the form Note that for large signal-to-noise, the likelihood overwhelms the prior, and in the limit SNR 1 the posterior distribution peaks at  Figure 1. The upper panel shows the SNR prior for the glitch wavelets (in red) and the signal wavelets (in blue) for the choice SNR * = 4. The lower panel shows the fractional bias in the inferred SNR caused by using these non-uniform SNR priors. The shaded regions fold in the 1-sigma statistical error in the SNR estimates, which are always far greater than the systematic biases from the non-uniform priors.

Signal amplitude prior
For gravitational wave signals we expect the sources to be distributed roughly uniformly in volume, which implies a prior on the distance D that scales as p(D) ∼ D 2 . Since the SNR of a signal scales as SNR ∼ 1/D, the corresponding prior on SNR scales as p(SNR) ∼ SNR −4 . One possibility is to use this as the prior on the individual wavelet SNRs, but such a prior is improper (can not be normalized), and it favors adding large numbers of undetectable wavelets to the signal model. While we are not big fans of so-called "Malmquist priors" that seek to account for observational selection effects (to us this smacks of double counting what we already learn from the likelihood), there is no great harm done by cutting off the distribution at some minimum SNR as the natural balance between model fidelity and model complexity will mean that only wavelets with SNR > 3 will play any significant role in the signal model. Rather than adopt a hard SNR cut-off, we instead choose a function that smoothly goes to zero as SNR → 0, and approaches the SNR −4 scaling for high SNR. One such choice is This distribution peaks at SNR = SNR * . The SNR of a signal wavelet can be estimated using (9), but with the power spectral density of a single detector S n (f 0 ) replaced by the network averagē Thus our SNR prior for the signal model is actually a joint prior on A, f 0 , Q, , ψ, θ, φ. As with the glitch model, there will be a bias away from the maximum likelihood solution for the inferred SNR caused by the non-uniform prior on SNR. Figure 1 compares the glitch and signal SNR priors for the choice SNR * = 4, and shows the fractional bias in the inferred SNR caused by using a non-uniform prior. In both cases, the bias is tiny compared to the statistical uncertainty.

Cluster Priors
In addition to the individual wavelet priors described above, it is possible to introduce non-local priors that apply to collections of wavelets. On physical grounds we expect GW signals and glitches to produce coherent clusters of power in a time-frequency representation. This suggests imposing priors that favor concentrated clusters of wavelets over isolated wavelets. One approach that we considered uses a distance measure to decide which wavelets are nearby in the time-frequency plane, then groups neighboring wavelets into a cluster. Treating the wavelets as Gaussian packets in time-frequency space with density where ∆f = f 0 − f , ∆t = t 0 − t, σ t = τ and σ f = 1/πτ , and computing the overlap between two wavelets yields the distance measure We then set a distance threshold such that all wavelets within ds < δ of each other were considered to be in a cluster, and used the Flood Fill algorithm to identify clusters. From this we were able to extract the number of clusters, the number of wavelets in each cluster, and the signal-to-noise ratio of each cluster. We experimented with a number of priors that used these quantities to favor groupings of wavelets, but each attempt ran into various difficulties. In general there was a strong tendency to form large clusters, even in Gaussian noise, and these clusters would typically have wavelets with extreme Q's producing an unphysical mix of both very tall and very long pixels that linked the clusters together. The tendency to over-cluster could be balanced to a certain extent by adding a prior that disfavored configurations with a large number of wavelets, but the bias towards extreme wavelet shapes remained. The approach we settled on employs a "proximity prior" based on a sum of probability densities surrounding each active wavelet. Since each wavelet covers a region of the time-frequency plane defined by the density (14), we expect neighboring wavelets to be close, but not too close. This motivates a probability density described by a "hollowed-out" Gaussian of the form where α > β sets the overall width of the distribution and β sets the size of the hollowedout region. We expect that choosing α in the range 2-4 and β in the range 0.5-1 will result in an effective proximity prior for most signal and glitch morphologies, but we have yet to conduct an extensive study of how these choices affect the performance. The full proximity prior is given by the normalized sum of densities of the form (16) surrounding each active wavelet, to which we also add a uniform density component to allow exploration of the entire time-frequency volume. We set the uniform fraction at γ, and the hollowed-out Gaussian fraction is 1−γ. Rather than choosing a fixed fraction γ, we instead choose a fixed density contrast between the two components, which results in γ taking different values depending on the time-frequency volume being analyzed. The proximity prior favors states with large numbers of nearby wavelets. This can be seen most clearly in prior-only studies where the likelihood is set equal to a constant. We find that the proximity prior leads to an approximately power law distribution in the number of active wavelets, as opposed to the uniform distribution in the number of active wavelets found when using a uniform time-frequency prior. The precise form of the wavelet number distribution depends on the choice of α, β, γ and the size of the time-frequency volume being analyzed. Figure 2 shows the distribution of the number of active wavelets in prior-only studies with (α, β) = (4, 1) and (α, β) = (2, 0.5). In each case the range in the number of active wavelets was set at N = [0, 19], and the analysis covered the frequency range f = [10, 512] Hz and a time interval of 8 seconds, resulting in a time-frequency volume of TFV = 4016. The density contrast in the proximity prior was set such that γ = 1/(4016/TFV + 1) = 0.5. The more aggressive proximity prior with (α, β) = (2, 0.5) leads to a steeper trend in the number of active wavelets being used. Since our preferred physical prior on the number of active wavelets is uniform, we include a prior of the number of wavelets that compensates for the over-clustering caused by the proximity prior. The obvious choice is to use the inverse of the distribution shown in Figure 2, which results in a uniform distribution in the number of active wavelets when the likelihood is set equal to a constant. Figure 3 shows the log of the proximity prior density for a cluster of wavelets randomly chosen from a MCMC analysis. The upper left panel shows what the prior would look like when proposing to add an additional wavelet to the cluster. The other panels show what the prior looks like for each individual wavelet in the cluster i.e. what the prior looks like when proposing to add the third wavelet (shown in white) to the cluster. Note that the dynamic range in the proximity prior between the peak value and the minimum is of order ∆ log p ≈ 6, which means that the improvement in the log likelihood needed to add a new wavelet near the cluster is significantly lower than elsewhere in the time-frequency plane. Figure 4 shows randomly chosen snapshots of the the wavelet distributions for states with 30 wavelets from runs where the log likelihood was set equal to a constant. The wavelets are grouped into clusters using the Flood Fill algorithm with a distance tolerance of ds < 4 and assigned a common color. Note that since there are only a finite number of colors in the color table being used to make the plot, some wavelets get plotted in the same color even though they are not in the same cluster. The upper panel is with a uniform time-frequency prior, while the lower panel is with the proximity  prior. As expected, the proximity prior shows a distinct preference for forming clusters composed of many wavelets. Figure 5 shows quantiles of the match between an injected black hole merger signal h and the wavelet reconstruction h found by BayesWave for a flat time-frequency prior and a proximity prior with α = 4, β = 1, γ = 0.5. The match, or overlap, is defined: and can range between -1 and 1. The black hole merger signal included inspiral, merger and ringdown for a system with spins aligned with the orbital angular momentum, and spin magnitudes χ 1 = 0. we plan to perform a more extensive study of how the proximity prior performs with different choices of α, β, γ for a variety of simulated signals. For the remainder of the paper we will adopt the fiducial values α = 4, β = 1, γ = 0.5. Figure 6 shows quantiles of the match as a function of signal-to-noise ratio for a simulated Gaussian enveloped white noise burst with central frequency f = 180 Hz and spreads σ f = 20 Hz, σ t = 0.1 s. Once again the proximity prior significantly outperforms the uniform time frequency prior in terms of the accuracy of the waveform recovery.

RJMCMC Implementation
The BayesWave algorithm uses a trans-dimensional RJMCMC technique [11] to explore the models and their many parameters. Transitions to new states of the chain are drawn from a proposal distribution, and to achieve efficient sampling the proposals need to be well adapted to the problem. The most challenging transitions are those that involve adding or removing a wavelet (the trans-dimensional steps), and birth/death moves that simultaneously remove one wavelet and add another. To ensure a high acceptance rate, BayesWave uses a mixture of custom tailored proposal distributions and other tricks:

Maximization during Burn-in
The RJMCMC algorithm does a very good job of exploring the region around the peaks of the posterior distribution, but it fares less well at finding the peaks. To overcome this limitation, standard FFT based tricks are used to maximize over t 0 , φ 0 and A for each Match SNR Figure 6. The median match for the wavelet model as a function of signal-to-noise ratio for a simulated white noise burst. The red line shows the median match and the symmetric 90% credible interval for a uniform time-frequency prior, while the blue line and shaded region shows the same quantities for the proximity prior with α = 4, β = 1, γ = 0.5 wavelet during the "burn-in" stage. Samples from this part of the chain are discarded since they do not obey the detailed balance condition required for Markov chains. The maximization allows signals and glitches to be identified very quickly.

Parallel Tempering
Convergence to the target distribution is improved by using the Parallel Tempering approach [35], where multiple chains are run in parallel at different "temperatures" T , where the likelihood function is modified such that p(s|h) → p(s|h) 1/T . Hot chains explore the entire prior volume, while cooler chains map the peaks of the posterior distribution. Exchanges between chains at different temperatures help ensure that the entire posterior distribution is thoroughly explored.

Fisher matrix proposals
A Fisher information matrix Γ can be computed for each wavelet, and proposals for updating the parameters of that wavelet can be drawn from the multivariate Gaussian distribution where the ∆x i = x i − y i denote displacements in the five intrinsic parameters that describe each wavelet. To speed calculation we use an analytic approximation to the Fisher matrix that can be derived along the same lines as the SNR in (9). Dropping terms down by factors of e −Q 2 relative to leading order, the Fisher matrix for a single wavelet using the parameters {t 0 , f 0 , Q, ln A, φ 0 } is given by To a good approximation, two of the eigenvectors of this matrix lie along the parameter directions t 0 (with eigenvalue (2πf 0 SNR) 2 ) and φ 0 (with eigenvalue (SNR/Q) 2 ). The other eigenvectors are complicated combinations of unit vectors in the {f 0 , Q, ln A} parameter directions.

Extrinsic parameter proposals
Changes to the extrinsic parameters which describe the location and orientation of the source are proposed using a cocktail of three distributions. The most heavily employed proposal uses the Fisher information matrix of the Geocenter waveform similar to what is used for the intrinsic parameters. There is no simple analytic expression for derivatives of the signal model with respect to extrinsic parameters, so the elements of the Fisher matrix are computed by finite differencing. Numerical differentiation is computationally costly so the Fisher matrix is treated as constant over several chain iterations before being updated. The Fisher matrix proposal assumes a multivariate Gaussian likelihood which is a poor approximation for parameters that encode the position and orientation of the signal. Ground-based gravitational wave detector networks acquire most of the positional information by differences in the arrival time of the signal at each observatory. For two detector networks (the minimum required for confident detection) timing-only considerations produce a degenerate ring on the sky of constant time delay. This degeneracy is broken by the different phase and amplitude measured in each detector but the resulting distributions are not remotely Gaussian (for example, see Figure 15).
We incorporate this known feature of gravitational wave posteriors into a proposal which rotates the current sky-location uniformly along the ring of constant time delay. Just changing the sky location is not enough because of the degeneracy-breaking correlations between position and orientation parameters. The response in each detector has the form Under a change of sky location we additionally demand that the change in the polarization phase φ p = arctan( F × /F + ) is the same for all the detectors in the array. This is always possible for a two detector network, but only sometimes for networks with three or more detectors. In general, we demand that the change in polarization phase is as similar as possible in all the detectors. This overall phase change can then be absorbed by shifting the phase φ 0 of each wavelet, which ensures that the signal in each detector for the new set of extrinsic parameters always closely matches that for the current set of extrinsic parameters, leading to a high acceptance rate for extrinsic parameter proposals. Rounding out the extrinsic parameter proposals is a uniform draw from the prior on all parameters. Proposals which make such large changes to the extrinsic parameters are seldom, if ever, accepted by the cold chains but are necessary for the hotter chains in the parallel tempering scheme to ensure that they fully sample the prior. The full prior exploration is important to help the cold chains fend off getting trapped in local features of the posterior and to get an accurate estimate of the expectation value of the likelihood as a function of temperature that is used in the thermodynamic integration.

Trans-dimensional and birth-death moves
When proposing to add a new wavelet, or swapping out an existing wavelet, the choice of time-frequency location for the proposed wavelet is key to getting the move accepted. BayesWave employs two kinds of time-frequency proposal to help facilitate these moves. The first is based on a time-frequency map which favors adding wavelets in regions with excess power. The second is a proximity proposal, which preferentially seeks to add new wavelets near to, but not on top of, existing wavelets.
The time-frequency map can be constructed in a variety of ways. One way is to use a normalized discrete scalogram of the whitened data in each detector. The glitch model in each detector uses the whitened scalogram for that detector, while the signal model uses the average of the scalograms for each detector in the network. Another way of constructing a time-frequency map is to perform a pilot run for each model (signal and glitch) and form a histogram of where the wavelets get placed in time-frequency space. Both methods for constructing time-frequency maps have proven to be successful in testing. The current implementation of the code uses the histogram based approach.
The proximity proposal distribution takes the same functional form as the proximity prior, though possibly with different values for the parameters (α, β, γ). The idea behind this proposal is that new wavelets are more likely to be needed in the vicinity of an existing cluster.
Once the central time t 0 and frequency f 0 have been drawn from either the timefrequency map or proximity proposal, values for the other wavelet parameters need to be proposed. The quality factor Q and phase φ 0 are drawn from their uniform prior range, while the amplitude A of the new wavelet is proposed by drawing a value for the proposed SNR from the amplitude prior (8), and solving for A by inverting the expression in (9). For wavelets in the signal model, the extrinsic parameters (sky location, polarization angle and ellipticity) are inherited from the existing signal model wavelets, or if the signal model currently has no wavelets, these are drawn from their prior distributions. Care must be taken to keep track of the proposal density factors in the Metropolis-Hastings ratio. The implementation is tested by setting the likelihood equal to a constant and seeing if the algorithm recovers the prior distributions. It is hard to over-stress how useful and important this test is when developing a code as complex as BayesWave.

Bayesian Evidence
In order to compare various models for the data we need to be able to compute the posterior probability for each model. The (un-normalized) posterior probability for a particular model M is given by the product of its prior odds p(M ) and its marginal likelihood or evidence p(s|M ). The marginal likelihood for a model is given by an integral of the model likelihood weighted by the prior distribution of the model parameters: This integral is notoriously difficult to evaluate when the model has a large number of parameters. Several ingenious techniques have been developed to estimate the evidence, including nested sampling and thermodynamic integration. An alternative approach is to avoid a direct calculation of the evidence and instead use a trans-dimensional RJMCMC algorithm to simultaneously explore both the space of models and the parameters of each model. The probability for each model can then be estimated from the relative frequency of visits to each model. BayesWave uses a RJMCMC routine to explore the full range of models (GW signals, glitches and Gaussian-noise), so it would seem that we get the model probabilities for free. We can define four disjoint composite models: (1) Gaussiannoise only; (2) Gaussian-noise and glitches -states with one or more glitch wavelets in use across the network; (3) Gaussian-noise and a GW signal -with one or more signal wavelets in use; (4) Gaussian-noise, glitches and a GW signal -with one or more signal wavelets and one or more glitch wavelets in use; and use the fraction of iterations spent in each model as an estimate of the posterior probability for each model. The challenge is to get the RJMCMC algorithm to efficiently explore each composite model and to transition between models. It is worth emphasizing that the composite models are of variable dimensionality -the Gaussian noise model allows for a variable number of spline control points and Lorentzian line features, and the glitch and GW models allow for a variable number of wavelets. Thus the RJMCMC algorithm has to marginalize over the internal model degrees of freedom in addition to exploring the composite model space. In situations where there are near-coincident loud glitches in multiple detectors or a loud gravitational wave signal, it becomes very difficult to come up with proposals that allow for efficient transitions between the all-glitch and all-signal models. One way around this is to do pairwise comparisons between a sub-set of models, such as Gaussian-noise versus Gaussian-noise and glitches or Gaussian-noise versus Gaussian-noise and gravitational waves. The ratio of the Bayes factors from these cases can the be used to compute the glitch/gravitational wave Bayes factor. However, even the pairwise model comparisons become challenging for loud events as the chains spend very little time in the Gaussian-noise model, and make few transitions between the models, which results in large statistical errors in the Bayes factors. To overcome this problem we adopted a hybrid approach that uses a RJMCMC algorithm to marginalize over the parameters in each composite model (including the number of spline points, wavelets etc), and Thermodynamic Integration [34] to compute the evidence for the composite model via an integral of the average log likelihood as a function of "inverse temperature" β = 1/T : Here the expectation value for the log likelihood for a chain with inverse temperature β, E β [ln p(s| λ, M )], averages over both the model parameters and the model dimension of the composite model M . To compute the integral (22) we first change variables to log β and discretize using a uniform spacing in log β. Since we already use parallel tempering to promote mixing of the Markov chains, we have a ready made temperature ladder with which to compute the evidence. Figure 7 compares Bayes factors computed using thermodynamic integration and from RJMCMC transitions for a simple test case with a simulated Sine-Gaussian "glitch" in simulated Gaussian noise. The two methods agree very well for low SNR glitches, but beyond SNR = 7.5 our implementation of the RJMCMC technique fails as the Markov chains rarely transition away from the glitch model. It is possible to push the RJMCMC technique to higher SNR and higher Bayes factors by introducing pseudo priors on the models so that the chains spend roughly equal time in each model [33], but even then it is difficult to find proposals that allow for frequent transitions between the models.

Evidence Error Estimation
When using Bayesian evidence to perform model selection, such as deciding if a certain non-Gaussian feature in the data is an instrumental artifact or a gravitational wave signal, it is important to have reliable estimates of the evidence and the error in the evidence estimate. One standard technique is to repeat the analysis multiple times using different random number seeds and compute the variance of the evidence estimates, but this can be extremely costly, especially if we demand better than ∼ 30% accuracy in the error estimates. A much cheaper alternative is to develop internal error estimates that can be computed from a single run. The Nested Sampling technique [15] comes with an internal error estimate, but we were unable to find similar estimates for the RJMCMC and thermodynamic integration techniques in the literature, so we derive them here.
where p ij is the transition probability between states i → j. The left eigenvalues of the transition matrix are λ 1 = 1 and λ 2 = 1 − (p 12 + p 21 ) < 1. As we iterate the Markov chain the components of the initial state that lie along the eigenvector corresponding to the smaller eigenalue λ 2 decay away, and we are left with the stationary state which is given by the λ 1 eigenvector with components The Bayes factor B 12 is simply the relative probability for the two models in the stationary state: B 12 = p 21 /p 12 . In practice we can only estimate the transition probabilities from number of counts in each model N ii and the number of transitions between models N ij . The joint likelihood of observing the collection of transitions d = (N 11 , N 12 , N 21 , N 22 ) is given by with N i = N i1 + N i2 . Setting ∂ p ij p(d|p ij ) = 0 yields the maximum likelihood estimatê The maximum likelihood estimate for the Bayes factor is given byB 12 =p 21 /p 12 = N 1 /N 2 since in the large N limit N 12 = N 21 . The covariance in these estimates can be found from the inverse, C µν , of the Fisher information matrix Γ µν = −∂ µ ∂ ν ln p(d|µ).
Here we have used the shorthand µ = p ij . The covariance matrix has the form and using the identity it follows that We have verified the accuracy of this error estimate using extensive Monte Carlo studies of two state Markov chains. We also find very good agreement between this error estimate and the brute-force error estimate found from performing multiple runs (see Figure 8). The one caveat is that the marginalization over the model parameters for each model needs to be efficient, otherwise our starting assumption that the transitions can be modeled by a two state Markov chain is violated. If the within-model mixing is poor, the estimate (29) provides a lower bound on the error. For poorly mixed chains the error estimate may need to be inflated by as much as a factor of 2.

Thermodynamic Integration Evidence Error Estimate
The log evidence for model M given data s can be computed: In practice the integral is evaluated by discretizing ln β and estimating the expectation value from the average log likelihood for chains running at the various inverse temperatures. Numerical error will result from both the discretization of the integral and from the statistical error in the estimate of the average log likelihoods. We have developed a numerical technique for estimating both sources of error that is described below. It is also useful to derive an analytic estimate for the statistical error that can be compared to the numerical estimate. Defining x = ln β, y = β E β [ln p(s| λ, M )], and ∆x i = x i − x i−1 , we can develop a trapezoid discrectization of the integral: Now suppose that our estimate for y i is centered atȳ i and has variance σ 2 i . Let us further assume that the errors at different temperatures are uncorrelated: E[(y i −ȳ i )(y j −ȳ j )] = σ 2 i δ ij . Then and The only subtle point to remember when computing this error estimate for thermodynamic integration is that the quantity y i is found as the average of a collection of samples. Thus we need the error on the mean, which for a large number of independent samples N is given by the central limit theorem as σ 2 i = Var[y i ]/N . Since we use a RJMCMC algorithm, our samples are correlated, and to estimate N we first have to thin the chains by an amount determined by their auto-correlation length.
Depending on the spacing of the chains and the length of the runs, the discretization error may be larger or smaller than the statistical error. While it is possible to estimate the discretization error analytically, the expression involves derivatives of the average likelihood which are hard to compute in a reliable fashion. Instead we have developed a novel approach that simultaneously accounts for the discretization error and the statistical error. Taking as input the estimatesȳ i and their variances σ 2 i , we fit a smooth curve to the data and compute the integral using the smooth curve. The curve is defined by a cubic spline, where the number of spline control points and their location can be varied. We define a likelihood ∼ e −χ 2 /2 with where c i are the values of the cubic spline curve at x i . A RJMCMC routine (adapted from the BayesLine algorithm) marginalizes over the number and location of the spline points, and at each iteration an estimate for the integral (30) is computed using the smooth spline curve. The expectation value and variance of these evidence estimates provide the central value and error estimate for the thermodynamic integration approach. We concede that there is an element of irony in the fact that we use a RJMCMC algorithm to estimate error in the thermodynamic integration evidence, but the method is very robust. The error estimate was tested in many ways, including the multiple-run test shown in Figure 8, and tests on known functions. For example, integrating y = 1 + tanh(x) in the interval x = [−1, 2] with 10 evenly spaced samples and errors of σ = 0.01 yielded I = 3.890496 ± 0.006605 ± 0.0032 using a trapezoid integration, where the first error is the statistical error (33) and the second is the discretization error, and I = 3.891174 ± 0.006745 for the RJMCMC approach with smooth splines. Since spline integration is much higher order than trapezoid integration, discretization error is a tiny fraction of the total error in the RJMCMC approach. The statistical error is in very good agreement with the estimate from (33). Note that both estimates include the exact value for the integral, I = 3.891222, within their error ranges. In multiple tests on different functions the RJMCMC approach consistently returned more accurate values than a simple trapezoid integration.

Model Selection Examples
The evidence ratio, or Bayes Factor, between competing models tells us how the prior odds ratio (betting odds) have been changed by the observed data. The Bayes Factor between the signal or glitch model and Gaussian noise is expected to follow the relation [38] where FF is the fitting factor, which is found by maximizing the match (17), and O is the Occam factor -the ratio of the prior to posterior volume, the latter being defined as the region where some large fraction of the posterior weight is concentrated. To leading order, the term ∆ ln O is expected to scale as ∆N log(FF), where ∆N is the difference in the number of parameters between the two models [39]. Empirically we find that the BayesWave fitting factors scale roughly as 1 − FF ∼ 1/SNR, so we expect the log Bayes factor to scale linearly with the signal-to-noise ratio. Figure 9 shows the Bayes factor between the signal model and the Gaussian noise models for a simulated non-spinning black hole merger with masses m 1 = m 2 = 20M , with the expected linear scaling in SNR for SNR > 18. Distinguishing non-Gaussian features in the data, be they instrument glitches or gravitational waves, from the Gaussian noise model is relatively easy. Far more challenging is distinguishing glitches from gravitational waves. While it is extremely rare that glitches occur in every detector in the network at the same time, at the same frequency, and with a similar morphology, the possibility exists, as demonstrated by the near coherent triggers found in time-slides of the LIGO/Virgo data. In the absence of waveform templates, BayesWave uses similar models for glitches and gravitational waves. Both are built from a sum of Morlet-Gabor wavelets, the only difference being that the signal model is defined at the Geocenter, and the signals in each detector are found by projecting onto the network, while the glitch model uses independent sets of wavelets in each detector. When operating on data containing a gravitational wave signal, both models are able to produce similar fitting factors. The signal model achieves a slightly higher fitting factor since both the match and fitting factor increase with SNR, and the network SNR always exceeds the SNR in any one instrument. In most instances, the signal model has a lower dimensionality, and pays a lower Occam penalty as a result. As the number of wavelets needed to reconstruct the waveform increases, or as the number of detectors in the network is increased, the signal model gains a decided advantage over the glitch model. The worst case scenario for the signal model is when there are just two detectors, the signal has a sky location and polarization that puts most of the SNR in a single detector, and a large fraction of the signal power can be captured by a single wavelet. In that case the glitch model can achieve a good fitting factor using 5 parameters, while the signal model is forced to use 9 parameters. BayesWave invariably classifies these signals as glitches. Figure 10 shows signal/glitch and signal/Gaussian noise Bayes factors for a simulated Sine-Gaussian signal using different detector combinations. The simulated signal had single detector signal-to-noise ratios of SNR H1 = 12.3, SNR L1 = 9.2, SNR V = 9.6. In each case, there is clear evidence for a non-Gaussian feature in the data. In the single detector analyses the signal and glitch models have comparable evidence, despite the fact that the signal model carries 4 additional extrinsic parameters. But the signal model pays no penalty for the extra parameters since they are completely unconstrained by data from a single detector and the posterior distribution follows the prior distribution. The signal model is heavily favored over the glitch model in all three two-detector network combinations, and the evidence for the signal model is larger still for the full three-detector network.

Parameter Estimation
Following the pioneering work of Christensen and Meyer [40], the application of Bayesian inference and parameter estimation to gravitational wave astronomy [41,42,43,44,45,46,47,48,49] has mostly been focused on fully-modeled gravitational wave signals described by waveform templates h( λ) that depend on some collection of physical parameters λ. Very recently Bayesian inference has been applied to partially modeled sources such as supernovae [50,51] and Neutron star r-modes [52], but parameter estimation and model selection for completely un-modeled signals has not been tackled until now (Bayesian inference for sky localization of burst signals has been considered using single sine-Gaussian waveforms [53]). When templates are available the posterior distribution for the gravitational wave signal p(h|s) gets mapped to posterior distributions for the model parameters p( λ|s), allowing for a straightforward interpretation of the results, such as producing credible intervals for the masses and spins of a binary system. But what does "parameter estimation" mean in the context of glitches and un-modeled gravitational wave signals? BayesWave provides posterior distributions for the wavelet parameters used to reconstruct the signals and glitches, but these have limited physical meaning.
BayesWave also provides posterior distributions for the reconstructed signals h(t) and glitches g(t), along with posterior distributions for the sky location and polarization of a gravitational wave source. The posterior distributions for the signal and glitch waveforms contain the complete information about the waveforms and their physical properties, but in a form that can be hard to interpret. Figure 11 shows a stretch of simulated whitened data containing a Gaussian-enveloped white-noise burst. The upper panel shows the simulated signal and noise and the median glitch waveform reconstruction. The lower panel shows the posterior distribution for the glitch model, where the colors correspond to the log posterior density ln p(g(t)|s). While providing a useful visualization, plots of the waveform posteriors do not provide a good quantitative understanding of the signal. At each step in the Markov chain, BayesWave can produce the reconstructed waveforms as a function of time, or frequency, or time-frequency. The waveforms can either be whitened, which is what the likelihood actually "sees", or un-whitened, which corresponds to the physical signal. While the un-whitened reconstructions would seem to be the most useful, they tend to be very noisy as the components of the reconstruction with frequencies outside of the "bucket" where the instruments are most sensitive are poorly constrained. Using the reconstructed waveforms it is possible to compute any physical quantity of interest, such as the energy flux, rise and decay times, duration etc at each step of the chain, and produce probability distributions for these quantities. The choice of summary quantities to extract from the waveforms is virtually unlimited, but they will have varying utility. The most useful choices will allow us to connect the reconstructed waveforms to astrophysical models for the signals, or for glitches, to physical models of the detector. It may be that higher dimensional quantities are needed, such as posterior distributions for Bezier curves that track the development of the signal in time-frequency. To illustrate the process of producing posterior distributions for parameters that summarize the physical content of the waveforms, we can use the "waveform moments" proposed by Patrick Sutton. The idea is to use the unit normalized densities where to define moments: The central time t c = t and central frequency f c = f are given by the first moments, while the duration and bandwidth are defined in terms of the variance: These quantities can be computed for both the whitened and un-whitened waveforms. Similar expressions can be defined for glitches. Figure 10 shows posterior distributions for the central time, central frequency, duration and bandwidth (computed using whitened waveforms) for a Gaussian enveloped white noise burst and a black hole merger using a uniform time-frequency prior and the proximity prior. In each case, the proximity prior improves the parameter estimation, which is consistent with the improvements seen in the fitting factors for these systems shown in Figures 5 and 6. These examples illustrate a general trend that we have seen toward systematic biases in the recovery of certain parameters, most notably the duration of the signal. In general we find that BayesWave under-estimates the duration and energy content of the signals. These biases are to be expected given that the match between the recovered and simulated signals is always less than unity. It may be possible to compensate for this bias since the cause is known. On the other hand, the waveform moments are particularly susceptible to such biases, and we are probably better off using more robust summary statistics, such as defining the duration and bandwidth as the intervals containing 90% of the energy.

Applications using LIGO/Virgo data
In addition to tests on simulated data, like those shown in the preceding sections, BayesWave has been extensively tested on data from the S5/6 LIGO and VSR2/3 Virgo science runs. These studies will be presented in full elsewhere, but as a prelude, we show a few highlights from these studies here. One application of BayesWave/BayesLine is to use the noise spectrum and glitch models in concert with standard gravitational wave template based searches and parameter estimation studies in the hope that removing glitches from the data will reduce the number of background events and lead to smaller biases in the extraction of the signal parameters when a detection is made. To illustrate the potential of this approach, we took a stretch of moderately glitchy data from the LIGO S6 science run and produced whitened scalograms with and without glitch fitting. Histograms of the wavelet amplitudes at various discrete wavelet resolutions are shown in Figure 13 (here we used the binary Meyer wavelet transformation developed by S. Klimenko [30]). The discrete wavelet transform has pixels with bandwidth ∆f , duration ∆t and area ∆f ∆t = 1. The histograms are shown for a variety of resolutions ∆f as glitches come with a variety of morphologies that are picked up more clearly in some resolutions than others [36]. Note that the BayesWave residuals, with glitches regressed, are Gaussian distributed at all resolutions. In addition to performing parameter estimation for signals and glitches, BayesWave can be used to help decide if candidate events are gravitational wave signals or (un)lucky coincidences between noise events in different detectors. The main burst search algorithm used by the LIGO/Virgo collaborations, coherent WaveBurst [32], performs a constrained maximum likelihood reconstruction of gravitational wave signals, supplemented by additional cuts on the "network correlation coefficient", the "effective correlated amplitude" and the "network energy disbalance", which are designed to reject accidental coincidences from glitches. These cuts are tuned using simulated signals and time-slides of the data, the former to measure the impact on the false dismisal rate, the latter to measure the the impact on the false alarm rate. The most significant events found in time-slides after cuts and data quality vetoes have been applied determine the overall sensitivity level for the coherent WaveBurst search. Figure 14 shows examples of the waveform reconstruction for the 4 km LIGO detectors at Hanford and Livingston taken from randomly chosen iterations of the Markov chains. In each case the signal and glitch reconstructions are very similar, which is why these spurious coincidences masqueraded as signals in the coherent WaveBurst analysis. Going from the top panel to the bottom panel, coherent WaveBurst reported network signal-to-noise ratios and "network correlation coefficients": Upper -SNR = 64.4, cc = 0.84; Middle -SNR = 20.4, cc = 0.93; Lower -SNR = 18.7, cc = 0.88. The BayesWave algorithm does not compute quite the same quantities, so instead we quote the maximum signalto-noise ratios and Bayes Factors between the signal, glitch and noise models: Upper -SNR max = 88.7, ln BF s/n = 1872.9 ± 0.5, ln BF g/s = 24.5 ± 3; Middle -SNR max = 19.5, ln BF s/n = 63.8±0.5, ln BF g/s = −6.9±2; Lower -SNR max = 18.9, ln BF s/n = 53.7±0.8, ln BF g/s = −2.3 ± 3. In each case BayesWave found overwhelming evidence for a non-Gaussian feature in the data. For the event in the upper panel the evidence decisively favored the glitch model over the signal model, while for the event in the middle panel the evidence slightly favored the signal model over the glitch model. For the event in the lower panel the evidence is inconclusive. The fact that BayesWave found the event in the middle panel to be most signal-like is consistent with this event having the highest network correlation coefficient in the coherent WaveBurst analysis -sometimes sufficiently similar coincident glitches in two detectors can fool any Burst search. As a final example, we show in Figure 15 sky maps for a notable blind signal injection made on September 16, 2010 [54], that has variously been dubbed "GW100916" for the date of discovery [55] or "the Big Dog", as the initial sky reconstruction placed the event in the direction of Canis Major. The panel on the left shows the log posterior density derived using inspiral templates, while the panel on the right used the BayesWave wavelet model. The two sky maps are broadly consistent. The BayesWave sky map is more diffuse than the template derived sky map, as is to be expected since the wavelet model captures less of the signal power than the template model.

Conclusion and Outlook
We have presented a new approach to gravitational wave data analysis that puts as much emphasis on modeling the instrument noise as it does on modeling the signals. Our approach follows the mantra "model everything and let the data sort it out", where crucially, the data gets to decide on the complexity of the model. The BayesWave algorithm provides a flexible framework for analyzing gravitational wave data. The signal and glitch models can be changed by modifying the priors. Targeted analyses can easily be developed for particular classes of signals that build in even partial knowledge about the waveforms. Work is currently underway to develop targeted analyses for spinning binary mergers and high eccentricity binaries that use parameterized priors based on the time-frequency evolution of these systems. As more detectors join the world-wide network we will be able to do away with the restriction that the signals are elliptically polarized and independently infer the two polarization states h + and h × .
Currently, parameter estimation for gravitational wave bursts and instrument glitches is in its infancy. There is much room for improvement in the development of quantities that meaningfully summarize the physical characteristics of the reconstructed waveforms and relate these to models for the astrophysical or electro-mechanical processes that gave rise to the signals or glitches.