Density estimation with Gaussian processes for gravitational wave posteriors

The properties of black hole and neutron-star binaries are extracted from gravitational waves (GW) signals using Bayesian inference. This involves evaluating a multidimensional posterior probability function with stochastic sampling. The marginal probability distributions of the samples are sometimes interpolated with methods such as kernel density estimators. Since most post-processing analysis within the ﬁeld is based on these parameter estimation products, interpolation accuracy of the marginals is essential. In this work, we propose a new method combining histograms and Gaussian processes (GPs) as an alternative technique to ﬁt arbitrary combinations of samples from the source parameters. This method comes with several advantages such as ﬂexible interpolation of non-Gaussian correlations, Bayesian estimate of uncertainty, and efﬁcient resampling with Hamiltonian Monte Carlo.


INTRODUCTION
The first detection of gravitational waves (GW) in 2015 (Abbott et al. 2016a) sparked a new era of Astronomy.Several years on from that event the number of detected gravitational waves keeps increasing and within this decade we expect to observe O(10 3 ) signals (Abbott et al. 2020b) from compact binaries coalescences (CBCs).This huge progress brings with it the challenge of efficiently analysing a large number of events.To address these computational challenges, machine learning techniques are being increasingly investigated within the field of gravitational-wave physics (Cuoco et al. 2020).Many studies have focused on speeding parameter estimation of the source parameters of the signals with various techniques, such as deep learning (George & Huerta 2018), variational autoencoders (Gabbard et al. 2019) and autoregressive neural flows (Green et al. 2020).Other work has focused on combining detection and parameter estimation with deep neural networks (Fan et al. 2019) as well as using neural networks to rapidly generate surrogate waveforms (Khan & Green 2020;Chua et al. 2019).
While the research efforts to speed up or completely revolutionise parameter estimation are ongoing, the issue of how to effectively deal with a large number of results from different events remains.In particular, how to streamline the analysis of the results, while maintaining accuracy.In this work, we demonstrate the efficiency and usefulness of using Gaussian processes (GP) for post-processing parameter estimation results of CBCs.Applications of GPs in the field of gravitational waves span a wide range of use-cases, such as marginalising waveform errors (Moore et al. 2016), regression of analytical waveforms (Setyawati et al. 2020), predictions of population synthesis simulations (Barrett et al. 2016), hierarchical population inference (Taylor & Gerosa 2018) and Equation of State (EOS) calculations (Landry & Essick 2019).They have also been exploited for the development of fast parameter estimation with RIFT sampler (Lange et al. 2018).
Here we exploit GPs to estimate probability density functions (PDFs) from parameter estimation of gravitational-wave signals.Non-parametric density estimation from a finite set of samples is an active research field in machine learning and statistics (Murray et al. 2008;Papamakarios et al. 2017;Wang & Scott 2019).
For most GW analysis, histograms are usually the preferred estimators to visualise the marginal posterior PDFs and to avoid over-smoothing sharp features, but often are not convenient for post-processing analyses such as population inference.These sorts of analyses either re-weight the posterior samples directly (Abbott et al. 2020a) or need to estimate a continuous representation of the gravitationalwave posterior density surface.Several density estimation methods such as Dirichlet process' (Del Pozzo et al. 2018), Gaussian Mixture Models (Talbot & Thrane 2020) and others have been employed to address this problem specifically for gravitational-waves.As well as these, A closely related method to GPs (Kanagawa et al. 2018), Gaussian Kernel Density Estimators (KDEs) are sometimes employed in gravitational waves' analyses (Pitkin et al. 2018;Pang et al. 2020;Lynch et al. 2017).
These KDEs are often effective but they assume correlations between parameters to be linear and smooth, making this method sometimes limited in flexibility.There exist many variations of the KDE algorithm to take into account specific interpolations problems, but there isn't a single implementation that is guaranteed to be robust against all pos-sibilities.A specific KDE implementation might solve an issue in one case and be the cause of some inaccuracies in another Wand & Jones (1994).
We implement a single technique that can interpolate arbitrary multi-dimensional slices in parameter space, which can handle both simple and difficult space morphology, such as sharp bounds and non-Gaussian correlations.Our modelling tool is based on the histogram density estimate, combining the histogram's accurate treatment of the samples' features with the predictive capabilities of GPs.An additional advantage of this technique is that it can provide a Bayesian measure of uncertainty from the finite (and sometimes small) number of samples for post-processing analysis.This measure of model uncertainty could then be incorporated into any analysis where the marginalised posterior density is used.
In Sec. 2 we describe our density estimation technique in the context of gravitational-wave parameter estimation and machine learning.We propose a series of example applications in Sec. 3, which allows us to discuss the advantageous features of our method.Finally, in Sec. 4 we summarise our findings and suggest future extensions of this work.

METHODS
In this section, we introduce the mathematical framework of the techniques discussed.In subsection 2.1 we discuss the Bayesian inference problem for gravitational waves and the density estimation techniques currently employed in the field.In subsection 2.2, we outline the fundamentals of GPs and their interpretation for interpolating a posterior density surface.We then describe the details of our GP implementation and how to model probability densities from parameter estimation.

Bayesian inference and density estimators
We describe the gravitational-wave data in the detector d as the sum of a waveform model h( θ) and a combination of instrumental noise, which we assume to be Gaussian.The probability of observing data parametrised by θ = {θ1, θ2, ..., θN }, can be defined as the probability that r = d − h is the realisation of the instrumental noise.This likelihood can be written as : where a|b denotes the inner product between two waveforms a and b and is defined as (Cutler & Flanagan 1994), where Sn(f ) is the one-sided power-spectral density (PSD) and ã denotes the Fourier transform of the gravitational waveform a.By choosing astrophysically motivated priors over the model parameters, we can use the Bayesian framework to calculate the posterior probability distribution for the source parameters (Thrane & Talbot 2019) The posterior probability is generally a 15-dimensional surface for a circular binary black hole (BBH) merger but can be 17-dimensional in the case of a binary neutron star (BNS) merger.The dimensionality depends on the physical parameters describing the signals.Generally, these are distinguished between extrinsic parameters, such as sky localisation, and intrinsic parameters, such as the masses of the sources.
Posterior distributions for specific parameters can then be found by marginalising over all other parameters, p(θi|d) = p( θ|d)dθ1...dθi−1dθi+1...dθN . (4) The posterior is however generally intractable and therefore must be evaluated via stochastic methods such as Markov Chain Monte Carlo (MCMC) and nested sampling, these are implemented (and specifically tuned for the gravitational wave problem) in Bayesian inference packages such as LAL-Inference (Veitch et al. 2015) and bilby (Ashton et al. 2019).

Definition and interpretation
GPs are interpolation methods with a probabilistic interpretation, they are built on a Bayesian philosophy, which allows you to update your beliefs based on new observations.The process can be understood as an infinite distribution over functions, such that any finite collection of these has a Gaussian distribution.As data is observed, the GP is conditioned and the range of possible functions that can explain the observations is constrained.As such a GP is defined by a mean, which represents the expectation value for the best fitting function, and by a covariance matrix, called a kernel, which measures the correlations between observations (Williams & Rasmussen 2006).In the absence of observations, the GP predictions will revert to a prior mean function, which is usually chosen to be zero, and which properties are determined by the kernel architecture.Mathematically this is written as: where the mean and covariance are denoted as: We can then model a surface y conditioned on our observations as: where, in this application, the dimensionality of f will depend on how many parameters p(θi|d) has been marginalised over.
The non-parametric nature of GPs makes this technique flexible, but it can be computationally expensive as the whole training set needs to be taken into account at each prediction.The standard implementation has O(N 3 ) computations and O(N 2 ) storage, this then becomes prohibitive for ∼ 10k data observations or more.To tackle this issue it is common to use sparse inference methods, which approximate the conditioning of the GP over a set of M << N 'inducing' points.The evaluation over the inducing points M is then much cheaper than for an 'exact' GP resulting in O(N M 2 ) computations rather than O(N 3 ) (Quiñonero-Candela & Rasmussen 2005; Hensman et al. 2013).As well as sparse methods one can exploit multi-GPU parallelization and methods like linear conjugate gradients to distribute the kernel matrix evaluations which then allows for exact inference to be performed on a short time scale (Wang et al. 2019).In this work, however, we find that sparse approximations were accurate enough to effectively model the marginalised posterior surfaces that we were interested in.Moreover, once a GP has been 'trained' over the data, it is possible to draw infinitely many function realisations from it without recomputing the expensive covariance matrix.
A recognised advantage of GPs is reliable uncertainty estimate when making predictions over unseen data.In this application, we are not interested in predicting the value of the posterior in unexplored regions of the parameter space, but only in generating a faithful model where we have posterior samples.In regions within the space of parameters, the GP variance depends on our choice of training points, which is useful to assess the accuracy of our density estimation.In terms of uncertainty estimation this can be explained as our model having very low epistemic uncertainty everywhere, we then seek to estimate the aleatoric uncertainty due to our model fit around the random fluctuations in the histogram densities which are used to train the GP.

Model construction
In this application, we want to use a GP to estimate the marginalised posterior density for any subset of parameters.We train our GP using the normalised histogram counts over a grid of points, i.e. the centroids of the histogram bins, that cover the marginalised parameter space.We then fit our GP to this discrete set of points to generate a continuous representation of the surface.
We employ TensorFlow and GPFlow to implement our GP modeling infrastructure, which includes two inference schemes: exact inference for 1-2 dimensional problems (O ∼ 1000 samples) and sparse inference for higher dimensionality due to computational costs.As well as a difference in the inference scheme, when extending this method to higher dimensions, our choice of training data changes.When creating the grid over four dimensions we find that the typical set has a volume of O(1%),(this is a common problem associated with the curse of dimensionality).We, therefore, choose to discard the empty bins and encode our knowledge of these points through the choice of prior over our GP.
Since the model is constructed with converged posterior samples, there is no probability support where the histogram bins are empty.To encode this, we set the mean of the GP to be equal to zero, such that far away from the training data the model will have a high variance but a mean of zero.
To estimate the density for a given region of parameter space we then simply evaluate the GP at those parameters, i.e.
Due to bounded priors (e.g at mass ratio m2/m1 := q = 1), the posterior surface often presents sharp discontinuities and therefore the surface is only piece-wise continuous.GPs are in principle flexible enough to model any surface including piece-wise continuous ones, however, we found in practice that it is more favourable to decompose our density function into two components, one smooth, continuous function, and one step function.We do this by multiplying the density and our GP estimate by a step function, which is zero at any discontinuities and 1 otherwise.
Multiplying by this step function is then analogous to imposing a prior over our posterior surface, i.e. it allows us to rewrite the equation 7 as We are free to encode our knowledge in this way and perform the decomposition as we do not change the original posterior surface that we would like to model in any way.This enhances the robustness of the model against all discontinuities, including artificial cuts in parameter space that might be required for post-processing analysis.
The variance of the GP depends on the kernel, but also on the noise variance parameter of the likelihood.Usually, the noise variance is given by a single number, i.e. homoskedastic noise, which reflects the random fluctuations of the posterior samples.In low-dimensional examples, where we employ an exact inference scheme, we can assign multiple values to the noise variance, i.e. heteroskedastic noise (McHutchon & Rasmussen 2011).In such instances, we are then able to propagate the error from the histogram on the density estimate, which is simply given by the Poisson noise in each bin σ bin ∼ √ Ncounts.Incorporating heteroskedastic errors within a sparse inference scheme is an area of current research in the field of machine learning (Liu et al. 2020).
It is common practice to build an interpolation of a posterior surface in order to draw more samples from it.As our model is implemented in TensorFlow we can quickly draw more samples from the marginalised posteriors using the many samplers available in the package library, such as Hamiltonian Monte Carlo (HMC) (Betancourt 2017).Further technical details which had a significant impact on our model accuracy, such as our data pre-processing scheme and our choice of kernel design are included in Appendix A.2.

RESULTS
In this section, we present our model and a series of example applications for gravitational waves.In Sec.3.1 we illustrate the method on a simple 1-dimensional analytical example.In Sec.3.2 we show examples of common post-processing applications for our density estimation tool.Finally, we discuss our treatment of GP model uncertainty and how we propagate it to produce uncertainty on the marginalised posterior distributions.

Analytical 1D example
Our proposed GP modelling technique is by construction flexible and robust against all distribution morphologies.To illustrate this, we construct a non-trivial 1dimensional example: an inverse gamma function f (x, α) = x −α−1 Γ(a) exp(− 1 x ), with α = 2 and a sharp bound at x = 0.75.In Figure 1 we show our GP model mean prediction and uncertainty, compared to a Gaussian KDE from scipy.statsVirtanen et al. (2020) and two KDE transformations implemented in PESummary (Hoy & Raymond 2020), a commonly used post-processing package in gravitationalwave astronomy.The reflection and transform KDEs, are examples of augmentations on the standard (Gaussian) KDE, and are generally used to model difficult features introduced at the boundaries of posterior distributions.Both of these improvements to the standard KDE apply a transformation at the boundary which implicitly assumes some distributional features (see (Hoy & Raymond 2020) for more details).A Gaussian Process on the other hand makes no assumptions about the distributional shape and can in principle fit any distribution.
We show an example in Fig 1 where our GP is able to well model the posterior and the reflection KDE provides a better fit than the other KDE methods.The transform KDE is more sensitive to noisy features in the samples and can present artifacts, while the Gaussian KDE over-smooths the sharp cut at 0.75.Following this illustrative example, there are others where the reflection KDE is less appropriate.This example was chosen to highlight a case where the choice of KDE is important to fit the distribution well.While synthetic and not representative, it does illustrate features that can and do happen in gravitational-wave astronomy when analysing posteriors.In examples such as this our GP model provides an alternative method to KDEs, requires less handtuning, and also provides a Bayesian estimate of the error on the density estimate, as propagated from the histogram errors.

GW Applications
We now look at a few important post-processing problems in gravitational-wave astrophysics.The training time required to generate the models presented in this section is of the order O(2mins), with variations due to the dimensionality of the surface and to the inference scheme employed.To assess the quality of the model in more than one dimension we decide to re-sample the surrogate surface and compare the new samples to the original set, part of which has been used for training.All samples used in the following sections are taken from the Bilby GWTC-1 catalog (Romero-Shaw et al. 2020).

Catalogue of gravitational-wave properties
Gravitational-wave detection parameters can be distinguished between those intrinsic to the sources, such as the component masses, and those extrinsic to them, such as the sky location.Interpolating the marginal posteriors of combinations of these parameters is often necessary for postprocessing.The following example illustrates a simple case where one can use a GP to interpolate the intrinsic parameters for a given detection.In practice, this could then be repeated for entire GW catalogues so that these interpolated posterior surfaces are then combined for population inferences on the sources of GWs.
For this example, we interpolate the marginal posterior distribution of the intrinsic parameters of the first BBH detection GW150914 (Abbott et al. 2016b), parametrised as follows: chirp mass M, mass ratio q = m2/m1 (where m1 > m2), effective inspiral spin component χ eff and effective precession spin χp, defined by the spin components that lie in the orbital plane (Schmidt et al. 2015).In Figure 2 we compare the marginal distributions sampled from our GP model to the original PE samples.We can visually assess that the correlations between parameters are accurately reconstructed as the 50% and 90% contour lines overlap for each pair of parameters.In Table 1 we report the mean and 90% confidence intervals of the samples drawn from our model and which we find in agreement to the values from the original samples within the expected uncertainty.

Accurate interpolation for conditional integrals
Many astrophysical inquiries in gravitational-wave astronomy require evaluating conditional integrals across parameter space, which in turns require sampling additional posterior points constrained to a hyperplane.This is for instance the case when estimating the Equation of State (EOS) from BNS collisions, an important post-processing analysis that allows us to probe extreme conditions of matter (Abbott et al. 2018).This is possible because the compactness of the objects is imprinted in the gravitational waveform and can be measured by the tidal deformability parameters.

GP samples PE samples
Chirp mass M/M 30.95The EOS integral involves evaluating the marginal posterior distribution over the masses (M, η) and tidal parameters ( Λ, δ Λ), subject to constraints between those parameters as parametrised by the EOS.
There are instances where the marginal posterior for these parameters contain non-linear correlations, as is the case for the first BNS event GW170817 (Abbott et al. 2017a).We test our interpolation model over this 4dimensional surface.In Figure 3 we compare the marginal distributions sampled from our GP model to the original PE samples.We see that our GP is able to faithfully represent the marginalised posterior surface, in particular, we see that there is good agreement between the 90% credible intervals.When looking at the 2d contours see that the 50% and 90% levels agree very well and that the GP model is able to capture degenerate features and bi-modalities.Finally, our interpolation of the surface can be re-sampled efficiently and for this example, we obtained 750k samples in a few O(5mins) (depending on hardware) using an HMC sampler.Hence this method can be advantageous over traditional methods, where the interpolation is generally performed with a Gaussian KDE by transforming the symmet-

Propagating GP uncertainty
GPs provide a fully Bayesian estimation of the uncertainty over model predictions, as the full covariance matrix between posterior samples is computed.In each of the GW applications shown so far we have utilised the mean prediction of the GP function.This uncertainty measurement can be very important in many cases, however here we illustrate with a single example how one can extract the uncertainty from the modelling.Accurate localisation of a gravitational signal can be of fundamental importance for multimessenger astronomy (Grover et al. 2014;Abbott et al. 2017b) and for measurements of cosmological parameters with dark sirens (Soares-Santos et al. 2019).As the localisation accuracy decreases, the marginal posteriors for the sky location parameters can look degenerate and non-Gaussian.We build an interpolation of the sky location parameters, right ascension (ra) and declination (dec), of GW150914.This event was observed by only two detectors, so albeit its high SNR, its sky location presents a typical ring-like shape.
The uncertainty measure produced by the GP is a Gaussian distribution about any given point on the surface, when considering the entire surface the combination of these Gaussians can be interpreted as a range of plausible density surfaces for any given confidence level (e.g.2σ).The uncertainty on the 1D marginal distributions can then be calculated by calculating an upper and lower bound for each point in the surface and then marginalising these across one of the dimensions to obtain an uncertainty estimate about the mean 1D predicted posterior density.
In 2D and especially when considering sky localisation, we are also interested in the contours that enclose a given volume of probability density (usually the 50% and 90% contours), to plan optimal observation strategies in the search of electromagnetic counterparts.We propagate the uncertainty estimate produced by the GP (in the space of all realisations from the GP) to the physical parameter space to obtain the difference in the area enclosed by a given contour integral between an upper and lower bound on the density surface.To calculate the uncertainty for a given confidence level, e.g.90%, we evaluate the upper and lower contour surfaces that correspond to the minimum probability density for the mean prediction contour which encloses 90% of the volume.
Geometrically, we are building an uncertainty envelope around the mean GP prediction.The uncertainty on the contour is then given by the location in physical parameter space where the edges of that envelope intersect the plane defined by the mean prediction contour.
In the central panel of Fig 4 we show the samples used to construct the model as well as the 50% and 90%, contours of the GP interpolation in 2D with their respective 2σ uncertainty (the shaded regions).The top and left panels of Fig. 4 show the mean prediction and its 2σ uncertainty marginalised over each parameter by a simple integration of the density over its projection.
The inclusion of the uncertainty highlights several features.On the central inset in Fig 4 we see that the lower bound on the 50% contour is composed of three islands which correspond to peaks, while for both the mean and the upper bound these islands are connected to obtain a smooth surface at this contour level.For the outer 90% contour we see that the differences mainly manifest in the tails, where as expected the upper bound follows the well known ring around the sky slightly further.This matches our intuition that there is possibly more density around the ring than around the edges of the contour in the middle of the plot.

CONCLUSIONS
We have presented an alternative method for density estimation of marginal PDFs for gravitational-wave parameters.Our method combines the desirable features of histograms to the extrapolation capabilities of KDEs, within a Bayesian framework.The choice of histogram binning determines the resolution of the PDF, while the kernel of the GP allows the interpolation to be flexible over non-Gaussian correlations and yet smooth.The noise variance parameter of the GP ensures that sources of stochastic noise from the histogram density estimation are taken into account.In cases where we employ an exact inference scheme, this noise variance can be evaluated for each histogram bin and it is equivalent to heteroskedastic errors over the density estimation.This allows to fully propagate the uncertainty from the PE samples.We plan to extend this method and fully incorporate uncertainties, as we showed in this work for the sky localisation example, over higher-dimensional posterior surfaces in future work.
This method may be preferable to other methods such as KDEs, a closely related method which is sometimes adopted in the field, depending upon the use-case requirements.It comes with three main advantages: it is suitable for most interpolation problems commonly encountered for gravitational-wave marginal posteriors; it provides a Bayesian measure of uncertainty over the model predic- tions; it allows to quickly re-sample the interpolation using HMC and other samplers available in TensorFlow .We presented a series of examples where we know the accuracy of the interpolation is important, such as EOS calculations and sky localisation.As the number of events will increase in the next observing run (O4), we need reliable tools to post-process the large volume of results.
This work has highlighted the power of GPs to fit a gravitational-wave posterior surface, a natural extension of this work is to generate a surrogate for the entire likelihood surface, similar to what was done by the authors of (Vivanco et al. 2019) using a random forest regressor.Such use of GPs has been already investigated in the field of cosmology to model the Planck18 posterior distribution (McClintock & Rozo 2019).This work has laid the foundation for us to apply a similar methodology to the gravitational-wave problem in a future work which is currently in preparation (d'Emilio et al. 2021).This has applications such as active sampling, efficient jump proposals (Graff et al. 2012;Farr et al. 2020) and more general use of the GP variance to guide the sampling process.The surface learned by the GP can be evaluated directly for a given set of parameters, therefore, avoiding the need to compute expensive waveforms.An example where such likelihood surrogates could be exploited is fast resampling with new astrophysical priors.This could replace an often difficult re-weighting procedure, especially when a prior assumption limits the number of available samples in a region of interest (Mandel & Fragos 2020).

A TECHNICAL DETAILS OF THE GP MODEL
A.1 Data pre-processing Data pre-processing, often referred to as data-set standardisation, is a common practice within the realm of machine learning and it can have a very high impact on the accuracy of the model.Our posterior samples have a wide range of values, some having bounds [−1, 1] and some reaching O(10 3 ).We re-scale our posterior samples such that each parameter ranges between [0, 1] by using the following transformation: where θ d is the vector of transformed samples and the min and max are evaluated for each parameter (i.e. each dimension of the posterior samples vector).The approximate marginalised posterior is scaled according to the z-score, such that it has zero mean and unit variance: where p(θi|d) is the transformed marginalised posterior, µ p(θ i |d) and σ p(θ i |d) are respectively the mean and standard deviation of the marginalised posterior points.All pre-processing in this work is performed using Scikit-Learn (Pedregosa et al. 2011).

A.2 Kernel design
The kernel is defined as the prior covariance between any two function values.Our prior knowledge about the morphology of the posterior can be encoded via this covariance, as it determines the space of functions that the GP sample paths live in.The radial basis function (RBF) or squared exponential kernel is the most basic kernel and it's given as: where the Euclidian distance between (x, x ) is scaled by the length-scale parameter (measure of deviations between points) and the overall variance is denoted by σ 2 (average distance of the function away from its mean).Functions drawn from a GP with this kernel are infinitely differentiable.
For our application, a more complex kernel architecture that can capture the correlations between parameters is needed.We need smoothness over small scale features, such that we don't model random noise fluctuations of samples, and flexibility over the large scale characteristics of the posterior.For this purpose we employ a combination of RBF and Matern, which is a generalisation of the RBF kernel with an additional smoothness parameter ν.The smaller ν, the less smooth the approximated function is: We choose ν = ( 1 2 , 5 2 ) depending on the specific morphology of the posterior, as this kernel is responsible for encoding its overall shape such as sharp boundary features.The resulting kernel equation is given by: The kernel multiplication is equivalent to an AND operation, as it corresponds to an element-wise multiplication of their corresponding covariance matrices.This means that the resulting covariance matrix will only have a high value if both covariances have a high value.We also apply automatic relevance determination (ARD), which modifies the kernel such that for each dimension an appropriate length scale is chosen (Neal 2012).

Figure 1 .
Figure 1.Interpolation of a bounded one-dimensional inverse gamma density function (in solid black) with our GP-based method (in solid orange).The histogram points used to generate the model and its uncertainty are shown as black points with error bars.Alternative KDE methods are shown for comparison as coloured dashed lines.

Figure 2 .
Figure 2. Corner plot of the intrinsic parameters of GW150914, drawn from our GP surrogate (in orange) compared to the original PE samples (in black).

Figure 3 .
Figure 3. Corner plot of the mass and tidal parameters of GW170817, drawn from our GP model (in orange), compared to original PE samples in black.
ric mass ratio parameter to be log(0.25 − η) (Pang et al. 2020) and there is no measure of uncertainty over the fit.

Figure 4 .
Figure 4. Central panel: contours of the 2D sky-location of GW150914, the GP model mean prediction and uncertainty (in orange) is compared to the points used to construct the fit (black crosses).Top and left panels show the GP model projections in 1D, compared to the original PE samples.All plots show the 2σ uncertainty around the density estimate as a shaded band

Table 1 .
Source properties of the intrinsic parameters of GW150914, original samples and samples from the GP interpolation.