Fitting a model to data with both x and y errors with Bilby

Usually when we fit a model to data with a Gaussian Likelihood we assume that we know x values exactly. This is almost never the case. Here we show how to fit a model with errors in both x and y.

import bilby
import matplotlib.pyplot as plt
import numpy as np
from bilby.core.utils import random

#sets seed of bilby's generator "rng" to "123"
random.seed(123)

%matplotlib inline

Simulate data

First we create the data and plot it

# define our model, a line
def model(x, m, c, **kwargs):
    y = m * x + c
    return y


# make a function to create and plot our data
def make_data(points, m, c, xerr, yerr, seed):
    xtrue = np.linspace(0, 100, points)
    ytrue = model(x=xtrue, m=m, c=c)

    xerr_vals = xerr * random.rng.standard_normal(points)
    yerr_vals = yerr * random.rng.standard_normal(points)
    xobs = xtrue + xerr_vals
    yobs = ytrue + yerr_vals

    plt.errorbar(xobs, yobs, xerr=xerr, yerr=yerr, fmt="x")
    plt.errorbar(xtrue, ytrue, yerr=yerr, color="black", alpha=0.5)
    plt.xlim(0, 100)
    plt.show()
    plt.close()

    data = {
        "xtrue": xtrue,
        "ytrue": ytrue,
        "xobs": xobs,
        "yobs": yobs,
        "xerr": xerr,
        "yerr": yerr,
    }

    return data


data = make_data(points=30, m=5, c=10, xerr=5, yerr=5, seed=123)

Define our prior and sampler settings

Now lets set up the prior and bilby output directory/sampler settings

# setting up bilby priors
priors = dict(
    m=bilby.core.prior.Uniform(0, 30, "m"), c=bilby.core.prior.Uniform(0, 30, "c")
)

sampler_kwargs = dict(priors=priors, sampler="bilby_mcmc", nsamples=1000, printdt=5, outdir="outdir", verbose=False, clean=True)

Fit with exactly known x-values

Our first step is to recover the straight line using a simple Gaussian Likelihood that only takes into account the y errors. Under the assumption we know x exactly. In this case, we pass in xtrue for x

known_x = bilby.core.likelihood.GaussianLikelihood(
    x=data["xtrue"], y=data["yobs"], func=model, sigma=data["yerr"]
)
result_known_x = bilby.run_sampler(
    likelihood=known_x,
    label="known_x",
    **sampler_kwargs,
)

/home/runner/work/bilby/bilby/bilby/core/utils/io.py:64: UserWarning: Wswiglal-redir-stdio:

SWIGLAL standard output/error redirection is enabled in IPython.
This may lead to performance penalties. To disable locally, use:

with lal.no_swig_redirect_standard_output_error():
    ...

To disable globally, use:

lal.swig_redirect_standard_output_error(False)

Note however that this will likely lead to error messages from
LAL functions being either misdirected or lost when called from
Jupyter notebooks.

To suppress this warning, use:

import warnings
warnings.filterwarnings("ignore", "Wswiglal-redir-stdio")
import lal

  import lal

_ = result_known_x.plot_corner(truth=dict(m=5, c=10), titles=True, save=False)
plt.show()
plt.close()

Fit with unmodeled uncertainty in the x-values

As expected this is easy to recover and the sampler does a good job. However this was made too easy - by passing in the ‘true’ values of x. Lets see what happens when we pass in the observed values of x

incorrect_x = bilby.core.likelihood.GaussianLikelihood(
    x=data["xobs"], y=data["yobs"], func=model, sigma=data["yerr"]
)
result_incorrect_x = bilby.run_sampler(
    likelihood=incorrect_x,
    label="incorrect_x",
    **sampler_kwargs,
)

_ = result_incorrect_x.plot_corner(truth=dict(m=5, c=10), titles=True, save=False)
plt.show()
plt.close()

Fit with modeled uncertainty in x-values

This is not good as there is unmodelled uncertainty in our x values. Getting around this requires marginalisation of the true x values or sampling over them. See discussion in section 7 of https://arxiv.org/pdf/1008.4686.pdf.

For this, we will have to define a new likelihood class. By subclassing the base bilby.core.likelihood.Likelihood class we can do this fairly simply.

class GaussianLikelihoodUncertainX(bilby.core.likelihood.Likelihood):
    def __init__(self, xobs, yobs, xerr, yerr, function):
        """

        Parameters
        ----------
        xobs, yobs: array_like
            The data to analyse
        xerr, yerr: array_like
            The standard deviation of the noise
        function:
            The python function to fit to the data
        """
        super(GaussianLikelihoodUncertainX, self).__init__()
        self.xobs = xobs
        self.yobs = yobs
        self.yerr = yerr
        self.xerr = xerr
        self.function = function

    def log_likelihood(self, parameters):
        variance = (self.xerr * parameters["m"]) ** 2 + self.yerr**2
        model_y = self.function(self.xobs, **parameters)
        residual = self.yobs - model_y

        return -0.5 * np.sum(residual**2 / variance + np.log(variance))

gaussian_unknown_x = GaussianLikelihoodUncertainX(
    xobs=data["xobs"],
    yobs=data["yobs"],
    xerr=data["xerr"],
    yerr=data["yerr"],
    function=model,
)
result_unknown_x = bilby.run_sampler(
    likelihood=gaussian_unknown_x,
    label="unknown_x",
    **sampler_kwargs,
)

_ = result_unknown_x.plot_corner(truth=dict(m=5, c=10), titles=True, save=False)
plt.show()
plt.close()

Success! The inferred posterior is consistent with the true values.