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:

[1]:
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:

[2]:
# 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:

[3]:
# 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

model.add_independent_variable(time)

# 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

model.link(spectrum.K, time, time_po)

# 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.00-0.51-0.730.49
-0.511.000.03-0.02
-0.730.031.00-0.88
0.49-0.02-0.881.00

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
[ ]: