Time-energy fit¶

3ML allows the possibility to model a time-varying source by explicitly fitting the time-dependent part of the model. Let’s see this with an example.

First we import what we need:

:
from threeML import *

import matplotlib.pyplot as plt

%matplotlib notebook
Configuration read from /home/giacomov/.threeML/threeML_config.yml
Plotter is MatPlotlib

Then we generate a simulated dataset for a source with a cutoff powerlaw spectrum with a constant photon index and cutoff but with a normalization that changes with time following a powerlaw:

:
# Let's generate our dataset of 4 spectra with a normalization that follows
# a powerlaw in time

def generate_one(K):

# Let's generate some data with y = Powerlaw(x)

gen_function = Cutoff_powerlaw()
gen_function.K = K

# Generate a dataset using the power law, and a
# constant 30% error

x = np.logspace(0, 2, 50)

xyl_generator = XYLike.from_function("sim_data", function = gen_function,
x = x,
yerr = 0.3 * gen_function(x))

y = xyl_generator.y
y_err = xyl_generator.yerr

return x, y, y_err

# These are the times at which the simulated spectra have been observed
time_tags = np.array([1.0, 2.0, 5.0, 10.0])

# This describes the time-varying normalization. If everything works as
# it should, we should recover from the fit a normalization of 0.23 and a
# index of -1.2 for the time law
normalizations = 0.23 * time_tags**(-1.2)

# Generate the datasets
datasets = map(generate_one, normalizations)
Using Gaussian statistic (equivalent to chi^2) with the provided errors.
Using Gaussian statistic (equivalent to chi^2) with the provided errors.
Using Gaussian statistic (equivalent to chi^2) with the provided errors.
Using Gaussian statistic (equivalent to chi^2) with the provided errors.

Now that we have our data, let’s model them with 3ML:

:
# Now set up the fit and fit it

# First we need to tell 3ML that we are going to fit using an
# independent variable (time in this case). We init it to 1.0
# and set the unit to seconds
time = IndependentVariable("time", 1.0, u.s)

# Then we load the data that we have generated, tagging them
# with their time of observation
plugins = []

for i, dataset in enumerate(datasets):

x, y, y_err = dataset

xyl = XYLike("data%i" % i, x, y, y_err)

# This is the important part: we need to tag the instance of the
# plugin so that 3ML will know that this instance corresponds to the
# given tag (a time coordinate in this case). If instead of giving
# one time coordinate we give two time coordinates, then 3ML will
# take the average of the model between the two time coordinates
# (computed as the integral of the model between t1 and t2 divided
# by t2-t1)

xyl.tag = (time, time_tags[i])

# To access the tag we have just set we can use:

independent_variable, start, end = xyl.tag

# NOTE: xyl.tag will return 3 things: the independent variable, the start and the
# end. If like in this case you do not specify an end when assigning the tag, end
# will be None

plugins.append(xyl)

# Generate the datalist as usual

data = DataList(*plugins)

# Now let's generate the spectral model, in this case a point source
# with a cutoff powerlaw spectrum

spectrum = Cutoff_powerlaw()

src = PointSource("test", ra=0.0, dec=0.0, spectral_shape=spectrum)

model = Model(src)

# Now we need to tell 3ML that we are going to use the time
# coordinate to specify a time dependence for some of the
# parameters of the model

# Now let's specify the time-dependence (a powerlaw) for the normalization
# of the powerlaw spectrum

time_po = Powerlaw()
time_po.K.bounds = (0.01, 1000)

# Link the normalization of the cutoff powerlaw spectrum with time through the
# time law we have just generated

# Now let's fit as usual

jl = JointLikelihood(model, data)

best_fit_parameters, likelihood_values = jl.fit()
Using Gaussian statistic (equivalent to chi^2) with the provided errors.
Using Gaussian statistic (equivalent to chi^2) with the provided errors.
Using Gaussian statistic (equivalent to chi^2) with the provided errors.
Using Gaussian statistic (equivalent to chi^2) with the provided errors.
Best fit values:

result unit
parameter
test.spectrum.main.Cutoff_powerlaw.K.Powerlaw.K (2.260 +/- 0.13) x 10^-1 1 / (cm2 keV s)
test...index -1.180 +/- 0.024
test.spectrum.main.Cutoff_powerlaw.index -1.995 +/- 0.031
test.spectrum.main.Cutoff_powerlaw.xc (1.007 +/- 0.016) x 10 keV

Correlation matrix:

 1 -0.51 -0.73 0.49 -0.51 1 0.03 -0.02 -0.73 0.03 1 -0.88 0.49 -0.02 -0.88 1

Values of -log(likelihood) at the minimum:

-log(likelihood)
data0 22.888946
data1 25.173887
data2 29.160302
data3 26.919687
total 104.142822

Values of statistical measures:

statistical measures
AIC 216.490772
BIC 229.478914
[ ]: