Novel Method for Incorporating Model Uncertainties into Gravitational Wave Parameter Estimates
Posterior distributions on parameters computed from experimental data using Bayesian techniques are only as accurate as the models used to construct them. In many applications these models are incomplete, which both reduces the prospects of detection and leads to a systematic error in the parameter estimates. In the analysis of data from gravitational wave detectors, for example, accurate waveform templates can be computed using numerical methods, but the prohibitive cost of these simulations means this can only be done for a small handful of parameters. In this work a novel method to fold model uncertainties into data analysis is proposed; the waveform uncertainty is analytically marginalised over using with a prior distribution constructed by using Gaussian process regression to interpolate the waveform difference from a small training set of accurate templates. The method is well motivated, easy to implement, and no more computationally expensive than standard techniques. The new method is shown to perform extremely well when applied to a toy problem. While we use the application to gravitational wave data analysis to motivate and illustrate the technique, it can be applied in any context where model uncertainties exist.
Bayesian techniques are widely used to draw inferences from experimental data. These rely on having a model for the signal and inferences are only as accurate as the underlying model. There are two potential problems with inaccurate models — detection and parameter estimation. If the model is a poor representation of any signal present then the integral of the posterior over parameter space (the evidence, a common model selection metric) will be underestimated. This decreases the probability that the signal will be detected and potentially increases the chance of incorrect model selection; this is the detection problem. The parameter estimation problem is that even when a detection is made the position of the peak of the posterior evaluated using an inaccurate model may be shifted relative to the true position, possibly by much more than the statistical error arising from random noise. This shift means that there will be a systematic error in the parameter estimates. The statistical error decreases with increasing signal amplitude, but the systematic error remains constant Cutler and Vallisneri (2007) and is therefore most important for the loudest sources. In this work we describe how model inaccuracies can be included and marginalised over within a Bayesian framework using Gaussian process regression (GPR). This approach is straightforward to implement, incurs no additional computational cost and naturally solves both of these problems.
One situation in which model inaccuracies are known to be present is in the analysis of data from gravitational wave (GW) detectors. Over the next few years kilometre-scale laser interferometers (LIGO, Virgo) will start to take data in Advanced configurations which are expected to make detections of GWs from merging compact binaries at the rate of a few tens of events per year Abadie et al. (2010). Observations of nanohertz and millihertz GW sources using pulsar timing arrays and space-based interferometers should follow in the coming decades. These observations have the potential to transform our understanding of the astrophysics of compact objects, but only if the source parameters, e.g., the masses and spins of the compact objects, can be accurately estimated. Inference will use templates of the GWs emitted in these systems. Recent advances in numerical relativity (NR; Pretorius (2005)) have made it possible to accurately compute GW templates. These simulations are prohibitively computationally expensive, however, and inference will therefore rely on approximate templates, including post-Newtonian (PN) models (for a recent review see Blanchet (2014)), numerical “kludge” techniques Babak et al. (2007) or the effective-one-body approach (see, for example, Damour (2012)). Use of these approximations will introduce systematic errors of the type described above and these have been shown to be potentially significant for both ground (Canitrot, 2001) and space-based detectors (Cutler and Vallisneri, 2007). These systematic errors must be accounted for to enable correct astrophysical inference from near-future GW observations.
The best current GW models, such as EOBNR Pan et al. (2011) or IMRPhenomC Santamaría et al. (2010), are constructed by using NR simulations to fit the values of free parameters in extended analytic models. Systematic errors are still present for these models, albeit reduced relative to pure analytic models. The approach we take to accounting for model uncertainties also exploits the fact that the model is known more precisely (although not necessarily perfectly) where NR simulations have been performed, but without needing to include arbitrary free parameters in the model.
The guiding philosophy of our approach is that of Bayesian analysis; the uncertainty in the waveform model is marginalised over, folding the information available from the NR waveforms into the prior probability distribution. This is accomplished by using GPR to interpolate the difference between the accurate waveform (here taken to be the NR waveform at the points where these are available, although using EOBNR or IMRPhenomC would also be possible) and the approximate waveform across parameter space. GPR returns a probability distribution for the waveform difference. This distribution is used as a prior to analytically marginalise the likelihood over the waveform uncertainty. The marginalised likelihood can be used in place of the usual likelihood in any existing search algorithm. The peak of the marginalised likelihood is shifted and its width broadened to account for systematic errors, thus solving the parameter estimation problem. In addition, using the evidence associated with the marginalised likelihood as a detection statistic solves the detection problem.
In Section II we demonstrate how to use GPR to construct a prior probability distribution on the waveform difference and how to marginalise the likelihood analytically over this distribution. In Section III we demonstrate the efficacy of the marginalised likelihood with a toy numerical problem — the extraction of the chirp mass of a high order PN waveform using a lower order PN waveform. We conclude with a discussion in Section IV.
Ii Modified likelihood
Assume that the GW source is fully characterised by parameters . Assume further that we have the ability to calculate the accurate (although not necessarily exact) waveform templates, (for notational simplicity the dependence of the templates on time is suppressed throughout). These exact templates are referred to as NR templates, and are computationally expensive to produce. Additionally, we can compute approximate waveform templates, . These approximate waveform templates are referred to as PN templates, and are computationally cheap. The templates are related by the waveform difference, :
To perform parameter estimation we must calculate the posterior distribution, , from the observed data . From Bayes theorem this is proportional to the product of the likelihood, , and the prior, . For a detector with stationary, Gaussian noise with power spectral density , the true likelihood is
For simplicity (although it is not necessary), we will assume throughout that is flat within the relevant region of parameter space. The posterior is proportional to the likelihood and the evidence becomes
In practice it is impossible to sample from the distribution in Eq. 2 because it is prohibitively expensive to calculate the NR waveforms; instead, we must use the PN waveforms to calculate an approximate likelihood,
The natural way to reduce the error in Eq. 5 is to construct improved (and inevitably more expensive) approximants with smaller . Instead, the proposal of this paper is to replace with a new likelihood which accounts for the uncertainty in the waveforms. The alternative likelihood is
This new likelihood has marginalised over the uncertainty in the waveform difference using the (as yet unspecified) prior . The prior on the waveform difference should include the information available from the limited number of available NR waveforms and must encode our expectation that the PN waveforms are accurate at early times (or equivalently low frequencies) when the orbiting bodies are well separated, but gradually become inaccurate as the bodies inspiral. At most points in parameter space a NR waveform will not be available, and so it is necessary to interpolate the waveform difference across parameter space and simultaneously account for the error this introduces. It would seem that the problem rapidly becomes complex, and even if a suitable prior could be constructed the computational time needed to evaluate would make it impractical in most contexts. Fortunately, the technique of GPR provides a natural way to interpolate the waveform differences across parameter space, incorporating all necessary prior information. GPR also has the additional property, arising from its construction, that it returns an expression for which is a Gaussian in . Since the exponential factor in Eq. II is also Gaussian the integral may be evaluated analytically. This gives an expression for which can be evaluated in the same computational time as . Henceforth, for brevity, the will be suppressed in all likelihoods; e.g. .
ii.1 Gaussian process regression
We will briefly summarise the key results from GPR here and refer the reader to standard texts for details (MacKay, 2003; Rasmussen and Williams, 2006). A useful way of thinking about the technique of GPR is as a probabilistic interpolation algorithm. Assume that we have access to NR waveforms at a few values of the parameters, and can cheaply compute PN waveforms at the same parameter values. Our training set are the waveform differences
Assume that all waveform differences in the , and one additional difference at parameters , are drawn from a multivariate Gaussian with zero mean and covariance :
where the matrix/vector/scalar
are defined in terms of a covariance function, . Specifying the covariance function is central to GPR as it encodes our expectations about the properties of the function being interpolated. Here we use the squared exponential covariance function (sum over )
This is a very common choice of covariance function (Rasmussen and Williams, 2006) and defines a stationary, smooth Gaussian Process (GP). Other covariance functions could be considered and their applicability verified by reserving a subset of to check against the interpolation. In Eq. 10, a scale and a (constant) metric defining a modulus in parameter space has been defined; these hyperparameters are determined by maximising the evidence for , which is discussed in Sec. III. These choices provide a simple, data-driven and quick-to-evaluate method which is ideal for parameter estimation, but the (dis)advantages of using other covariance functions should be explored in future work. If the available NR waveforms contain some uncertainty then this may be included by adding a diagonal matrix to Eq. 10, where the element is the fractional uncertainty in the NR simulation at times the signal-to-noise. The choice of a zero mean GP reflects our prior belief in the validity of the PN waveforms. Our uncertainty in the PN waveforms is encoded in the matrix via the choice of the covariance function. The conditional probability of the unknown waveform difference given the known differences in is then
ii.2 Marginalised Likelihood
Furnished with , the marginalised likelihood may now be evaluated. The integrand in Eq. II is the product of two Gaussians and may be evaluated analytically to give
The best-fit waveform has shifted by ; this is the best estimate of the waveform difference returned by the GPR. The likelihood distribution is also broadened, since .
The likelihood in Eq. II.2 takes no appreciable extra time to evaluate compared to the likelihood in Eq. 5, because the most expensive step is evaluating the inner product, which is common to both. Computing and requires inverting the matrix ; however this only depends on quantities in and so may be performed offline (i.e. before the evaluation of the posterior begins) and saved in memory for subsequent evaluations.
Iii Numerical implementation
To illustrate this method a simple search using non-spinning PN waveforms was performed. In place of the NR waveforms 3.5PN order waveforms were used, while the approximant was a 3PN waveform. A phase shift was included in the approximate waveforms to match them with the exact waveforms at the start, to illustrate the fact that GPR naturally accounts for waveform differences that are small at early times and then grow. The search was restricted to one parameter, the chirp mass . A signal was with intrinsic parameters and extrinsic parameters , and a short sample of the waveform between and was analysed using an analytic fit to the LIGO noise curve (from table-1 of Sathyaprakash and Schutz (2009)). For this toy problem no noise was injected into the mock data; this was so the peak of the exact likelihood was located on the true parameters and the systematic error which is the subject of this paper is clearly visible as a shift from this position. If noise was added the peak would shift by an amount consistent with the width of the posteriors but the systematic error would remain unchanged.
For this toy problem, was taken to consist of the waveform difference computed at five chirp masses, . The matrix from Eq. 8 was calculated an inverted offline for use during later likelihood evaluations.
was maximised for and where . Fixing the hyperparameters by maximising the evidence is a convenient heuristic, but other approaches could be considered. In general the hyperparameters could be marginalised over, treating them as nuisance parameters in a Bayesian search. Our approach has the advantage of performing all the calculation at the “offline” stage, adding no extra cost to the “online” search.
The PN approach breaks down near merger and hence there is poor agreement between the different PN orders; however, this just serves here to make the systematic errors which are the subject of this paper more apparent. Using an inspiral-merger-ringdown approximant, it would be expected that much smaller waveform differences, and hence smaller , would be obtained.
GPR can now be used to interpolate the waveform difference across the range of used in the search; plotted in Fig. 1 are examples of the interpolated waveform difference, and plotted in Fig. 2 is the error in the GPR analysis as a function of . It can be seen from Fig. 1 that well outside the GPR returns the zero function for , whilst within the range of it returns a smooth interpolation of the waveform difference. From Fig. 2 it can be seen that the GPR returns a large error outside of and a small error near a point in .
Since we were working in only one dimension, the search was performed using a template grid. The resulting likelihood surfaces for , and are shown in Fig. 3. The left-hand panel of Fig. 3 was the most computationally expensive to produce due to the higher order PN waveforms; as anticipated in Sec. II.2, there was no measurable difference in computational time between the centre and right-hand panels. Fig. 3 shows that the approximate likelihood is shifted substantially with respect to the exact likelihood (the parameter estimation problem) and the evidence is dramatically reduced (the detection problem). The marginalised likelihood in the right-hand panel has addressed both of these issues.
The method outlined here does not assume that the waveform differences, , are small. In addition, although we have ignored errors in the training set waveforms for our example, the method can account for errors in the accurate waveforms. We expect the best performance when both the waveform difference and the errors in the accurate waveforms are small. However, even if this were not the case we would still expect the modified likelihood to outperform the approximate likelihood.
In this paper, a novel approach to tackling the twin problems of GW detection and parameter estimation with imperfect templates is proposed. The technique involves replacing the likelihood used in standard approaches with a modified likelihood marginalised over the uncertainty in the waveform. GPR is used to interpolate the waveform difference from a training set of NR templates, and this provides the prior for the marginalisation. There have been previous attempts to improve the prospects for GW parameter estimation by interpolating accurate waveforms Smith et al. (2013); the current approach of using GPs has the advantages of being non-parametric and automatically accounting for errors introduced by the interpolation. The resulting likelihood may be used in the same manner as the standard likelihood, and it overcomes both detection and parameter estimation problems. The new technique is exceedingly easy to implement into existing algorithms as only simple modifications to the likelihood are required. The GPR allowed us to incorporate prior beliefs about the accuracy of the templates in a non-parametric manner, and perform the marginalisation analytically. There is also ongoing effort into overcoming model uncertainties using parametric techniques Berry et al. (2014). The method proposed here marginalises over an unknown part of the signal to obtain the likelihood, this is similar to the marginalisation of the timing model in the context of pulsar timing array data analysis (see van Haasteren and Vallisneri (2014) for a recent treatment using GPs).
In our example the marginalised likelihood was shown to perform significantly better than the usual likelihood for both detection and parameter estimation. Crucially, the marginalised likelihood is no slower to evaluate and only a modest amount of offline calculation is required. For these reasons it is anticipated that the marginalised likelihood will be useful for future GW searches.
In this paper the focus has been on inaccuracies that arise from difficulties in building accurate models; however the method in this paper could be adapted to marginalise over any source of uncertainty in the signal assuming enough prior information is available to form a training set. Examples of errors which could be addressed in this manner include the calibration error (a frequency dependent amplitude and phase error in the signal returned from the detector Vitale et al. (2012)), the stealth bias (a systematic error associated with a deviation from GR, before the deviation becomes detectable Vallisneri and Yunes (2013)), and the error from neglecting certain physical phenomena (e.g. the presence of an accretion disk in a compact binary).
This approach can also be used to identify local maxima of the GPR error estimate, which could guide the NR community to regions of parameter space where new simulations are most needed. Moreover, it should have applications beyond GW data analysis. GPs are already commonly used in engineering, e.g., Ko and Fox (2009); Schwaighofer et al. (2004), but these techniques apply to any problem where construction of detailed models is expensive and inference relies on approximations.
The next step is to implement this technique in a higher dimensional parameter space using more realistic waveform models. This will not only provide an assessment of how significant the systematic errors from standard searches of near future advanced detector data will be, but will provide a marginalised likelihood that will give correct parameter inferences no significant additional computational cost.
Acknowledgements.We thank Christopher Berry for helpful discussions. CM is supported by the STFC. JG’s work is supported by the Royal Society.
- Cutler and Vallisneri (2007) C. Cutler and M. Vallisneri, Phys. Rev. D 76, 104018 (2007), URL http://link.aps.org/doi/10.1103/PhysRevD.76.104018.
- Abadie et al. (2010) J. Abadie, B. P. Abbott, R. Abbott, M. Abernathy, T. Accadia, F. Acernese, C. Adams, R. Adhikari, P. Ajith, B. Allen, et al., Classical and Quantum Gravity 27, 173001 (2010), eprint 1003.2480.
- Pretorius (2005) F. Pretorius, Phys. Rev. Lett. 95, 121101 (2005), URL http://link.aps.org/doi/10.1103/PhysRevLett.95.121101.
- Blanchet (2014) L. Blanchet, Living Reviews in Relativity 17 (2014), URL http://www.livingreviews.org/lrr-2014-2.
- Babak et al. (2007) S. Babak, H. Fang, J. R. Gair, K. Glampedakis, and S. A. Hughes, Phys. Rev. D 75, 024005 (2007), eprint gr-qc/0607007.
- Damour (2012) T. Damour, ArXiv e-prints (2012), eprint 1212.3169.
- Canitrot (2001) P. Canitrot, Phys. Rev. D 63, 082005 (2001).
- Pan et al. (2011) Y. Pan, A. Buonanno, M. Boyle, L. T. Buchman, L. E. Kidder, H. P. Pfeiffer, and M. A. Scheel, Phys. Rev. D 84, 124052 (2011), eprint 1106.1021.
- Santamaría et al. (2010) L. Santamaría, F. Ohme, P. Ajith, B. Brügmann, N. Dorband, M. Hannam, S. Husa, P. Mösta, D. Pollney, C. Reisswig, et al., Phys. Rev. D 82, 064016 (2010), eprint 1005.3306.
- MacKay (2003) D. MacKay, Information Theory, Inference, and Learning Algorithms (Cambridge University Press, 2003).
- Rasmussen and Williams (2006) C. Rasmussen and C. Williams, Gaussian Processes for Machine Learning (The MIT Press, 2006).
- Sathyaprakash and Schutz (2009) B. Sathyaprakash and B. F. Schutz, Living Reviews in Relativity 12 (2009), URL http://www.livingreviews.org/lrr-2009-2.
- Smith et al. (2013) R. J. E. Smith, K. Cannon, C. Hanna, D. Keppel, and I. Mandel, Phys. Rev. D 87, 122002 (2013), eprint 1211.1254.
- Berry et al. (2014) C. P. L. Berry, W. M. Farr, J. R. Gair, and I. Mandel, in preperation (2014).
- van Haasteren and Vallisneri (2014) R. van Haasteren and M. Vallisneri, ArXiv e-prints (2014), eprint 1407.1838.
- Vitale et al. (2012) S. Vitale, W. Del Pozzo, T. G. F. Li, C. Van Den Broeck, I. Mandel, B. Aylott, and J. Veitch, Phys. Rev. D 85, 064034 (2012), URL http://link.aps.org/doi/10.1103/PhysRevD.85.064034.
- Vallisneri and Yunes (2013) M. Vallisneri and N. Yunes, Phys. Rev. D 87, 102002 (2013), eprint 1301.2627.
- Ko and Fox (2009) J. Ko and D. Fox, Autonomous Robots 27, 75 (2009), ISSN 0929-5593, URL http://dx.doi.org/10.1007/s10514-009-9119-x.
- Schwaighofer et al. (2004) A. Schwaighofer, M. Grigoraş, V. Tresp, and C. Hoffmann, in Advances in Neural Information Processing Systems 16 (MIT Press, 2004), pp. 579–586, URL http://research.microsoft.com/apps/pubs/default.aspx?id=74491%.