PyCBC Inference: A Python-based parameter estimation toolkit for compact binary coalescence signals

Introduction

The observations of six binary black hole mergers , and the binary neutron star merger GW170817 by Advanced LIGO and Virgo have established the field of gravitational-wave astronomy. Understanding the origin, evolution, and physics of gravitational-wave sources requires accurately measuring the properties of detected events. In practice, this is performed using Bayesian inference . Bayesian inference allows us to determine the signal model that is best supported by observations and to obtain posterior probability densities for a model’s parameters, hence inferring the properties of the source. In this paper, we present PyCBC Inference; a set of Python modules that implement Bayesian inference in the PyCBC open-source toolkit for gravitational-wave astronomy . PyCBC Inference has been used to perform Bayesian inference for several astrophysical problems, including: testing the black hole area increase law ; combining multi-messenger obervations of GW170817 to constrain the viewing angle of the binary ; and determining that the gravitational-wave observations of GW170817 favor a model where both compact objects have the same equation of state, and measuring the tidal deformabilities and radii of the neutron stars .

We provide a comprehensive description of the methods and code implemented in PyCBC Inference. We then demonstrate that PyCBC Inference can produce unbiased estimates of the parameters of a simulated population of binary black holes. We show that PyCBC Inference can recover posterior probability distributions that are in good agreement with the published measurements of the binary black holes detected in the first LIGO-Virgo observing run . This paper is organized as follows: Sec. 2 gives an overview of the Bayesian inference methods used in gravitational-wave astronomy for compact-object binary mergers. We provide an overview of the waveform models used; the likelihood function for a known signal in stationary, Gaussian noise; the sampling methods used to estimate the posterior probability densities and the evidence; the selection of independent samples; and the estimation of parameter values from posterior probabilities. Sec. 3 describes the design of the PyCBC Inference software and how the methods described in Sec. 2 are implemented in the code. Sec. 4 uses a simulated population of binary black holes and the black-hole mergers detected in the first LIGO-Virgo observing run to demonstrate the use of PyCBC Inference. We provide the posterior probability densities for the events GW150914, GW151226, and LVT151012, and the command lines and configurations to reproduce these results as supplemental materials . Finally, we summarize the status of the code and possible future developments in Sec. 5.

Bayesian Inference for Binary Mergers

In gravitational-wave astronomy, Bayesian methods are used to infer the properties of detected astrophysical sources . Given the observed data $\vec{d}(t)$—here this is data from a gravitational-wave detector network in which a search has identified a signal —Bayes’ theorem states that for a hypothesis $H$,

In our case, hypothesis $H$ is the model of the gravitational-wave signal and $\vec{\vartheta}$ are the parameters of this model. Together, these describe the properties of the astrophysical source of the gravitational waves. In Eq. (1), the prior probability density $p(\vec{\vartheta}|H)$ describes our knowledge about the parameters before considering the observed data $\vec{d}(t)$, and the likelihood $p(\vec{d}(t)|\vec{\vartheta},H)$ is the probability of obtaining the observation $\vec{d}(t)$ given the waveform model $H$ with parameters $\vec{\vartheta}$.

Often we are only interested in a subset of the parameters $\vec{\vartheta}$. To obtain a probability distribution on one or a few parameters, we marginalize the posterior probability by integrating $p(\vec{d}(t)|\vec{\vartheta},H)p(\vec{\vartheta}|H)$ over the unwanted parameters. Marginalizing over all parameters yields the evidence, $p(\vec{d}(t)|H)$, which is the denominator in Eq. (1). The evidence serves as a normalization constant of the posterior probability for the given model $H$. If we have two competing models $H_{A}$ and $H_{B}$, the evidence can be used to determine which model is favored by the data via the Bayes factor ,

If $\mathcal{B}$ is greater than 1 then model $H_{A}$ is favored over $H_{B}$, with the magnitude of $\mathcal{B}$ indicating the degree of belief.

PyCBC Inference can compute Bayes factors and produce marginalized posterior probability densities given the data from a network of gravitational-wave observatories with $N$ detectors $\vec{d}(t)=\{d_{i}(t);1<i<N\}$, and a model $H$ that describes the astrophysical source. In the remainder of this section, we review the methods used to compute these quantities.

Binary mergers present a challenging problem for Bayesian inference, as the dimensionality of the signal parameter space is large. This is further complicated by correlations between the signal’s parameters. For example, at leading order the gravitational waveform depends on the chirp mass $\mathcal{M}$ . The mass ratio enters the waveform at higher orders and is more difficult to measure. This results in an amplitude-dependent degeneracy between the component masses . Similarly, the binary’s mass ratio can be degenerate with its spin , although this degeneracy can be broken if the binary is precessing. Much of the effort of parameter estimation in gravitational-wave astronomy has focused on developing computationally feasible ways to explore this signal space, and on extracting physically interesting parameters (or combinations of parameters) from the large, degenerate parameter space (see e.g. Ref. and references therein). However, in many problems of interest, we are not concerned with the full parameter space described above. For example, field binaries are expected to have negligible eccentricity when they are observed by LIGO and Virgo , and so eccentricity is neglected in the waveform models. Simplifying assumptions can be made about the compact object’s spins (e.g. the spins are aligned with the binary’s orbital angular momentum), reducing the dimensionality of the waveform parameters space.

Given a set of parameters $\vec{\vartheta}$, one can obtain a model of the gravitational-wave signal from a binary merger using a variety of different methods, including: post-Newtonian theory (see e.g. Ref. and references therein), analytic models calibrated against numerical simulations , perturbation theory , and full numerical solution of the Einstein equations (see e.g. Ref. and references therein). Obtaining posterior probabilities and evidences can require calculating $\mathcal{O}\left(10^{9}\right)$ template waveforms, which restricts us to models that are computationally efficient to calculate. The cost of full numerical simulations makes them prohibitively expensive at present. Even some analytic models are too costly to be used, and surrogate models have been developed that capture the features of these waveforms at reduced computational cost .

The specific choice of the waveform model $H$ for an analysis depends on the physics that we wish to explore, computational cost limitations, and the level of accuracy desired in the model. A variety of waveform models are available for use in PyCBC Inference, either directly implemented in PyCBC or via calls to the LIGO Algorithm Library (LAL) . We refer to the PyCBC and LAL documentation, and references therein, for detailed descriptions of these models. In this paper, we demonstrate the use of PyCBC Inference using the IMRPhenomPv2 waveform model for binary black hole mergers. This model captures the inspiral-merger-ringdown physics of spinning, precessing binaries and parameterizing spin effects using a spin magnitude $a_{j}$, an azimuthal angle $\theta_{j}^{a}$, and a polar angle $\theta_{j}^{p}$ for each of the two compact objects. Examples of using PyCBC Inference with different waveform models include the analysis of Ref. that used the TaylorF2 post-Newtonian waveform model with tidal corrections, and Ref. that used a ringdown-only waveform that models the quasi-normal modes of the remnant black hole.

2 Likelihood Function

The data observed by the gravitational-wave detector network enters Bayes’ theorem through the likelihood $p(\vec{d}(t)|\vec{\vartheta},H)$ in Eq. (1). Currently, PyCBC Inference assumes that the each detector produces stationary, Gaussian noise $n_{i}(t)$ that is uncorrelated between the detectors in the network. The observed data is then $d_{i}(t)=n_{i}(t)+s_{i}(t)$, where $s_{i}(t)$ is the gravitational waveform observed in the $i$-th detector. For detectors that are not identical and co-located (as in the case of the LIGO-Virgo network), each detector observers a slightly different waveform due to their different antennae patterns, however the signal in the $i$-th detector can be calculated given the subset of the parameters $\vec{\vartheta}$ that describes the location of the binary.

Under these assumptions, the appropriate form of $p(\vec{d}(t)|\vec{\vartheta},H)$ is the well-known likelihood for a signal of known morphology in Gaussian noise (see e.g. Ref. for its derivation), which is given by

In general, gravitational-wave signals consist of a superposition of harmonic modes. However, in many cases it is sufficient to model only the most dominant mode, since the sub-dominant harmonics are too weak to be measured. In this case, the signal observed in all detectors has the same simple dependence on the fiducial phase $\phi$,

The posterior probability $p(\vec{\vartheta}|\vec{d}(t),H)$ can be analytically marginalized over $\phi$ for such models . Assuming a uniform prior on $\phi\in[0,2\pi)$, the marginalized posterior is

and $I_{0}$ is the modified Bessel function of the first kind.

We have found that analytically marginalizing over $\phi$ in this manner reduces the computational cost of the analysis by a factor of 2 – 3. The IMRPhenomPv2 model that we use here is a simplified model of precession that allows for this analytic marginalization . Since fiducial phase is generally a nuisance parameter, we use this form of the likelihood function in Secs. 4.1 and 4.2.

3 Sampling Methods

Stochastic sampling techniques, and in particular Markov-chain Monte Carlo (MCMC) methods , have been used to numerically sample the posterior probability density function of astrophysical parameters for binary-merger signals . Ensemble MCMC algorithms use multiple Markov chains to sample the parameter space. A simple choice to initialize the $k$-th Markov chain in the ensemble is to draw a set of parameters $\vec{\vartheta}_{1}^{(k)}$ from the prior probability density function. The Markov chains move around the parameter space according to the following set of rules. At iteration $l$, the $k$-th Markov chain has the set of parameters $\vec{\vartheta}_{l}^{(k)}$. The sampling algorithm chooses a new proposed set of parameters $\vec{\vartheta}_{l^{\prime}}^{(k)}$ with probability $Q(\vec{\vartheta}_{l}^{(k)},\vec{\vartheta}_{l^{\prime}}^{(k)})$. When a new set of parameters is proposed, the sampler computes an acceptance probability $\gamma$ which determines if the Markov chain should move to the proposed parameter set $\vec{\vartheta}_{l^{\prime}}^{(k)}$ such that $\vec{\vartheta}_{l+1}^{(k)}=\vec{\vartheta}_{l^{\prime}}^{(k)}$. If $\vec{\vartheta}_{l^{\prime}}^{(i)}$ is rejected, then $\vec{\vartheta}_{l+1}^{(k)}=\vec{\vartheta}_{l}^{(k)}$. After a sufficient number of iterations, the ensemble converges to a distribution that is proportional to a sampling of the posterior probability density function. The true astrophysical parameters $\vec{\vartheta}$ can then be estimated from histograms of the position of the Markov chains in the parameter space. Different ensemble sampling algorithms make particular choices for the proposal probability $Q(\vec{\vartheta}_{l}^{(k)},\vec{\vartheta}_{l^{\prime}}^{(k)})$ and acceptance probability $\gamma$.

The open-source community has several well-developed software packages that implement algorithms for sampling the posterior probability density function. PyCBC Inference leverages these developments, and we have designed a flexible framework that allows the user to choose from multiple ensemble sampling algorithms. Currently, PyCBC Inference supports the open-source ensemble sampler emcee , its parallel-tempered version emcee_pt , and the kombine sampler. All of the three are ensemble MCMC samplers. The sampling algorithm advances the positions of the walkers based on their previous positions and provides PyCBC Inference the positions of the walkers along the Markov chain.

The emcee_pt sampler is a parallel-tempered sampler which advances multiple ensembles based on the tempering or the “temperatures” used to explore the posterior probability density function. The posterior probability density function for a particular temperature $T$ is modified such that

The emcee_pt sampler uses several temperatures in parallel, and the position of Markov chains are swapped between temperatures using an acceptance criteria described in Ref. . Mixing of Markov chains from the different temperatures makes parallel-tempered samplers suitable for sampling posterior probability density functions with widely separated modes in the parameter space . The emcee sampler performs the sampling using one temperature where $T=1$.

The kombine sampler on the other hand uses clustered kernel-density estimates to construct its proposal distribution, and proposals are accepted using the Metropolis–Hastings condition . The kombine sampler has been included in PyCBC Inference due to its efficient sampling which significantly lowers the computational cost of an analysis relative to the emcee_pt sampler. However, in Sec. 4.1, we found that the nominal configuration of the kombine sampler produced biased estimates of parameters for binary black holes.

4 Selection of Independent Samples

The output returned by the sampling algorithms discussed in Sec. 2.3 are Markov chains. Successive states of these chains are not independent, as Markov processes depend on the previous state . The autocorrelation length $\tau_{K}$ of a Markov chain is a measure of the number of iterations required to produce independent samples of the posterior probability density function . The autocorrelation length of the $k$-th Markov chain $X_{l}^{(k)}=\{\vec{\vartheta}_{g}^{(k)};1<g<l\}$ of length $l$ obtained from the sampling algorithm is defined as

where $K$ is the first iteration along the Markov chain the condition $m\tau_{K}\leq K$ is true, $m$ being a parameter which in PyCBC Inference is set to $5$ . The autocorrelation function $\hat{R}_{i}$ is defined as

where $X_{t}$ are the samples of $X_{l}^{(k)}$ between the 0-th and the $t$-th iteration, $X_{t+i}$ are the samples of $X_{l}^{(k)}$ between the 0-th and the $(t+1)$-th iterations. Here, $\mu$ and $\sigma^{2}$ are the mean and variance of $X_{t}$, respectively.

The initial positions of the Markov chains influence their subsequent positions. The length of the Markov chains before they are considered to have lost any memory of the initial positions is called the “burn-in” period. It is a common practice in MCMC analyses to discard samples from the burn-in period to prevent any bias introduced by the initial positions of the Markov chains on the estimates of the parameters from the MCMC. PyCBC Inference has several methods to determine when the Markov chains are past the burn-in period. Here, we describe two methods, max_posterior and n_acl, which we have found to work well with the kombine and emcee_pt samplers used in Sections 4.1 and 4.2.

The max_posterior algorithm is an implementation of the burn-in test used for the MCMC sampler in Ref. . In this method, the $k$-th Markov chain is considered to be past the burn-in period at the first iteration $l$ for which

where $\mathcal{L}$ is the prior-weighted likelihood

and $N_{p}$ is the number of dimensions in the parameter space. The maximization $\max_{k,l}\log\mathcal{L}$ is carried out over all Markov chains and iterations. The ensemble is considered to be past the burn-in period at the first iteration where all chains pass this test. We have found this test works well with the kombine sampler if the network signal-to-noise ratio of the signal is $\gtrsim 5$.

While the max_posterior test works well with the kombine sampler, we have found that it underestimates the burn-in period when used with the emcee_pt sampler. Instead we use the n_acl test with the emcee_pt sampler. This test posits that the sampler is past the burn-in period if the length of the chains exceed $10$ times the autocorrelation length. The autocorrelation length is calculated using samples from the second half of the Markov chains. If the test is satisfied, the sampler is considered to be past the burn-in period at the midway point of the Markov chains.

Correlations between the neighboring samples after the burn-in period are removed by “thinning” or drawing samples from the Markov chains with an interval of the autocorrelation length . This is done so that the samples used to estimate the posterior probability density function are independent. Therefore, the number of independent samples of the posterior probability density function is equal to the number of Markov chains used in the ensemble times the number of iterations after the burn-in period divided by the autocorrelation length. PyCBC Inference will run until it has obtained the desired number of independent samples after the burn-in period.

5 Credible Intervals

After discarding samples from the burn-in period and thinning the remaining samples of the Markov chains, the product is the set of independent samples as described in Sec. 2.4. Typically we summarize the measurement of a given parameter using a credible interval. The $x$% credible interval is an interval where the true parameter value lies with a probability of $x$%. PyCBC Inference provides the capability to calculate credible intervals based on percentile values. In the percentile method, the $x$% credible interval of a parameter value is written as $A_{-B}^{+C}$ where $A$ is typically the 50-th percentile (median) of the marginalized histograms. The values $A-B$ and $A+C$ represent the lower and upper boundaries of the $x$-th percentile respectively.

An alternative method of calculating a credible interval estimate is the Highest Posterior Density (HPD) method. An $x\%$ HPD interval is the shortest interval that contains $x\%$ of the probability. The percentile method explained above imposes a non-zero lower boundary to the interval being measured. This can be perceived as a limitation in cases where the weight of histogram at the $\sim$ 0-th percentile is not significantly different from the weight at the lower boundary of the credible interval. Intervals constructed using the HPD method may be preferred in such cases. Previous studies have noted that HPD intervals may be useful when the posterior distribution is not symmetric . PyCBC Inference uses HPD to construct confidence contours for two-dimensional marginal distributions, but HPD is not used in the contruction of one-dimensional credible intervals for a single parameter. This functionality will be included in a future release of PyCBC Inference.

The PyCBC Inference Toolkit

In this section we describe the implementation of PyCBC Inference within the broader PyCBC toolkit. PyCBC provides both modules for developing code and executables for performing specific tasks with these modules. The code is available on the public GitHub repository at https://github.com/gwastro/pycbc, with executables located in the directory bin/inference and the modules in the directory pycbc/inference. PyCBC Inference provides an executable called pycbc_inference that is the main engine for performing Bayesian inference with PyCBC. A call graph of pycbc_inference is shown in Figure 1. In this section, we review the structure of the main engine and the Python objects used to build pycbc_inference.

The methods presented in Sec. 2.1, 2.2, 2.3, 2.4, and 2.5 are used to build the executable pycbc_inference. For faster performances, pycbc_inference can be run on high-throughput computing frameworks such as HTCondor and the processes for running the sampler can be parallelized over multiple compute nodes using MPI . The execution of the likelihood computation and the PSD estimation are done using either single-threaded or parallel FFT engines, such as FFTW or the Intel Math Kernel Library (MKL). For maximum flexibility in heterogeneous computing environments, the processing scheme to be used is specified at runtime as a command line option to pycbc_inference.

The input to pycbc_inference is a configuration file which contains up to seven types of sections. The variable_args section specifies the parameters that are to be varied in the MCMC. There is a prior section for each parameter in the variable_args section which contains arguments to initialize the prior probability density function for that parameter. There is a static_args section specifying any parameter for waveform generation along with its assigned value that should be fixed in the ensemble MCMC. Optionally, the configuration file may also include a constraint section(s) containing any conditions that constrain the prior probability density functions of the parameters. For efficient convergence of a Markov chain, it may be desirable to sample the prior probability density function in a different coordinate system than the parameters defined in the variable_args section or the parameters inputted to the waveform generation functions. Therefore, the configuration file may contain a sampling_parameters and sampling_transform section(s) that specifies the transformations between parameters in the variable_args sections and the parameters evaluated in the prior probability density function. Finally, the waveform generation functions recognize only a specific set of input parameters. The waveform_transforms section(s) may be provided which maps parameters in the variable_args section to parameters understood by the waveform generation functions. More details on application of constraints and execution of coordinate transformations are provided in Sec. 3.3 and 3.4 respectively.

The location of the configuration file, gravitational-wave detector data files, data conditioning settings, and settings for the ensemble MCMC are supplied on the command line interface to pycbc_inference. The results from running pycbc_inference are stored in a

HDF file whose location is provided on the command line to pycbc_inference as well. The main results of interest are stored under the HDF groups [‘samples’] and [‘likelihood_stats’]. The [‘samples’] group contains the history of the Markov chains as separate datasets for each of the variable parameters. The [‘likelihood_stats’] group contains a dataset of the natural logarithm of the Jacobian which is needed to transform from the variable parameters to sampling parameters, a dataset containing natural logarithm of the likelihood ratio $\log p(\vec{d}(t)|\vec{\vartheta},H)/p(\vec{d}(t)|\vec{n})$ and a dataset containing the natural logarithm of the prior probabilities. The natural logarithm of the noise likelihood $\log p(\vec{d}(t)|\vec{n})$ is stored as an attribute in the output file, and the likelihood is the summation of this quantity with the natural logarithm of the likelihood ratio. Each of the datasets under the [‘samples’] group and the [‘likelihood_stats’] group has shape nwalkers $\times$ niterations if the sampling algorithm used in the analysis did not include parallel tempering, and has shape ntemps $\times$ nwalkers $\times$ niterations for parallel-tempered samplers. Here, nwalkers is the number of Markov chains, niterations is the number of iterations, and ntemps is the number of temperatures.

pycbc_inference has checkpointing implemented which allows users to resume an analysis from the last set of Markov chains positions written to the output file. It is computationally expensive to obtain the desired number of independent samples using ensemble MCMC methods, and the pycbc_inference processes may terminate early due to problems on distributed-computing networks. Therefore, the samples from the Markov chains should be written at regular intervals so pycbc_inference can resume the ensemble MCMC from the position of the Markov chains near the state the process was terminated. The frequency pycbc_inference writes the samples from Markov chains and the state of the random number generator to the output file and a backup file is specified by the user on the command line. A backup file is written by pycbc_inference because the output file from pycbc_inference may be corrupted. For example, if the process is aborted while writing to the output file, then the output file may be corrupted. In that case, samples and the state of the random number generator are loaded from the backup file, and the backup file is copied to the output file. This ensures that the pycbc_inference process can always be resumed.

For analyses that use the emcee_pt sampler, the likelihood can be used to compute the natural logarithm of the evidence using the emcee_pt sampler’s thermodynamic_integration_log_evidence function . Then, the evidences from two analyses can be used to compute the Bayes factor $\mathcal{B}$ for the comparison of two waveform models.

We provide example configuration files and run scripts for the analysis of the binary black hole mergers detected in the Advanced LIGO’s first observing run in Ref. . These examples can be used with the open-source datasets provided by the LIGO Open Science Center . The results of these analyses are presented in Sec. 4.2.

2 Sampler objects

The PyCBC Inference modules provide a set of Sampler objects which execute the Bayesian sampling methods. These objects provide classes and functions for using open-source samplers such as emcee , emcee_pt or kombine . This acts as an interface between PyCBC Inference and the external sampler package. The executable pycbc_inference initializes, executes, and saves the output from the Sampler objects. A particular Sampler object is chosen on the command line of pycbc_inference with the --sampler option. The Sampler object provides the external sampler package the positions of the walkers in the parameter space, the natural logarithm of the posterior probabilities at the current iteration, the current “state” determined from a random number generator, and the number of iterations that the sampler is requested to run starting from the current iteration. After running for the given number of iterations, the sampler returns the updated positions of the Markov chains, the natural logarithm of the posterior probabilities, and the new state.

3 Transform objects

The Transform objects in PyCBC Inference are used to perform transformations between different coordinate systems. Currently, the Transform objects are used in two cases: sampling transforms and waveform transforms.

Sampling transforms are used for transforming parameters that are varied in the ensemble MCMC to a different coordinate system before evaluating the prior probaility density function. Since there exists degeneracies between several parameters in a waveform model it is useful to parameterize the waveform using a preferred set of parameters which could minimize the correlations. This leads to more efficient sampling, and therefore, it leads to faster convergence of the Markov chains. One example of a sampling transformation is the transformation between the component masses $m_{1}$ and $m_{2}$ to chirp mass $\mathcal{M}$ and mass ratio $q$. The convention adopted for $q$ in PyCBC Inference is $q=m_{1}/m_{2}$, where $m_{1}$ and $m_{2}$ are the component masses with $m_{1}>m_{2}$. The chirp mass $\mathcal{M}$ is the most accurately measured parameter in a waveform model because it is in the leading order term of the post-Newtonian expression of the waveform model. In contrast, the degeneracies of the mass ratio with spin introduces uncertainties in measurements of the component masses. Therefore, sampling in $(\mathcal{M}$ and $q)$ proves to be more efficient than $m_{1}$ and $m_{2}$ . In the GW150914, LVT151012, and GW151226 configuration files in , we demonstrate how to allow the Sampler object to provide priors in the $(m_{1},m_{2})$ coordinates, and specify sampling transformations to the $(\mathcal{M},q)$ coordinates.

4 LikelihoodEvaluator object

The LikelihoodEvaluator object computes the natural logarithm of the prior-weighted likelihood given by the numerator of Eq. 1. Since the evidence is constant for a given waveform model, then the prior-weighted likelihood is proportional to the posterior probability density function, and it can be used in sampling algorithms to compute the acceptance probability $\gamma$ instead of the full posterior probability density function. The prior-weighted likelihood is computed for each new set of parameters as the Sampler objects advance the Markov chains through the parameter space.

5 Distribution objects

The LikelihoodEvaluator object must compute the prior probability density function $p(\vec{\vartheta}|H)$. There exists several Distribution objects that provide functions for evaluating the prior probability density function to use for each parameter, and for drawing random samples from these distributions. Currently, PyCBC Inference provides the following Distributions:

Arbitrary : Reads a set of samples stored in a HDF format file and uses Gaussian kernel-density estimation to construct the distribution.

Gaussian : A multivariate Gaussian distribution.

Uniform : A multidimensional uniform distribution.

UniformAngle : A uniform distribution between 0 and 2$\pi$.

UniformLog : A multidimensional distribution that is uniform in its logarithm.

UniformPowerLaw : A multidimensional distribution that is uniform in a power law.

UniformSky : A two-dimensional isotropic distribution.

UniformSolidAngle : A two-dimensional distribution that is uniform in solid angle.

Multiple Distribution objects are needed to define the prior probability density function for all parameters. The JointDistribution object combines the individual prior probability density functions, providing a single interface for the LikelihoodEvaluator to evaluate the prior probability density function for all parameters. As the sampling algorithm advances the positions of the Markov chains, the JointDistribution computes the product of the prior probability density functions for the proposed new set of points in the parameter space. The JointDistribution can apply additional constraints on the prior probability density functions of parameters and it renormalizes the prior probability density function accordingly. If multiple constraints are provided, then the union of all constraints are applied. We demonstrate how to apply a cut on $\mathcal{M}$ and $q$ obtained from the $m_{1}$ and $m_{2}$ prior probability density functions in the GW151226 configuration file in Ref. .

6 Generator objects

Validation of the Toolkit

PyCBC Inference includes tools for visualizing the results of parameter estimation, and several analytic functions that can be used to test the generation of known posterior probabilities. Two common ways to visualize results are a scatter plot matrix of the independent samples of the Markov chains, and a marginalized one-dimensional histograms showing the bounds of each parameter’s credible interval. Analytic likelihood functions available to validate the code include: the multivariate normal, Rosenbrock, eggbox, and volcano functions. An example showing the visualization of results from the multi-variate Gaussian test is shown in Fig. 2. This figure was generated using the executable pycbc_inference_plot_posterior which make extensive use of tools from the open-source packages Matplotlib and SciPy .

We can also validate the performance of PyCBC Inference by: (i) determining if the inferred parameters of a population of simulated signals agrees with known the parameters that population, and (ii) comparing PyCBC Inference’s parameter credible intervals astrophysical signals to the published LIGO-Virgo results that used a different inference code. In this section, we first check that the credible intervals match the probability of finding the simulated signal parameters in that interval, that is, that $x\%$ of signals should have parameter values in the $x\%$ credible interval. We then compare the recovered parameters of the binary black hole mergers GW150914, GW151226, and LVT151012 to those published in Ref. .

To test the performance of PyCBC Inference, we generate 100 realizations of stationary Gaussian noise colored by power-spectral densities representative of the sensitivity of Advanced LIGO detectors at the time of the detection of GW150914 . To each realization of noise we add a simulated signal whose parameters were drawn from the same prior probability density function used in the analysis of GW150914 , with an additional cut placed on distance to avoid having too many injections with low matched-filter SNR. The resulting injections have matched-filter SNRs between $5$ and $160$, with the majority between $\sim 10$ and $\sim 40$. We then perform a parameter estimation analysis on each signal to obtain credible intervals on all parameters.

We perform this test using both the emcee_pt and kombine samplers. For the emcee_pt sampler we use 200 walkers and 20 temperatures. We run the sampler until we obtain at least 2000 independent samples after the burn-in period as determined using the n_acl burn-in test. For the kombine sampler, we use 5000 walkers and the max_posterior burn-in test. As a result, we need only to run the kombine sampler until the burn-in test is satisfied, at which point we immediately have 5000 independent samples of the posterior probability density function.

Both simulated signals and the waveforms in the likelihood computation are generated using IMRPhenomPv2 . This waveform model has 15 parameters. To reduce computational cost, we analytically marginalize over the fiducial phase $\phi$ by using Eq. (6) for the posterior probability, thereby reducing the number of sampled parameters to 14. For each parameter, we count the number of times the simulated parameter falls within the measured credible interval.

Figure 3 summarizes the result of this test using the emcee_pt and kombine samplers. For each of the parameters we plot the fraction of signals whose true parameter value fall within a credible interval as a function of credible interval (this is referred to as a percentile-percentile plot). We expect the former to equal the latter for all parameters, though some fluctuation is expected due to noise. We see that all parameters follow a 1-to-1 relation, though the results from the kombine sampler have greater variance then the emcee_pt sampler.

To quantify the deviations seen in Fig. 3, we perform a Kolmogorov–Smirnov (KS) test on each parameter to see whether the percentile-percentile curves match the expected 1-to-1 relation. If the samplers and code are performing as expected, then these p-values should in turn follow a uniform distribution. We therefore perform another KS test on the collection of p-values, obtaining a two-tailed p-value of $0.50$ for emcee_pt and $0.03$ for kombine. In other words, if emcee_pt provides an unbiased estimate of the parameters, then there is a $50\%$ chance that we would obtain a collection of percentile-percentile curves more extreme than seen in Fig. 3. For the kombine sampler, the probability of obtaining a more extreme collection of curves than that seen in Fig. 3 is only $3\%$.

Based on these results, we conclude that PyCBC Inference does indeed provide unbiased estimates of binary black hole parameters when used with emcee_pt with the above settings. The kombine sampler does not appear to provide unbiased parameter estimates when used to sample the full parameter space of precessing binary black holes with the settings we have used.

2 Astrophysical Events

In this section, we present PyCBC Inference measurements of properties of the binary black hole sources of the two gravitational-wave signals GW150914 and GW151226, and the third gravitational-wave signal LVT151012 consistent with the properties of a binary black hole source from Advanced LIGO’s first observing run . We perform the parameter estimation analysis on the Advanced LIGO data available for these events at the LIGO Open Science Center . We use the emcee_pt sampler for these analyses. For computing the likelihood, we analyze the gravitational-wave dataset $\vec{d}(t)$ from the Hanford and Livingston detectors. $\vec{d}(t)$ in our analyses are taken from GPS time intervals 1126259452 to 1126259468 for GW150914, 1135136340 to 1135136356 for GW151226, and 1128678874 to 1128678906 for LVT151012. Detection of gravitational waves from the search pipeline gives initial estimates of the mass, and hence estimates of the length of the signal. From results of the search, LVT151012 was a longer signal with more cycles than the other two events, and LVT151012 had characteristics which were in agreement with a lower mass source than GW150914 and GW151226. Therefore, more data is required for the analysis of LVT151012. The PSD used in the likelihood is constructed using the median PSD estimation method described in Ref. with 8 s Hann-windowed segments ( overlapped by 4 s ) taken from GPS times 1126258940 to 1126259980 for GW150914, 1135136238 to 1135137278 for GW151226, and 1128678362 to 1128679418 for LVT151012. The PSD estimate is truncated to 4 s in the time-domain using the method described in Ref. . The dataset is sampled at 2048 Hz, and the likelihood is evaluated between a low frequency cutoff of 20 Hz and 1024 Hz.

We assume uniform prior distributions for the binary component masses $m_{1,2}\in$ M⊙ for GW150914, $m_{1,2}\in$ M⊙ for LVT151012, and $m_{1,2}$ corresponding to chirp mass $\mathcal{M}\in$ [9.5, 10.5] M⊙ and mass ratio $q\in$ for GW151226. We use uniform priors on the spin magnitudes $a_{1,2}\in$ [0.0, 0.99]. We use a uniform solid angle prior, where $\theta_{1,2}^{a}$ is a uniform distribution $\theta_{1,2}^{a}\in[0,2\pi)$ and $\theta_{1,2}^{p}$ is a sine-angle distribution. For the luminosity distance, we use a uniform in volume prior with $d_{L}\in$ Mpc for GW150914, $d_{L}\in$ Mpc for GW151226, and $d_{L}\in$ Mpc for LVT151012. We use uniform priors for the arrival time $t_{c}\in[t_{s}-0.2~{}s,t_{s}+0.2~{}s]$ where $t_{s}$ is the trigger time for the particular event obtained from the gravitational-wave search . For the sky location parameters, we use a uniform distribution prior for $\alpha\in[0,2\pi)$ and a cosine-angle distribution prior for $\delta$. The priors described above are the same as those used in Ref. .

Conclusions

In this paper we have described PyCBC Inference, a Python-based toolkit with a simplified interface for parameter estimation studies of compact-object binary mergers. We have used this toolkit to estimate the parameters of the gravitational-wave events GW150914, GW151226, and LVT151012; our results are consistent with previously published values. In these analyses, we do not marginalize over calibration uncertainty of the measured strain in our results, which was included in prior work, for example Refs. . We will implement this in PyCBC Inference in the future. We have made the samples of the posterior probability density function from the PyCBC Inference analysis of all three events available in Ref. along with the instructions and configuration files needed to replicate these results. The source code and documentation for PyCBC Inference is available as part of the PyCBC software package at http://pycbc.org.

PyCBC Inference has already been used to produce several astrophysical results: (i) a test of the black hole area increase law , (ii) measuring the viewing angle of GW170817 with electromagnetic and gravitational-wave signals , and (iii) measuring the tidal deformabilities and radii of neutron stars from the observation of GW170817 . The results presented in this paper and in the studies above demonstrate the capability of PyCBC Inference to perform gravitational-wave parameter estimation analyses. Future developments under consideration are implementation of models to marginalize over calibration errors, generic algorithms to perform model selection, HPD to compute credible intervals, and methods for faster computation of the likelihood.

Acknowledgements

The authors would like to thank Will Farr and Ben Farr for valuable insights into the intricacies of ensemble MCMCs. We also thank Ian Harry, Christopher Berry, and Daniel Wysocki for helpful comments on the manuscript. This work was supported by NSF awards PHY-1404395 (DAB, CMB), PHY-1707954 (DAB, SD), and PHY-1607169 (SD). Computations were supported by Syracuse University and NSF award OAC-1541396. We also acknowledge the Max Planck Gesellschaft for support and the Atlas cluster computing team at AEI Hannover. DAB thanks the École de Physique des Houches for hospitality during the completion of this manuscript. The authors thank the LIGO Scientific Collaboration for access to the data and acknowledge the support of the United States National Science Foundation (NSF) for the construction and operation of the LIGO Laboratory and Advanced LIGO as well as the Science and Technology Facilities Council (STFC) of the United Kingdom, and the Max-Planck-Society (MPS) for support of the construction of Advanced LIGO. Additional support for Advanced LIGO was provided by the Australian Research Council. This research has made use of data obtained from the LIGO Open Science Center https://losc.ligo.org.

References

References