Skip to content
Snippets Groups Projects
Commit f5a4adaf authored by Colm Talbot's avatar Colm Talbot
Browse files

update linear regression example

parent 3ed50894
No related branches found
No related tags found
No related merge requests found
#!/bin/python
#!/usr/bin/env python
"""
An example of how to use tupak to perform paramater estimation for
non-gravitational wave data. In this case, fitting a linear function to
......@@ -9,7 +9,6 @@ from __future__ import division
import tupak
import numpy as np
import matplotlib.pyplot as plt
import inspect
# A few simple setup steps
label = 'linear_regression'
......@@ -26,14 +25,14 @@ def model(time, m, c):
injection_parameters = dict(m=0.5, c=0.2)
# For this example, we'll use standard Gaussian noise
sigma = 1
# These lines of code generate the fake data. Note the ** just unpacks the
# contents of the injection_paramsters when calling the model function.
# contents of the injection_parameters when calling the model function.
sampling_frequency = 10
time_duration = 10
time = np.arange(0, time_duration, 1/sampling_frequency)
time = np.arange(0, time_duration, 1 / sampling_frequency)
N = len(time)
sigma = np.random.normal(1, 0.01, N)
data = model(time, **injection_parameters) + np.random.normal(0, sigma, N)
# We quickly plot the data to check it looks sensible
......@@ -45,59 +44,19 @@ ax.set_ylabel('y')
ax.legend()
fig.savefig('{}/{}_data.png'.format(outdir, label))
# Parameter estimation: we now define a Gaussian Likelihood class relevant for
# our model.
class GaussianLikelihoodKnownNoise(tupak.Likelihood):
def __init__(self, x, y, sigma, function):
"""
A general Gaussian likelihood - the parameters are inferred from the
arguments of function
Parameters
----------
x, y: array_like
The data to analyse
sigma: float
The standard deviation of the noise
function:
The python function to fit to the data. Note, this must take the
dependent variable as its first argument. The other arguments are
will require a prior and will be sampled over (unless a fixed
value is given).
"""
self.x = x
self.y = y
self.sigma = sigma
self.N = len(x)
self.function = function
# These lines of code infer the parameters from the provided function
parameters = inspect.getargspec(function).args
parameters.pop(0)
self.parameters = dict.fromkeys(parameters)
def log_likelihood(self):
res = self.y - self.function(self.x, **self.parameters)
return -0.5 * (np.sum((res / self.sigma)**2)
+ self.N*np.log(2*np.pi*self.sigma**2))
# Now lets instantiate a version of our GaussianLikelihood, giving it
# the time, data and signal model
likelihood = GaussianLikelihoodKnownNoise(time, data, sigma, model)
likelihood = tupak.likelihood.GaussianLikelihood(time, data, model, sigma)
# From hereon, the syntax is exactly equivalent to other tupak examples
# We make a prior
priors = {}
priors = dict()
priors['m'] = tupak.core.prior.Uniform(0, 5, 'm')
priors['c'] = tupak.core.prior.Uniform(-2, 2, 'c')
# And run sampler
result = tupak.run_sampler(
likelihood=likelihood, priors=priors, sampler='dynesty', npoints=500,
walks=10, injection_parameters=injection_parameters, outdir=outdir,
likelihood=likelihood, priors=priors, sampler='dynesty', nlive=500,
sample='unif', injection_parameters=injection_parameters, outdir=outdir,
label=label)
result.plot_corner()
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