Bilby

Bayesian inference is the method by which we take physical observations and extract information about the Universe. It is a tool which is used ubiquitously in all areas of science and there are a plethora of methods for performing this inference. For a thorough introduction to Bayesian inference see this paper.

The most commonly used methods in gravitational-wave astronomy are stochastic sampling algorithms, e.g., Markov-chain Monte Carlo, and Nested Sampling. These are implemented in a wide variety of packages known as “samplers”, each suited to different problems. However, it is not clear when setting out on a project which is going to be best for the problem at hand, and these packages involve formulating the problem in a different way.

bilby is a unified framework for these packages, defining a standard for how to set up these problems. By doing this we allow users to seamlessly switch between samplers to validate robustness of our inferences.

An example reproduced from our git repository of how to perform a basic Bayesian inference problem in bilby is shown below. Many more examples for generic inference problems, compact binary coalescences, and population inference, can be found in our repository.

#!/usr/bin/env python
"""
An example of how to use bilby to perform paramater estimation for
non-gravitational wave data. In this case, fitting a linear function to
data with background Gaussian noise

"""
import numpy as np
import matplotlib.pyplot as plt

import bilby as bb

# First, we define our "signal model", in this case a simple linear function
# for a time series
def model(time, m, c):
    return time * m + c

# Now we define the injection parameters which we make simulated data with
injection_parameters = dict(m=0.5, c=0.2)

# These lines of code generate the fake data. Note the ** just unpacks the
# contents of the injection_parameters when calling the model function.
# We add some Gaussian noise to the data.
times = np.linspace(0, 10, 100)
sigma = 0.01
data = model(times, **injection_parameters) + np.random.normal(0, sigma, 100)

# We quickly plot the data to check it looks sensible
fig, ax = plt.subplots()
ax.plot(time, data, 'o', label='data')
ax.plot(time, model(time, **injection_parameters), '--r', label='signal')
ax.set_xlabel('time')
ax.set_ylabel('y')
ax.legend()
plt.show()

# Now create our likelihood, since the noise is Gaussian, we use a GaussianLikelihood
likelihood = bilby.likelihood.GaussianLikelihood(time, data, model, sigma)

# We must also define a prior, this is essentially a dictionary within bilby.
priors = bb.core.prior.PriorDict()
priors['m'] = bb.core.prior.Uniform(0, 5, 'm')
priors['c'] = bb.core.prior.Uniform(-2, 2, 'c')

# And run sampler
result = bb.run_sampler(
    likelihood=likelihood, priors=priors, sampler='dynesty', nlive=500,
    sample='unif', injection_parameters=injection_parameters,
    label='linear_regression')
result.plot_corner()