alternative_likelihoods.py 3.24 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
#!/bin/python
"""
An example of how to use tupak to perform paramater estimation for
non-gravitational wave data
"""
from __future__ import division
import tupak
import numpy as np
import matplotlib.pyplot as plt

# A few simple setup steps
tupak.utils.setup_logger(log_level="info")
label = 'test'
outdir = 'outdir'


# Here is minimum requirement for a Likelihood class needed to run tupak. In
# this case, we setup a GaussianLikelihood, which needs to have a
# log_likelihood and noise_log_likelihood method. Note, in this case we make
# use of the `tupak` waveform_generator to make the signal (more on this later)
# But, one could make this work without the waveform generator.

class GaussianLikelihood():
    def __init__(self, x, y, waveform_generator):
        self.x = x
        self.y = y
        self.N = len(x)
        self.waveform_generator = waveform_generator

    def log_likelihood(self):
        sigma = 1
        res = self.y - self.waveform_generator.time_domain_strain()
        return -0.5 * (np.sum((res / sigma)**2)
                       + self.N*np.log(2*np.pi*sigma**2))

    def noise_log_likelihood(self):
        sigma = 1
        return -0.5 * (np.sum((self.y / sigma)**2)
                       + self.N*np.log(2*np.pi*sigma**2))


# Here we define our signal model, in this case a very basic trig. function
def model(time, A, P):
    return A*np.sin(2*np.pi*time/P)


# Here we define the injection parameters which we make simulated data with
injection_parameters = dict(A=1.5, P=10)

# 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.
sampling_frequency = 10
time_duration = 100
time = np.arange(0, time_duration, 1/sampling_frequency)
N = len(time)
data = model(time, **injection_parameters) + np.random.normal(0, sigma, N)

# We quickly plot the data to check it looks sensible
fig, ax = plt.subplots()
ax.plot(time, data)
ax.plot(time, model(time, **injection_parameters), '--r')
fig.savefig('{}/data.png'.format(outdir))

# Here, we make a `tupak` waveform_generator. In this case, of course, the
# name doesn't make so much sense. But essentially this is an objects that
# can generate a signal. We give it information on how to make the time series
# and the model() we wrote earlier.
71
72
73
waveform_generator = tupak.waveform_generator.WaveformGenerator(time_duration=time_duration,
                                                                sampling_frequency=sampling_frequency,
                                                                time_domain_source_model=model)
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92


# Now lets instantiate a version of out Likelihood, giving it the time, data
# and waveform_generator
likelihood = GaussianLikelihood(time, data, waveform_generator)

# From hereon, the syntax is exactly equivalent to other tupak examples
# We make a prior
priors = {}
priors['A'] = tupak.prior.Uniform(0, 5, 'A')
priors['P'] = tupak.prior.Uniform(0, 20, 'P')

# And run sampler
result = tupak.sampler.run_sampler(
    likelihood=likelihood, priors=priors, sampler='dynesty', npoints=1000,
    walks=10, injection_parameters=injection_parameters, outdir=outdir,
    label=label, use_ratio=False)
result.plot_corner()
print(result)