PyCBC Inference: A Python-based parameter estimation toolkit for compact binary coalescence signals
C. M. Biwer Collin D. Capano Soumi De Miriam Cabero Duncan A. Brown Alexander H. Nitz V. Raymond




Los Alamos National Laboratory, Los Alamos, NM 87545, USADepartment of Physics, Syracuse University, Syracuse, NY 13244, USA
Albert-Einstein-Institut, Max-Planck-Institut für Gravitationsphysik, D-30167 Hannover, Germany
Department of Physics, Syracuse University, Syracuse, NY 13244, USA
Albert-Einstein-Institut, Max-Planck-Institut für Gravitationsphysik, D-30167 Hannover, Germany
Department of Physics, Syracuse University, Syracuse, NY 13244, USA
Albert-Einstein-Institut, Max-Planck-Institut für Gravitationsphysik, D-30167 Hannover, Germany
Albert-Einstein-Institut, Max-Planck-Institut für Gravitationsphysik, D-14476 Potsdam, GermanySchool of Physics and Astronomy, Cardiff University, Cardiff, CF243AA, Wales, UK


We introduce new modules in the open-source PyCBC gravitational-wave astronomy toolkit that implement Bayesian inference for compact-object binary mergers. We review the Bayesian inference methods implemented and describe the structure of the modules. 
We demonstrate that the PyCBC Inference modules produce unbiased estimates of the parameters of a simulated population of binary black hole mergers. We show that the posterior parameter distributions obtained used our new code agree well with the published estimates for binary  black  holes in the  first LIGO-Virgo observing run.


PyCBC Inference: A Python-based parameter estimation toolkit for compact binary coalescence signals
    V. Raymond
    December 30, 2023
===================================================================================================





§ INTRODUCTION


The observations of six binary black hole
mergers <cit.>,
and the binary neutron star merger GW170817 <cit.>
by Advanced LIGO <cit.> and
Virgo <cit.> 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 <cit.>. 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 <cit.>.  PyCBC Inference has been used to
perform Bayesian inference for several 
astrophysical problems, including: testing the black hole area increase
law <cit.>; combining multi-messenger
obervations of GW170817 to constrain the viewing angle of the
binary <cit.>; 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 <cit.>. 

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 <cit.>.
This paper is organized as follows: Sec. <ref> 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. <ref> describes the design of the PyCBC Inference software and
how the methods described in  Sec. <ref> are implemented in the
code.  Sec. <ref> 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 <cit.>.  Finally, we summarize the
status of the code and possible future developments in Sec. <ref>.



§ BAYESIAN INFERENCE FOR BINARY MERGERS


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

    =
    / .
 
In our case, hypothesis H is the model of the gravitational-wave signal 
and  are the parameters of
this model. Together, these describe the properties of the astrophysical
source of the gravitational waves.
In Eq. (<ref>), the prior
probability density  describes our knowledge about the
parameters before considering the observed data d⃗(t), and the
likelihood  is the probability of obtaining the observation
d⃗(t) given the waveform model H with parameters .

Often we are only interested in a subset of the parameters .
To obtain a probability distribution on one or a few parameters, we marginalize
the posterior probability by integrating  over the unwanted
parameters. Marginalizing over all parameters yields the evidence, ,
which is the denominator in Eq. (<ref>). 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 <cit.>,

    ℬ = p(d⃗(t)|H_A)/p(d⃗(t)|H_B)  .

If ℬ is greater than 1 then model H_A is
favored over H_B, with the magnitude of ℬ 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
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.



 §.§ Waveform Models


The gravitational waves radiated in  a binary merger's source frame  are described by
the component masses m_1,2, the three-dimensional spin vectors
s⃗_1,2 of the compact objects <cit.>, and the
binary's eccentricity e <cit.>.  A parameter ϕ describes
the phase of the binary at a fiducial reference time, although this is not
usually of physical interest. For binaries containing neutron stars,
additional parameters Λ_1,2 describe the star's tidal
deformabilities <cit.>, which depend on the
equation of state of the neutron stars.  The waveform observed by the
Earth-based detector network depends on seven additional parameters: the
signal's time of arrival t_c, the binary's luminosity distance
d_L, and four Euler angles that describe the transformation from the
binary's frame to the detector network frame <cit.>.  These
angles are typically written as the binary's right ascension α,
declination δ, a polarization angle Ψ, and the inclination angle
ι (the angle between the binary's angular momentum axis and the line of
sight).  

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
ℳ <cit.>. 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 <cit.>. Similarly, the binary's mass ratio can be
degenerate with its spin <cit.>, 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. <cit.> 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 <cit.>, 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 , 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. <cit.>
and references therein), analytic models calibrated against numerical
simulations <cit.>,
perturbation theory <cit.>, and full numerical
solution of the Einstein equations (see e.g. Ref. <cit.> and
references therein). Obtaining posterior probabilities and evidences can
require calculating 𝒪(10^9) 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 <cit.>.  

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) <cit.>. 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 <cit.> 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
θ_j^a, and a polar angle θ_j^p for each of the two
compact objects. Examples of using PyCBC Inference with different waveform models include the analysis of
Ref. <cit.> that used the TaylorF2 post-Newtonian waveform model with
tidal corrections, and Ref. <cit.> that used a ringdown-only
waveform that models the quasi-normal modes of the remnant black hole.



 §.§ Likelihood Function


The data observed by the gravitational-wave detector network enters Bayes'
theorem through the likelihood  in Eq. (<ref>).
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 
that describes the location of the binary.

Under these assumptions, the appropriate form of  is the
well-known likelihood for a signal of known morphology in Gaussian noise (see
e.g.  Ref. <cit.> for its derivation), which is
given by

    = exp[ -1/2∑_i = 1^N<ñ_i(f) | ñ_i(f)> ] 
     
    = exp[ -1/2∑_i = 1^N<d̃_i(f) - s̃_i(f, )| d̃_i(f) - s̃_i(f, )> ] ,

where N is the number of detectors in the network.
The inner product ⟨ã | b̃⟩ is

    <ã_i(f) | b̃_i(f)> = 4∫_0^∞ã_i(f) b̃_i(f)/S^(i)_n(f)df  ,

where S^(i)_n(f) is the power spectral density of the of the i-th
detector's noise.  Here, d̃_i(f) and ñ_i(f) are the
frequency-domain representations of the data and noise, obtained by a Fourier
transformation of d_i(t) and n_i(t), respectively.  The model waveform
s̃_i(f, ) may be computed directly in the frequency
domain, or in the time domain and then Fourier transformed to the frequency
domain.  There are several operations (e.g. Fourier transforms, noise power
spectral density estimation, and inner products) that are common between the
calculation of Eq. (<ref>) and the computation of the matched
filter signal-to-noise ratio (SNR) in PyCBC <cit.>. PyCBC
Inference uses these existing functions, where appropriate.

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 ϕ,

    s̃_i(f, , ϕ) = s̃_i^0(f, , 0) e^iϕ.

The posterior probability  can be analytically marginalized over ϕ
for such models <cit.>. Assuming a uniform prior on
ϕ∈ [0,2π), the marginalized posterior is

    log   ∝log +
            I_0(|∑_i O(s̃^0_i, d̃_i)|) 
            - 1/2∑_i[ <s̃^0_i, s̃^0_i> -
                                    <d̃_i, d̃_i> ],

where

    s̃_i^0 ≡s̃_i(f, , ϕ=0),


    O(s̃^0_i, d̃_i) ≡ 4 ∫_0^∞s̃_i^*(f; ,0)d̃_i(f)/S^(i)_n(f)df,

and I_0 is the modified Bessel function of the first kind.

We have found that analytically
marginalizing over ϕ 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 <cit.>. 
Since fiducial phase is generally a
nuisance parameter, we use this form of the likelihood function in
Secs. <ref> and <ref>.




 §.§ Sampling Methods


Stochastic sampling techniques, and in particular Markov-chain Monte Carlo (MCMC) methods <cit.>, have been used to numerically sample the posterior probability density function of astrophysical parameters for binary-merger signals <cit.>.
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 _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 _l^(k).
The sampling algorithm chooses a new proposed set of parameters _l'^(k) with probability Q(_l^(k), _l'^(k)).
When a new set of parameters is proposed, the sampler computes an acceptance probability γ which determines if the Markov chain should move to the proposed parameter set _l'^(k) such that _l+1^(k) = _l'^(k).
If _l'^(i) is rejected, then _l+1^(k) = _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  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(_l^(k), _l'^(k)) and acceptance probability γ.

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  <cit.>, its parallel-tempered version  <cit.>, and the  <cit.> 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  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

    p_T(|d⃗(t),H) = ^1/T/p(d⃗(t)|H) .

The  sampler uses several temperatures in parallel, and the position of Markov chains are swapped between temperatures using an acceptance criteria described in Ref. <cit.>.
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 <cit.>.
The  sampler performs the sampling using one temperature where T = 1.

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



 §.§ Selection of Independent Samples


The output returned by the sampling algorithms discussed in Sec. <ref> are Markov chains. Successive states of these chains are not independent, as Markov processes depend on the previous state <cit.>.
The autocorrelation length τ_K of a Markov chain is a measure of the number of iterations required to produce independent samples of the posterior probability density function <cit.>.
The autocorrelation length of the k-th Markov chain X_l^(k) = {_g^(k) ; 1 < g < l } of length l obtained from the sampling algorithm is defined as 

    τ_K = 1 + 2 ∑_i=1^KR̂_i,

where K is the first iteration along the Markov chain the condition m τ_K ≤ K is true, m being a parameter which in PyCBC  Inference is set to 5 <cit.>.
The autocorrelation function R̂_i is defined as

    R̂_i = 1/l σ^2∑_t=1^l-i( X_t - μ) ( X_t+i - μ),

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, μ and σ^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,  and , which we have found to work well with the  and  samplers used in Sections <ref> and <ref>.

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

    logℒ_l^(k)≥max_k,llogℒ - N_p/2 ,

where ℒ is the prior-weighted likelihood

    ℒ =  ,

and N_p is the number of dimensions in the parameter space. The maximization max_k,llogℒ 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  sampler if the network signal-to-noise ratio of the signal is ≳ 5. 

While the  test works well with the  sampler, we have found that it underestimates the burn-in period when used with the  sampler.
Instead we use the  test with the  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 <cit.>. 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.






 §.§ 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. <ref>.
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 ∼ 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 <cit.>.
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  and the modules in the directory
.  PyCBC Inference provides an executable called
 that is the main engine for performing Bayesian
inference with PyCBC.
A call graph of  is shown in
Figure <ref>. In this section, we review the structure
of the main engine and the
Python objects used to build .




 §.§  executable


The methods presented in Sec. <ref>, <ref>, <ref>, <ref>, and <ref> are used to build the executable .
For faster performances,  can be run on high-throughput computing frameworks such as HTCondor <cit.> and the processes for running the sampler can be parallelized over multiple compute nodes using MPI <cit.>.
The execution of the likelihood computation and the PSD estimation are done using either single-threaded or parallel FFT engines, such as FFTW <cit.> 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 .

The input to  is a configuration file which contains up to seven types of sections.
The  section specifies the parameters that are to be varied in the MCMC.
There is a  section for each parameter in the  section 
which contains arguments to initialize the prior probability density function for that parameter.
There is a  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  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  section or the parameters inputted to the waveform generation functions.
Therefore, the configuration file may contain a  and  section(s) that specifies the transformations between parameters in the  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  section(s) may be provided which maps parameters in the  section to parameters understood by the waveform generation functions.
More details on application of constraints and execution of coordinate transformations are provided in Sec. <ref> and <ref> 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 .
The results from running  are stored in a

0pt 



HDF <cit.> file whose location is provided on the command line to  as well.
The main results of interest are stored under the HDF groups  and . The  group contains the history of the Markov chains as separate datasets for each of the variable parameters.
The  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(d⃗(t)|n⃗) and a dataset containing the natural logarithm of the prior probabilities.
The natural logarithm of the noise likelihood log p(d⃗(t)|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  group and the  group has shape  ×  if the sampling algorithm used in the analysis did not include parallel tempering, and has shape  ×  ×  for parallel-tempered samplers.
Here,  is the number of Markov chains,  is the number of iterations, and  is the number of temperatures.

 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  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  can resume the ensemble MCMC from the position of the Markov chains 
near the state the process was terminated.
The frequency  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  because the output file from  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  process can always be resumed.

For analyses that use the  sampler, the likelihood can be used to compute the natural logarithm of the evidence using the  sampler's  function <cit.>.
Then, the evidences from two analyses can be used to compute the Bayes factor ℬ 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. <cit.>.
These examples can be used with the open-source datasets provided by the LIGO Open Science Center <cit.>.
The results of these analyses are presented in Sec. <ref>.





 §.§  objects


The PyCBC Inference modules provide a set of  objects which execute the Bayesian sampling methods.
These objects provide classes and functions for using open-source samplers such as  <cit.>,  <cit.> or  <cit.>.
This acts as an interface between PyCBC Inference and the external sampler package.
The executable  initializes, executes, and saves the output from the  objects.
A particular  object is chosen on the command line of  with the  option. The  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. 



 §.§  objects
 

The  objects in PyCBC Inference are used to perform transformations between different coordinate systems.
Currently, the  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
ℳ 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 ℳ
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
(ℳ and q) proves to be more efficient than m_1 and
m_2 <cit.>.  In the GW150914,
LVT151012, and GW151226 configuration files in <cit.>, we
demonstrate how to allow the  object to provide priors in the
(m_1, m_2) coordinates, and specify sampling transformations to the
(ℳ, q) coordinates.

Waveform transforms are used to transform any variable parameters in the ensemble MCMC that may not be understood by the waveform model functions.
In PyCBC, the waveform model functions accept the following parameters: component masses m_1 and m_2, d_L, ι, t_c, ϕ, and any additional spin parameters in Cartesian coordinates. The convention adopted in PyCBC for ι denotes ι = 0 as a
“face-on” binary (line of sight parallel to binary angular momentum), ι
= π/2 as an “edge-on” binary (line of sight perpendicular to
binary angular momentum), and ι = π as a “face-off” binary (line of
sight anti-parallel to binary angular momentum).



 §.§  object


The  object computes the natural logarithm of the prior-weighted likelihood given by the numerator of Eq. <ref>.
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 γ instead of the full posterior probability density function.
The prior-weighted likelihood is computed for each new set of parameters as the  objects advance the Markov chains through the parameter space.



 §.§  objects


The  object must compute the prior probability density function p(|H). 
There exists several  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 s:


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

  *  : A cosine distribution.

  *  : A sine distribution.

  *  : A multivariate Gaussian distribution.

  *  : A multidimensional uniform distribution.

  *  : A uniform distribution between 0 and 2π.

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

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

  *  : A two-dimensional isotropic distribution.

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

Multiple  objects are needed to define the prior probability density function for all parameters.
The  object combines the individual prior probability density functions, providing a single interface for the  to evaluate the prior probability density function for all parameters.
As the sampling algorithm advances the positions of the Markov chains, the  computes the product of the prior probability density functions for the proposed new set of points in the parameter space.
The  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 ℳ and q obtained from the m_1 and m_2 prior probability density functions in the GW151226 configuration file in Ref. <cit.>.



 §.§  objects


As part of the likelihood calculation described in Sec. <ref>, a waveform s̃_i(f, ) is generated from a waveform model H and set of parameters .
PyCBC Inference provides  objects that allow waveforms s̃_i(f, ) to be generated for waveform models described in Sec. <ref> using PyCBC's interface to LAL <cit.>.
There are also  objects provided for generating ringdown waveforms as used in Ref. <cit.>.
Given the waveform model provided in the configuration file,  will automatically select the associated  object.



§ 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. <ref>.  This
figure was generated using the executable
 which make extensive use of tools
from the open-source packages Matplotlib <cit.> and
SciPy <cit.>. 



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. <cit.>.



 §.§ Simulated Signals



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 <cit.>.
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 <cit.>, 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 ∼10 and ∼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  and  samplers.
For the  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  burn-in test. For the  sampler, we use 5000 walkers and the  burn-in test.
As a result, we need only to run the  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 <cit.>.
This waveform model has 15 parameters. To reduce computational cost, we analytically marginalize over the fiducial phase ϕ by using Eq. (<ref>) 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 <ref> summarizes the result of this test using the  and  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  sampler have greater variance then the  sampler.

To quantify the deviations seen in Fig. <ref>, 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  and 0.03 for . In other words, if  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. <ref>. For the  sampler, the probability of obtaining a more extreme collection of curves than that seen in Fig. <ref> 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  with the above settings.
The  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.





 §.§ 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 <cit.>. We perform the parameter estimation analysis on the Advanced LIGO data available for these events at the LIGO Open Science Center <cit.>. We use the  sampler for these analyses.
For computing the likelihood, we analyze the  gravitational-wave dataset d⃗(t) from the Hanford and Livingston detectors.
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 <cit.> 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. <cit.> 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. <cit.>.
The dataset is sampled at 2048 Hz, and the likelihood is evaluated between a low frequency cutoff of 20 Hz and 1024 Hz.

The waveforms s̃_i(f, ) used in the likelihood are generated using the IMRPhenomPv2 <cit.> model implemented in the LIGO Algorithm Library (LAL) <cit.>.
The parameters inferred for these three events are ={α, δ, ψ, m_1, m_2, d_L, ι, t_c, a_1, a_2, θ_1^a, θ_2^a, θ_1^p, θ_2^p}, and we analytically marginalize over the fiducial phase ϕ.
These parameters form the complete set of parameters to construct a waveform from a binary black hole merger, and are the same parameters that were inferred from the parameter estimation analyses in Ref. <cit.>.  
Since faster convergence of m_1 and m_2 can be obtained with mass parameterizations of the waveform in ℳ and q we perform the coordinate transformation from (m_1, m_2) to (ℳ, q) before evaluating the priors.

We assume uniform prior distributions for the binary component masses m_1,2∈ [10, 80] M_⊙ for GW150914, m_1,2∈ [5, 80] M_⊙ for LVT151012, and m_1,2 corresponding to chirp mass ℳ∈ [9.5, 10.5] M_⊙ and mass ratio q ∈ [1, 18] for GW151226.
We use uniform priors on the spin magnitudes a_1,2∈ [0.0, 0.99].
We use a uniform solid angle prior, where θ_1,2^a is a uniform distribution θ_1,2^a∈ [0, 2π) and θ_1,2^p is a sine-angle distribution.
For the luminosity distance, we use a uniform in volume prior with d_L ∈ [10, 1000] Mpc for GW150914, d_L ∈ [10, 1500] Mpc for GW151226, and d_L ∈ [10, 2500] Mpc for LVT151012.
We use uniform priors for the arrival time t_c ∈ [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 <cit.>.
For the sky location parameters, we use a uniform distribution prior for α∈ [0, 2π) and a cosine-angle distribution prior for δ.
The priors described above are the same as those used in Ref. <cit.> . 

The parameter estimation analysis produces distributions that are a sampling of the posterior probability density function for the variable parameters from the ensemble MCMC.
We map these distributions obtained directly from the analysis to obtain estimates of other parameters of interest such as the chirp mass ℳ, mass ratio q, effective spin χ_eff, and the precession spin χ_p <cit.> parameters.
We use d_L to relate the detector-frame masses obtained from the MCMC to the source-frame masses using the standard Λ-CDM cosmology <cit.>. 

Recorded in Table <ref>, is a summary of the median and 90% credible interval values calculated for GW150914, GW151226, and LVT151012 analyses.
Results for m_1^src - m_2^src, q - χ_eff, and d_L - ι are shown in  Figs. <ref>, <ref>, and <ref> for GW150914, GW151226, and LVT151012 respectively.
The two-dimensional plots in these figures show the 50% and 90% credible regions, and the one-dimensional marginal distributions show the median and 90% credible intervals.
Overlaid are the one-dimensional marginal distributions, median, and 90% credible intervals, as well as the 50% and 90% credible regions using the samples obtained from the LIGO Open Science Center <cit.> for the analyses of the three events reported in Ref. <cit.> using the IMRPhenomPv2 model for comparison.
The results show that GW150914 has the highest mass components among the three events. GW150914 has more support for equal mass ratios whereas the GW151226 and LVT151012 posteriors support more asymmetric mass ratios. Overall, there is preference for smaller spins, with GW151226 having the highest spins among the three events. While the inclination and luminosity distances are not very well constrained, with generally a support for both “face-on” (ι = 0, line of sight parallel to binary angular momentum) and “face-off” (ι = π, line of sight anti-parallel to binary angular momentum) systems for all the three events, and GW150914 seem to have more preference for face-off systems.
We also computed χ_p for each of the three events and found no significant measurements of precession.
Overall, our results are in agreement with those presented in Ref. <cit.> within the statistical errors of measurement.

















§ 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. <cit.>. 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. <cit.> 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 <cit.>, (ii) measuring the viewing angle of GW170817 with electromagnetic and gravitational-wave signals <cit.>, and (iii) measuring the tidal deformabilities and radii of neutron stars from the observation of GW170817 <cit.>.
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





10
url<#>1#1urlprefixURL 



TheLIGOScientific:2016pea
Abbott B P et al. (Virgo, LIGO Scientific) 2016 Phys. Rev. X6 041015 (Preprint 1606.04856)

Abbott:2017vtc
Abbott B P et al. (VIRGO, LIGO Scientific) 2017 Phys. Rev.
  Lett. 118 221101 (Preprint 1706.01812)

Abbott:2017gyy
Abbott B P et al. (Virgo, LIGO Scientific) 2017 Astrophys. J.
  851 L35 (Preprint 1711.05578)

Abbott:2017oio
Abbott B P et al. (Virgo, LIGO Scientific) 2017 Phys. Rev.
  Lett. 119 141101 (Preprint 1709.09660)

TheLIGOScientific:2017qsa
Abbott B et al. (Virgo, LIGO Scientific) 2017 Phys. Rev. Lett.
  119 161101 (Preprint 1710.05832)

TheLIGOScientific:2014jea
Aasi J et al. (LIGO Scientific) 2015 Class. Quant. Grav. 32 074001 (Preprint 1411.4547)

TheVirgo:2014hva
Acernese F et al. (VIRGO) 2015 Class. Quant. Grav. 32
  024001 (Preprint 1408.3978)

Bayes:1763
Bayes M and Price M 1763 Philosophical Transactions (1683-1775)
  <http://dx.doi.org/10.1098/rstl.1763.0053>

Jaynes:2003jaq
Jaynes E T 2003 Probability Theory: The Logic of Science (CUP) ISBN
  9780521592710
  <http://www.cambridge.org/au/academic/subjects/physics/theoretical-physics-and-mathematical-physics/probability-theory-logic-science?format=HB   isbn=9780521592710>

alex_nitz_2018_1208115
Nitz A et al. 2018 Pycbc v1.9.4 <https://github.com/gwastro/pycbc>
  <https://doi.org/10.5281/zenodo.1208115>

PhysRevD.97.124069
Cabero M, Capano C D, Fischer-Birnholtz O, Krishnan B, Nielsen A B, Nitz A H
  and Biwer C M 2018 Phys. Rev. D 97(12) 124069
  <https://link.aps.org/doi/10.1103/PhysRevD.97.124069>

Finstad:2018wid
Finstad D, De S, Brown D A, Berger E and Biwer C M 2018 Astrophys. J.
  Lett. 860 L2 (Preprint 1804.04179)

De:2018uhw
De S, Finstad D, Lattimer J M, Brown D A, Berger E and Biwer C M 2018
  (Preprint 1804.08583)

ResultsRepo
<https://github.com/gwastro/pycbc-inference-paper>

Finn:1992xs
Finn L S and Chernoff D F 1993 Phys. Rev. D47 2198–2219
  (Preprint gr-qc/9301003)

Cutler:1994ys
Cutler C and Flanagan E E 1994 Phys. Rev. D49 2658–2697
  (Preprint gr-qc/9402014)

Nicholson:1997qh
Nicholson D and Vecchio A 1998 Phys. Rev. D57 4588–4599
  (Preprint gr-qc/9705064)

Christensen:2001cr
Christensen N and Meyer R 2001 Phys. Rev. D64 022001
  (Preprint gr-qc/0102018)

Allen:2005fk
Allen B, Anderson W G, Brady P R, Brown D A and Creighton J D E 2012 Phys.
  Rev. D85 122006

Usman:2015kfa
Usman S A et al. 2016 Class. Quant. Grav. 33 215004
  (Preprint 1508.02357)

Nitz:2017svb
Nitz A H, Dent T, Dal Canton T, Fairhurst S and Brown D A 2017 Astrophys.
  J. 849 118 (Preprint 1705.01513)

KassRaftery:1995
Kass R E and Raftery A E 1995 Journal of the American Statistical
  Association 90 773–795 (Preprint
  https://www.tandfonline.com/doi/pdf/10.1080/01621459.1995.10476572)
  <https://www.tandfonline.com/doi/abs/10.1080/01621459.1995.10476572>

GelmanMeng:1998
 1998 Statistical Science 13 163–185 ISSN 08834237
  <http://www.jstor.org/stable/2676756>

skilling2006
Skilling J 2006 Bayesian Anal. 1 833–859
  <https://doi.org/10.1214/06-BA127>

Hawking:1987en
Thorne K S 1987 THREE HUNDRED YEARS OF GRAVITATION ed Hawking S W and Israel
  W

Peters:1964zz
Peters P C 1964 Phys. Rev. 136 B1224–B1232

Flanagan:2007ix
Flanagan E E and Hinderer T 2008 Phys. Rev. D77 021502

Hinderer:2007mb
Hinderer T 2008 Astrophys. J. 677 1216–1220

Wahlquist:1987rx
Wahlquist H 1987 Gen. Rel. Grav. 19 1101–1113

PhysRev.131.435
Peters P C and Mathews J 1963 Phys. Rev. 131(1) 435–440
  <https://link.aps.org/doi/10.1103/PhysRev.131.435>

Hannam:2013uu
Hannam M, Brown D A, Fairhurst S, Fryer C L and Harry I W 2013 Astrophys.
  J. 766 L14 (Preprint 1301.5616)

Veitch:2014wba
Veitch J et al. 2015 Phys. Rev. D91 042003
  (Preprint 1409.7215)

Blanchet2006
Blanchet L 2006 Living Reviews in Relativity 9 4 ISSN 1433-8351
  <https://doi.org/10.12942/lrr-2006-4>

Buonanno:1998gg
Buonanno A and Damour T 1999 Phys. Rev. D59 084006
  (Preprint gr-qc/9811091)

Buonanno:2000ef
Buonanno A and Damour T 2000 Phys. Rev. D62 064015
  (Preprint gr-qc/0001013)

Damour:2000we
Damour T, Jaranowski P and Schaefer G 2000 Phys. Rev. D62 084011
  (Preprint gr-qc/0005034)

Damour:2001tu
Damour T 2001 Phys. Rev. D64 124013 (Preprint
  gr-qc/0103018)

Ajith:2007qp
Ajith P et al. 2007 Class. Quant. Grav. 24 S689–S700
  (Preprint 0704.3764)

Ajith:2009bn
Ajith P et al. 2011 Phys. Rev. Lett. 106 241101
  (Preprint 0909.2867)

Santamaria:2010yb
Santamaria L et al. 2010 Phys. Rev. D82 064016
  (Preprint 1005.3306)

Teukolsky:1972my
Teukolsky S A 1972 Phys. Rev. Lett. 29 1114–1118

Berti:2009kk
Berti E, Cardoso V and Starinets A O 2009 Class. Quant. Grav. 26
  163001 (Preprint 0905.2975)

Cardoso:2014uka
Cardoso V, Gualtieri L, Herdeiro C and Sperhake U 2015 Living Rev.
  Relativity 18 1 (Preprint 1409.0014)

Purrer:2015tud
Pürrer M 2016 Phys. Rev. D93 064041 (Preprint
  1512.02248)

Lackey:2016krb
Lackey B D, Bernuzzi S, Galley C R, Meidam J and Van Den Broeck C 2017 Phys. Rev. D95 104036 (Preprint 1610.04742)

lal
Mercer R A et al. 2017 LIGO Algorithm Library
  https://git.ligo.org/lscsoft/lalsuite

Schmidt:2014iyl
Schmidt P, Ohme F and Hannam M 2015 Phys. Rev. D91 024043
  (Preprint 1408.1810)

Hannam:2013oca
Hannam M, Schmidt P, Bohé A, Haegel L, Husa S, Ohme F, Pratten G and Pürrer M
  2014 Phys. Rev. Lett. 113 151101 (Preprint
  1308.3271)

Wainstein:1962
Wainstein L A and Zubakov V D 1962 Extraction of signals from noise
  (Englewood Cliffs, NJ: Prentice-Hall)

metropolis1953
Metropolis N, Rosenbluth A W, Rosenbluth M N, Teller A H and Teller E
  1953 J. Chem. Phys. 21 1087–1092

geman:1984
Geman S and Geman D 1984 IEEE Transactions on Pattern Analysis and Machine
  Intelligence PAMI-6 721–741 ISSN 0162-8828

Gilks:1996
Gilks W R, Richardson S and Spiegelhalter D J Markov Chain Monte Carlo in
  Practice, Chapman and Hall, London, 1996

Gelman95bayesiandata
Gelman A, Robert C, Chopin N and Rousseau J 1995 Bayesian data analysis

Christensen:2004bc
Christensen N, Libson A and Meyer R 2004 Class. Quant. Grav. 21
  317–330

Rover:2006ni
Rover C, Meyer R and Christensen N 2006 Class. Quant. Grav. 23
  4895–4906 (Preprint gr-qc/0602067)

Rover:2006bb
Rover C, Meyer R and Christensen N 2007 Phys. Rev. D75 062004
  (Preprint gr-qc/0609131)

emcee
Foreman-Mackey D, Hogg D W, Lang D and Goodman J 2013 Publ.
  Astron. Soc. Pac. 125 306 (Preprint 1202.3665)

emcee_repo
Foreman-Mackey D et al. 2018 dfm/emcee: emcee v3.0rc1
  <https://doi.org/10.5281/zenodo.1297477>

ptemcee_repo
Vousden W et al. 2015
  <https://github.com/willvousden/ptemcee>

vousden:2016
Vousden W D, Farr W M and Mandel I 2016 Monthly Notices of the Royal
  Astronomical Society 455 1919–1937 (Preprint
  /oup/backfile/content_public/journal/mnras/455/2/10.1093_mnras_stv2422/2/stv2422.pdf)
  <http://dx.doi.org/10.1093/mnras/stv2422>

kombine
Farr B and Farr W M 2015  In prep

kombine_repo
Farr B et al. 2018 <https://github.com/bfarr/kombine>

Hastings:1970aa
Hastings W K 1970 Biometrika 57 97–109

Christensen:2004jm
Christensen N, Dupuis R J, Woan G and Meyer R 2004 Phys. Rev. D70
  022001 (Preprint gr-qc/0402038)

Madras:1988ei
Madras N and Sokal A D 1988 J. Statist. Phys. 50 109–186

Chen2000
Chen M H, Shao Q M and Ibrahim J G 2000 Computing Bayesian Credible and
  HPD Intervals (New York, NY: Springer New York) ISBN 978-1-4612-1276-8
  <https://doi.org/10.1007/978-1-4612-1276-8_7>

beowulfbook-condor
Tannenbaum T, Wright D, Miller K and Livny M 2001 Condor – a distributed job
  scheduler Beowulf Cluster Computing with Linux ed Sterling T (MIT
  Press)

condor-practice
Thain D, Tannenbaum T and Livny M 2005 Concurrency - Practice and
  Experience 17 323–356

dalcin2011parallel
Dalcin L D, Paz R R, Kler P A and Cosimo A 2011 Advances in Water
  Resources 34 1124–1139

Dalcin:2008
Dalcín, Lisandro, Paz R, Storti M and D'Elia J 2008 Journal of
  parallel and distributed computing 68 655–662

Dalcin:2005:MP:1196218.1196228
Dalcín L, Paz R and Storti M 2005 J. Parallel Distrib. Comput.
  65 1108–1115 ISSN 0743-7315
  <http://dx.doi.org/10.1016/j.jpdc.2005.03.010>

FFTW05
Frigo M and Johnson S G 2005 Proceedings of the IEEE 93 216–231
  special issue on “Program Generation, Optimization, and Platform
  Adaptation”

andrew_collette_2018_1246321
Collette A et al. 2018 h5py/h5py 2.8.0
  <https://doi.org/10.5281/zenodo.1246321>

Vallisneri:2014vxa
Vallisneri M, Kanner J, Williams R, Weinstein A and Stephens B 2015 J.
  Phys. Conf. Ser. 610 012021 (Preprint 1410.4839)

Farr:2015lna
Farr B et al. 2016 Astrophys. J. 825 116
  (Preprint 1508.05336)

scipy_kde
<https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gaussian_kde.html>

Hunter:2007
Hunter J D 2007 Computing In Science & Engineering 9 90–95

scipy
Jones E, Oliphant T, Peterson P et al. 2001– SciPy: Open source
  scientific tools for Python [Online; accessed <today>]
  <http://www.scipy.org/>

TheLIGOScientific:2016wfe
Abbott B P et al. (Virgo, LIGO Scientific) 2016 Phys. Rev.
  Lett. 116 241102 (Preprint 1602.03840)

Abbott:2016blz
Abbott B P et al. (Virgo, LIGO Scientific) 2016 Phys. Rev.
  Lett. 116 061102 (Preprint 1602.03837)

Canton:2014ena
Dal Canton T et al. 2014 Phys. Rev. D90 082004
  (Preprint 1405.6731)

TheLIGOScientific:2016qqj
Abbott B P et al. (Virgo, LIGO Scientific) 2016 Phys. Rev. D93 122003 (Preprint 1602.03839)

Schutz:1986gp
Schutz B F 1986 Nature 323 310–311

Abbott:2016izl
Abbott T D et al. (Virgo, LIGO Scientific) 2016 Phys. Rev. X6 041014 (Preprint 1606.01210)



