Skip to content
Snippets Groups Projects
Commit 84050660 authored by Nikhil Sarin's avatar Nikhil Sarin :mountain_snow: Committed by Gregory Ashton
Browse files

Added tutorial for fitting a line with x and y

parent 5c494930
No related branches found
No related tags found
No related merge requests found
%% Cell type:markdown id: tags:
# 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.
%% Cell type:code id: tags:
``` python
import bilby
import inspect
%pylab inline
```
%% Cell type:markdown id: tags:
### First we create the data and plot it
%% Cell type:code id: tags:
``` python
#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):
np.random.seed(int(seed))
xtrue = np.linspace(0,100,points)
ytrue = model(x = xtrue, m = m, c = c)
xerr = xerr * np.random.randn(points)
yerr = yerr * np.random.randn(points)
xobs = xtrue + xerr
yobs = ytrue + yerr
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()
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)
```
%% Cell type:markdown id: tags:
### Now lets set up the prior and bilby output directory
%% Cell type:code id: tags:
``` python
#setting up bilby priors
priors = dict(m=bilby.core.prior.Uniform(0, 30, 'm'),
c=bilby.core.prior.Uniform(0, 30, 'c'))
outdir = 'outdir'
livepoints = 100
walks = 100
```
%% Cell type:markdown id: tags:
### 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
%% Cell type:code id: tags:
``` python
OneDGaussian_knownx = bilby.core.likelihood.GaussianLikelihood(x = data['xtrue'], y = data['yobs'], func = model, sigma = data['yerr'])
result_1D_xtrue = bilby.run_sampler(
likelihood=OneDGaussian_knownx, priors=priors, sampler='dynesty', npoints=livepoints,
walks=walks, outdir=outdir, label='xtrue_1D_Gaussian')
```
%% Cell type:code id: tags:
``` python
result_1D_xtrue.plot_corner(truth=dict(m=5, c = 10), titles = True)
result_1D_xtrue.plot_with_data(model = model, x = data['xtrue'], y = data['yobs'], ndraws=1000, npoints=100)
plt.show()
```
%% Cell type:markdown id: tags:
### 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
%% Cell type:code id: tags:
``` python
OneDGaussian_unknownx = bilby.core.likelihood.GaussianLikelihood(x = data['xobs'], y = data['yobs'],
func = model, sigma = data['yerr'])
result_1D_xobs = bilby.run_sampler(
likelihood=OneDGaussian_unknownx, priors=priors, sampler='dynesty', npoints=livepoints,
walks=walks, outdir=outdir, label='xobs_1D_Gaussian')
```
%% Cell type:code id: tags:
``` python
result_1D_xobs.plot_corner(truth=dict(m=5, c = 10), titles = True)
result_1D_xobs.plot_with_data(model = model, x = data['xobs'], y = data['yobs'], ndraws=1000, npoints=100)
plt.show()
```
%% Cell type:markdown id: tags:
### As expected, this is significantly worse. Let us now define a new likelihood which takes into account x errors but you also pass in xtrue
%% Cell type:code id: tags:
``` python
class TwoDGaussianLikelihood_knownxtrue(bilby.Likelihood):
def __init__(self, xtrue, xobs, yobs, xerr, yerr, function):
"""
Parameters
----------
xtrue: array_like
The true injected x values
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
"""
self.xobs = xobs
self.xtrue = xtrue
self.yobs = yobs
self.yerr = yerr
self.xerr = xerr
self.function = function
parameters = inspect.getargspec(function).args
parameters.pop(0)
self.parameters = dict.fromkeys(parameters)
def log_likelihood(self):
resy = self.yobs - self.function(self.xtrue, **self.parameters)
resx = self.xobs - self.xtrue
return -0.5 * (np.sum(((resy) / self.yerr) ** 2) + np.sum(((resx) / self.xerr) ** 2))
```
%% Cell type:code id: tags:
``` python
TwoDGaussian_knownx = TwoDGaussianLikelihood_knownxtrue(xtrue = data['xtrue'], xobs = data['xobs'],
yobs = data['yobs'], xerr=data['xerr'],
yerr = data['yerr'], function=model)
result_2D_knownx = bilby.run_sampler(
likelihood=TwoDGaussian_knownx, priors=priors, sampler='dynesty', npoints=livepoints,
walks=walks, outdir=outdir, label='knownx_2D_Gaussian')
```
%% Cell type:code id: tags:
``` python
result_2D_knownx.plot_corner(truth=dict(m=5, c = 10), titles = True)
result_2D_knownx.plot_with_data(model = model, x = data['xobs'], y = data['yobs'], ndraws=1000, npoints=100)
plt.show()
```
%% Cell type:markdown id: tags:
### This works well, however it still is not realistic as one still needs to 'know' the true 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
%% Cell type:code id: tags:
``` python
class TwoDGaussianLikelihood_unknownx(bilby.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
"""
self.xobs = xobs
self.yobs = yobs
self.yerr = yerr
self.xerr = xerr
self.function = function
parameters = inspect.getargspec(function).args
parameters.pop(0)
self.parameters = dict.fromkeys(parameters)
def log_likelihood(self):
m = self.parameters['m']
v = np.array([-m, 1.0])
Sigma2 = (self.xerr*m)**2 + self.yerr**2
model_y = self.function(self.xobs, **self.parameters)
Delta = self.yobs - model_y
ll = -0.5 * np.sum(Delta**2 / Sigma2 + np.log(Sigma2))
return ll
```
%% Cell type:code id: tags:
``` python
TwoDGaussian_unknownx = TwoDGaussianLikelihood_unknownx(xobs = data['xobs'], yobs = data['yobs'],
xerr= data['xerr'], yerr = data['yerr'],
function=model)
result_2D_unknownx = bilby.run_sampler(
likelihood=TwoDGaussian_unknownx, priors=priors, sampler='dynesty', npoints=livepoints,
walks=walks, outdir=outdir, label='unknownx_2D_Gaussian')
```
%% Cell type:code id: tags:
``` python
result_2D_unknownx.plot_corner(truth=dict(m=5, c = 10), titles = True)
result_2D_unknownx.plot_with_data(model = model, x = data['xobs'], y = data['yobs'], ndraws=1000, npoints=100)
plt.show()
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment