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 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:
|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|
Values of -log(likelihood) at the minimum:
Values of statistical measures: