Time-energy fit

In [1]:

    

from threeML import *

import matplotlib.pyplot as plt

%matplotlib notebook


    

    




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

In [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.

In [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.73
0.49
-0.51
1.00
0.03
-0.02
-0.73
0.03
1.00
-0.88
0.49
-0.02
-0.88
1.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
    
  

In [ ]: