Analyzing GRB 080916C

Alt text (NASA/Swift/Cruz deWilde)

To demonstrate the capabilities and features of 3ML in, we will go through a time-integrated and time-resolved analysis. This example serves as a standard way to analyze Fermi-GBM data with 3ML as well as a template for how you can design your instrument’s analysis pipeline with 3ML if you have similar data.

3ML provides utilities to reduce time series data to plugins in a correct and statistically justified way (e.g., background fitting of Poisson data is done with a Poisson likelihood). The approach is generic and can be extended. For more details, see the time series documentation.

[1]:

import warnings

warnings.simplefilter("ignore")

[2]:

%%capture
import matplotlib.pyplot as plt
import numpy as np

np.seterr(all="ignore")


from threeML import *
from threeML.io.package_data import get_path_of_data_file

[3]:

silence_warnings()
%matplotlib inline
from jupyterthemes import jtplot

jtplot.style(context="talk", fscale=1, ticks=True, grid=False)
set_threeML_style()

Examining the catalog

As with Swift and Fermi-LAT, 3ML provides a simple interface to the on-line Fermi-GBM catalog. Let’s get the information for GRB 080916C.

[4]:

gbm_catalog = FermiGBMBurstCatalog()
gbm_catalog.query_sources("GRB080916009")

00:18:17 INFO      The cache for fermigbrst does not yet exist. We will try to    get_heasarc_table_as_pandas.py:64
                  build it

         INFO      Building cache for fermigbrst                                 get_heasarc_table_as_pandas.py:112

[4]:

Table length=1

name	ra	dec	trigger_time	t90
object	float64	float64	float64	float64
GRB080916009	119.800	-56.600	54725.0088613	62.977

To aid in quickly replicating the catalog analysis, and thanks to the tireless efforts of the Fermi-GBM team, we have added the ability to extract the analysis parameters from the catalog as well as build an astromodels model with the best fit parameters baked in. Using this information, one can quickly run through the catalog an replicate the entire analysis with a script. Let’s give it a try.

[5]:

grb_info = gbm_catalog.get_detector_information()["GRB080916009"]

gbm_detectors = grb_info["detectors"]
source_interval = grb_info["source"]["fluence"]
background_interval = grb_info["background"]["full"]
best_fit_model = grb_info["best fit model"]["fluence"]
model = gbm_catalog.get_model(best_fit_model, "fluence")["GRB080916009"]

[6]:

model

[6]:

Model summary:

	N
Point sources	1
Extended sources	0
Particle sources	0

Free parameters (5):

	value	min_value	max_value	unit
GRB080916009.spectrum.main.SmoothlyBrokenPowerLaw.K	0.012255	0.0	None	keV-1 s-1 cm-2
GRB080916009.spectrum.main.SmoothlyBrokenPowerLaw.alpha	-1.130424	-1.5	2.0
GRB080916009.spectrum.main.SmoothlyBrokenPowerLaw.break_energy	309.2031	10.0	None	keV
GRB080916009.spectrum.main.SmoothlyBrokenPowerLaw.break_scale	0.3	0.0	10.0
GRB080916009.spectrum.main.SmoothlyBrokenPowerLaw.beta	-2.096931	-5.0	-1.6

Fixed parameters (3):
(abridged. Use complete=True to see all fixed parameters)

Properties (0):

(none)

Linked parameters (0):

(none)

Independent variables:

(none)

Linked functions (0):

(none)

Downloading the data

We provide a simple interface to download the Fermi-GBM data. Using the information from the catalog that we have extracted, we can download just the data from the detectors that were used for the catalog analysis. This will download the CSPEC, TTE and instrument response files from the on-line database.

[7]:

dload = download_GBM_trigger_data("bn080916009", detectors=gbm_detectors)

Let’s first examine the catalog fluence fit. Using the TimeSeriesBuilder, we can fit the background, set the source interval, and create a 3ML plugin for the analysis. We will loop through the detectors, set their appropriate channel selections, and ensure there are enough counts in each bin to make the PGStat profile likelihood valid.

First we use the CSPEC data to fit the background using the background selections. We use CSPEC because it has a longer duration for fitting the background.
The background is saved to an HDF5 file that stores the polynomial coefficients and selections which we can read in to the TTE file later.
The light curve is plotted.
The source selection from the catalog is set and DispersionSpectrumLike plugin is created.
The plugin has the standard GBM channel selections for spectral analysis set.

[8]:

fluence_plugins = []
time_series = {}
for det in gbm_detectors:
    ts_cspec = TimeSeriesBuilder.from_gbm_cspec_or_ctime(
        det, cspec_or_ctime_file=dload[det]["cspec"], rsp_file=dload[det]["rsp"]
    )

    ts_cspec.set_background_interval(*background_interval.split(","))
    ts_cspec.save_background(f"{det}_bkg.h5", overwrite=True)

    ts_tte = TimeSeriesBuilder.from_gbm_tte(
        det,
        tte_file=dload[det]["tte"],
        rsp_file=dload[det]["rsp"],
        restore_background=f"{det}_bkg.h5",
    )

    time_series[det] = ts_tte

    ts_tte.set_active_time_interval(source_interval)

    ts_tte.view_lightcurve(-40, 100)

    fluence_plugin = ts_tte.to_spectrumlike()

    if det.startswith("b"):
        fluence_plugin.set_active_measurements("250-30000")

    else:
        fluence_plugin.set_active_measurements("9-900")

    fluence_plugin.rebin_on_background(1.0)

    fluence_plugins.append(fluence_plugin)

00:19:18 INFO      Auto-determined polynomial order: 0                                binned_spectrum_series.py:389

00:19:30 INFO      None 0-order polynomial fit with the mle method                               time_series.py:458

         INFO      Saved fitted background to n3_bkg.h5                                         time_series.py:1064

         INFO      Saved background to n3_bkg.h5                                         time_series_builder.py:471

         INFO      Successfully restored fit from n3_bkg.h5                              time_series_builder.py:171

         INFO      Interval set to 1.28-64.257 for n3                                    time_series_builder.py:290

         INFO      Auto-probed noise models:                                                    SpectrumLike.py:490

         INFO      - observation: poisson                                                       SpectrumLike.py:491

         INFO      - background: gaussian                                                       SpectrumLike.py:492

         INFO      Range 9-900 translates to channels 5-124                                    SpectrumLike.py:1247

00:19:34 INFO      Now using 120 bins                                                          SpectrumLike.py:1739

00:19:36 INFO      Auto-determined polynomial order: 1                                binned_spectrum_series.py:389

00:19:52 INFO      None 1-order polynomial fit with the mle method                               time_series.py:458

         INFO      Saved fitted background to n4_bkg.h5                                         time_series.py:1064

         INFO      Saved background to n4_bkg.h5                                         time_series_builder.py:471

00:19:53 INFO      Successfully restored fit from n4_bkg.h5                              time_series_builder.py:171

         INFO      Interval set to 1.28-64.257 for n4                                    time_series_builder.py:290

         INFO      Auto-probed noise models:                                                    SpectrumLike.py:490

         INFO      - observation: poisson                                                       SpectrumLike.py:491

         INFO      - background: gaussian                                                       SpectrumLike.py:492

         INFO      Range 9-900 translates to channels 5-123                                    SpectrumLike.py:1247

         INFO      Now using 119 bins                                                          SpectrumLike.py:1739

00:19:58 INFO      Auto-determined polynomial order: 1                                binned_spectrum_series.py:389

00:20:14 INFO      None 1-order polynomial fit with the mle method                               time_series.py:458

         INFO      Saved fitted background to b0_bkg.h5                                         time_series.py:1064

         INFO      Saved background to b0_bkg.h5                                         time_series_builder.py:471

00:20:15 INFO      Successfully restored fit from b0_bkg.h5                              time_series_builder.py:171

         INFO      Interval set to 1.28-64.257 for b0                                    time_series_builder.py:290

00:20:16 INFO      Auto-probed noise models:                                                    SpectrumLike.py:490

         INFO      - observation: poisson                                                       SpectrumLike.py:491

         INFO      - background: gaussian                                                       SpectrumLike.py:492

         INFO      Range 250-30000 translates to channels 1-119                                SpectrumLike.py:1247

         INFO      Now using 119 bins                                                          SpectrumLike.py:1739

../_images/notebooks_grb080916C_12_42.png

../_images/notebooks_grb080916C_12_43.png

../_images/notebooks_grb080916C_12_44.png

Setting up the fit

Let’s see if we can reproduce the results from the catalog.

Set priors for the model

We will fit the spectrum using Bayesian analysis, so we must set priors on the model parameters.

[9]:

model.GRB080916009.spectrum.main.shape.alpha.prior = Truncated_gaussian(
    lower_bound=-1.5, upper_bound=1, mu=-1, sigma=0.5
)
model.GRB080916009.spectrum.main.shape.beta.prior = Truncated_gaussian(
    lower_bound=-5, upper_bound=-1.6, mu=-2.25, sigma=0.5
)
model.GRB080916009.spectrum.main.shape.break_energy.prior = Log_normal(mu=2, sigma=1)
model.GRB080916009.spectrum.main.shape.break_energy.bounds = (None, None)
model.GRB080916009.spectrum.main.shape.K.prior = Log_uniform_prior(
    lower_bound=1e-3, upper_bound=1e1
)
model.GRB080916009.spectrum.main.shape.break_scale.prior = Log_uniform_prior(
    lower_bound=1e-4, upper_bound=10
)

Clone the model and setup the Bayesian analysis class

Next, we clone the model we built from the catalog so that we can look at the results later and fit the cloned model. We pass this model and the DataList of the plugins to a BayesianAnalysis class and set the sampler to MultiNest.

[10]:

new_model = clone_model(model)

bayes = BayesianAnalysis(new_model, DataList(*fluence_plugins))

# share spectrum gives a linear speed up when
# spectrumlike plugins have the same RSP input energies
bayes.set_sampler("multinest", share_spectrum=True)

         INFO      sampler set to multinest                                                bayesian_analysis.py:202

Examine at the catalog fitted model

We can quickly examine how well the catalog fit matches the data. There appears to be a discrepancy between the data and the model! Let’s refit to see if we can fix it.

[11]:

fig = display_spectrum_model_counts(bayes, min_rate=20, step=False)

../_images/notebooks_grb080916C_18_0.png

Run the sampler

We let MultiNest condition the model on the data

[12]:

bayes.sampler.setup(n_live_points=400)
bayes.sample()

 *****************************************************
 MultiNest v3.10
 Copyright Farhan Feroz & Mike Hobson
 Release Jul 2015

 no. of live points =  400
 dimensionality =    5
 *****************************************************
  analysing data from chains/fit-.txt ln(ev)=  -3100.9725311704797      +/-  0.22304480182573688
 Total Likelihood Evaluations:        25478
 Sampling finished. Exiting MultiNest

00:20:38 INFO      fit restored to maximum of posterior                                         sampler_base.py:178

         INFO      fit restored to maximum of posterior                                         sampler_base.py:178

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
GRB080916009...K	(1.464 -0.013 +0.024) x 10^-2	1 / (cm2 keV s)
GRB080916009...alpha	-1.097 +0.004 +0.04
GRB080916009...break_energy	(1.95 +0.07 +0.7) x 10^2	keV
GRB080916009...break_scale	(0.0 +1.7 +3.5) x 10^-1
GRB080916009...beta	-1.983 -0.23 -0.030

Values of -log(posterior) at the minimum:

	-log(posterior)
b0	-1050.585751
n3	-1019.624748
n4	-1010.768624
total	-3080.979123

Values of statistical measures:

	statistical measures
AIC	6172.128701
BIC	6191.360911
DIC	6180.202740
PDIC	4.449591
log(Z)	-1346.735259

Now our model seems to match much better with the data!

[13]:

bayes.restore_median_fit()
fig = display_spectrum_model_counts(bayes, min_rate=20)

         INFO      fit restored to median of posterior                                          sampler_base.py:164

../_images/notebooks_grb080916C_22_1.png

But how different are we from the catalog model? Let’s plot our fit along with the catalog model. Luckily, 3ML can handle all the units for is

[14]:

conversion = u.Unit("keV2/(cm2 s keV)").to("erg2/(cm2 s keV)")
energy_grid = np.logspace(1, 4, 100) * u.keV
vFv = (energy_grid**2 * model.get_point_source_fluxes(0, energy_grid)).to(
    "erg2/(cm2 s keV)"
)

[15]:

fig = plot_spectra(bayes.results, flux_unit="erg2/(cm2 s keV)")
ax = fig.get_axes()[0]
_ = ax.loglog(energy_grid, vFv, color="blue", label="catalog model")

../_images/notebooks_grb080916C_25_2.png

Time Resolved Analysis

Now that we have examined fluence fit, we can move to performing a time-resolved analysis.

Selecting a temporal binning

We first get the brightest NaI detector and create time bins via the Bayesian blocks algorithm. We can use the fitted background to make sure that our intervals are chosen in an unbiased way.

[16]:

n3 = time_series["n3"]

[17]:

n3.create_time_bins(0, 60, method="bayesblocks", use_background=True, p0=0.2)

00:23:01 INFO      Created 15 bins via bayesblocks                                       time_series_builder.py:708

Sometimes, glitches in the GBM data cause spikes in the data that the Bayesian blocks algorithm detects as fast changes in the count rate. We will have to remove those intervals manually.

Note: In the future, 3ML will provide an automated method to remove these unwanted spikes.

[18]:

fig = n3.view_lightcurve(use_binner=True)

../_images/notebooks_grb080916C_30_0.png

[19]:

bad_bins = []
for i, w in enumerate(n3.bins.widths):
    if w < 5e-2:
        bad_bins.append(i)


edges = [n3.bins.starts[0]]

for i, b in enumerate(n3.bins):
    if i not in bad_bins:
        edges.append(b.stop)

starts = edges[:-1]
stops = edges[1:]


n3.create_time_bins(starts, stops, method="custom")

00:23:02 INFO      Created 12 bins via custom                                            time_series_builder.py:708

Now our light curve looks much more acceptable.

[20]:

fig = n3.view_lightcurve(use_binner=True)

../_images/notebooks_grb080916C_33_0.png

The time series objects can read time bins from each other, so we will map these time bins onto the other detectors’ time series and create a list of time plugins for each detector and each time bin created above.

[21]:

time_resolved_plugins = {}

for k, v in time_series.items():
    v.read_bins(n3)
    time_resolved_plugins[k] = v.to_spectrumlike(from_bins=True)

         INFO      Created 12 bins via custom                                            time_series_builder.py:708

00:23:03 INFO      Interval set to 1.28-64.257 for n3                                    time_series_builder.py:290

         INFO      Created 12 bins via custom                                            time_series_builder.py:708

         INFO      Interval set to 1.28-64.257 for n4                                    time_series_builder.py:290

         INFO      Created 12 bins via custom                                            time_series_builder.py:708

00:23:05 INFO      Interval set to 1.28-64.257 for b0                                    time_series_builder.py:290

Setting up the model

For the time-resolved analysis, we will fit the classic Band function to the data. We will set some principled priors.

[22]:

band = Band()
band.alpha.prior = Truncated_gaussian(lower_bound=-1.5, upper_bound=1, mu=-1, sigma=0.5)
band.beta.prior = Truncated_gaussian(lower_bound=-5, upper_bound=-1.6, mu=-2, sigma=0.5)
band.xp.prior = Log_normal(mu=2, sigma=1)
band.xp.bounds = (None, None)
band.K.prior = Log_uniform_prior(lower_bound=1e-10, upper_bound=1e3)
ps = PointSource("grb", 0, 0, spectral_shape=band)
band_model = Model(ps)

Perform the fits

One way to perform Bayesian spectral fits to all the intervals is to loop through each one. There can are many ways to do this, so find an analysis pattern that works for you.

[23]:

models = []
results = []
analysis = []
for interval in range(12):
    # clone the model above so that we have a separate model
    # for each fit

    this_model = clone_model(band_model)

    # for each detector set up the plugin
    # for this time interval

    this_data_list = []
    for k, v in time_resolved_plugins.items():
        pi = v[interval]

        if k.startswith("b"):
            pi.set_active_measurements("250-30000")
        else:
            pi.set_active_measurements("9-900")

        pi.rebin_on_background(1.0)

        this_data_list.append(pi)

    # create a data list

    dlist = DataList(*this_data_list)

    # set up the sampler and fit

    bayes = BayesianAnalysis(this_model, dlist)

    # get some speed with share spectrum
    bayes.set_sampler("multinest", share_spectrum=True)
    bayes.sampler.setup(n_live_points=500)
    bayes.sample()

    # at this stage we coudl also
    # save the analysis result to
    # disk but we will simply hold
    # onto them in memory

    analysis.append(bayes)

         INFO      Range 9-900 translates to channels 5-124                                    SpectrumLike.py:1247

         INFO      Now using 120 bins                                                          SpectrumLike.py:1739

         INFO      Range 9-900 translates to channels 5-123                                    SpectrumLike.py:1247

         INFO      Now using 119 bins                                                          SpectrumLike.py:1739

         INFO      Range 250-30000 translates to channels 1-119                                SpectrumLike.py:1247

         INFO      Now using 107 bins                                                          SpectrumLike.py:1739

         INFO      sampler set to multinest                                                bayesian_analysis.py:202

 *****************************************************
 MultiNest v3.10
 Copyright Farhan Feroz & Mike Hobson
 Release Jul 2015

 no. of live points =  500
 dimensionality =    4
 *****************************************************
  analysing data from chains/fit-.txt ln(ev)=  -789.41846841416634      +/-  0.18048268178896379
 Total Likelihood Evaluations:        16084
 Sampling finished. Exiting MultiNest

00:23:18 INFO      fit restored to maximum of posterior                                         sampler_base.py:178

         INFO      fit restored to maximum of posterior                                         sampler_base.py:178

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
grb.spectrum.main.Band.K	(3.8 -0.9 +0.6) x 10^-2	1 / (cm2 keV s)
grb.spectrum.main.Band.alpha	(-5.1 -2.2 +1.0) x 10^-1
grb.spectrum.main.Band.xp	(3.0 -0.5 +1.3) x 10^2	keV
grb.spectrum.main.Band.beta	-2.06 -0.6 +0.19

Values of -log(posterior) at the minimum:

	-log(posterior)
b0_interval0	-285.593580
n3_interval0	-250.031211
n4_interval0	-268.019015
total	-803.643806

Values of statistical measures:

	statistical measures
AIC	1615.400925
BIC	1630.809743
DIC	1570.600468
PDIC	1.760517
log(Z)	-342.840085

         INFO      Range 9-900 translates to channels 5-124                                    SpectrumLike.py:1247

         INFO      Now using 120 bins                                                          SpectrumLike.py:1739

         INFO      Range 9-900 translates to channels 5-123                                    SpectrumLike.py:1247

         INFO      Now using 119 bins                                                          SpectrumLike.py:1739

         INFO      Range 250-30000 translates to channels 1-119                                SpectrumLike.py:1247

         INFO      Now using 119 bins                                                          SpectrumLike.py:1739

         INFO      sampler set to multinest                                                bayesian_analysis.py:202

 *****************************************************
 MultiNest v3.10
 Copyright Farhan Feroz & Mike Hobson
 Release Jul 2015

 no. of live points =  500
 dimensionality =    4
 *****************************************************
  analysing data from chains/fit-.txt ln(ev)=  -1963.8093317326559      +/-  0.22349390491567409
 Total Likelihood Evaluations:        21371
 Sampling finished. Exiting MultiNest

00:23:36 INFO      fit restored to maximum of posterior                                         sampler_base.py:178

         INFO      fit restored to maximum of posterior                                         sampler_base.py:178

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
grb.spectrum.main.Band.K	(5.192 +0.010 +0.08) x 10^-2	1 / (cm2 keV s)
grb.spectrum.main.Band.alpha	(-6.973 -0.21 +0.005) x 10^-1
grb.spectrum.main.Band.xp	(3.803 -0.10 +0.011) x 10^2	keV
grb.spectrum.main.Band.beta	-1.902 -0.021 +0.04

Values of -log(posterior) at the minimum:

	-log(posterior)
b0_interval1	-682.222264
n3_interval1	-646.169508
n4_interval1	-643.547594
total	-1971.939366

Values of statistical measures:

	statistical measures
AIC	3951.992046
BIC	3967.400864
DIC	3905.945476
PDIC	4.262979
log(Z)	-852.871556

         INFO      Range 9-900 translates to channels 5-124                                    SpectrumLike.py:1247

         INFO      Now using 120 bins                                                          SpectrumLike.py:1739

         INFO      Range 9-900 translates to channels 5-123                                    SpectrumLike.py:1247

00:23:37 INFO      Now using 119 bins                                                          SpectrumLike.py:1739

         INFO      Range 250-30000 translates to channels 1-119                                SpectrumLike.py:1247

         INFO      Now using 115 bins                                                          SpectrumLike.py:1739

         INFO      sampler set to multinest                                                bayesian_analysis.py:202

 *****************************************************
 MultiNest v3.10
 Copyright Farhan Feroz & Mike Hobson
 Release Jul 2015

 no. of live points =  500
 dimensionality =    4
 *****************************************************
  analysing data from chains/fit-.txt ln(ev)=  -910.81739038145872      +/-  0.19141821890680788
 Total Likelihood Evaluations:        19270
 Sampling finished. Exiting MultiNest

00:23:55 INFO      fit restored to maximum of posterior                                         sampler_base.py:178

         INFO      fit restored to maximum of posterior                                         sampler_base.py:178

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
grb.spectrum.main.Band.K	(2.86 +0.07 +0.8) x 10^-2	1 / (cm2 keV s)
grb.spectrum.main.Band.alpha	(-9.82 -0.26 +1.3) x 10^-1
grb.spectrum.main.Band.xp	(3.9 -1.4 -0.4) x 10^2	keV
grb.spectrum.main.Band.beta	-1.667 -0.04 +0.035

Values of -log(posterior) at the minimum:

	-log(posterior)
b0_interval2	-324.534988
n3_interval2	-288.070103
n4_interval2	-310.970190
total	-923.575281

Values of statistical measures:

	statistical measures
AIC	1855.263877
BIC	1870.672695
DIC	1807.386314
PDIC	1.876346
log(Z)	-395.562967

         INFO      Range 9-900 translates to channels 5-124                                    SpectrumLike.py:1247

         INFO      Now using 120 bins                                                          SpectrumLike.py:1739

         INFO      Range 9-900 translates to channels 5-123                                    SpectrumLike.py:1247

         INFO      Now using 119 bins                                                          SpectrumLike.py:1739

         INFO      Range 250-30000 translates to channels 1-119                                SpectrumLike.py:1247

         INFO      Now using 109 bins                                                          SpectrumLike.py:1739

         INFO      sampler set to multinest                                                bayesian_analysis.py:202

 *****************************************************
 MultiNest v3.10
 Copyright Farhan Feroz & Mike Hobson
 Release Jul 2015

 no. of live points =  500
 dimensionality =    4
 *****************************************************
  analysing data from chains/fit-.txt ln(ev)=  -788.82547684786982      +/-  0.17777153569790402
 Total Likelihood Evaluations:        17581
 Sampling finished. Exiting MultiNest

00:24:12 INFO      fit restored to maximum of posterior                                         sampler_base.py:178

         INFO      fit restored to maximum of posterior                                         sampler_base.py:178

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
grb.spectrum.main.Band.K	(2.85 -0.28 +0.4) x 10^-2	1 / (cm2 keV s)
grb.spectrum.main.Band.alpha	(-9.4 -0.8 +1.1) x 10^-1
grb.spectrum.main.Band.xp	(3.4 +/- 0.7) x 10^2	keV
grb.spectrum.main.Band.beta	-2.248 -0.012 +0.31

Values of -log(posterior) at the minimum:

	-log(posterior)
b0_interval3	-298.472349
n3_interval3	-242.578577
n4_interval3	-262.511897
total	-803.562822

Values of statistical measures:

	statistical measures
AIC	1615.238959
BIC	1630.647777
DIC	1569.942822
PDIC	2.625488
log(Z)	-342.582552

         INFO      Range 9-900 translates to channels 5-124                                    SpectrumLike.py:1247

         INFO      Now using 120 bins                                                          SpectrumLike.py:1739

         INFO      Range 9-900 translates to channels 5-123                                    SpectrumLike.py:1247

         INFO      Now using 119 bins                                                          SpectrumLike.py:1739

         INFO      Range 250-30000 translates to channels 1-119                                SpectrumLike.py:1247

         INFO      Now using 119 bins                                                          SpectrumLike.py:1739

         INFO      sampler set to multinest                                                bayesian_analysis.py:202

 *****************************************************
 MultiNest v3.10
 Copyright Farhan Feroz & Mike Hobson
 Release Jul 2015

 no. of live points =  500
 dimensionality =    4
 *****************************************************
  analysing data from chains/fit-.txt ln(ev)=  -2279.8055153561954      +/-  0.20854959376331531
 Total Likelihood Evaluations:        19692
 Sampling finished. Exiting MultiNest

00:24:34 INFO      fit restored to maximum of posterior                                         sampler_base.py:178

         INFO      fit restored to maximum of posterior                                         sampler_base.py:178

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
grb.spectrum.main.Band.K	(2.4608 +0.0013 +0.19) x 10^-2	1 / (cm2 keV s)
grb.spectrum.main.Band.alpha	(-8.68 +0.05 +0.4) x 10^-1
grb.spectrum.main.Band.xp	(2.817 -0.33 -0.024) x 10^2	keV
grb.spectrum.main.Band.beta	-1.792 -0.008 +0.018

Values of -log(posterior) at the minimum:

	-log(posterior)
b0_interval4	-779.708451
n3_interval4	-758.818561
n4_interval4	-745.911150
total	-2284.438163

Values of statistical measures:

	statistical measures
AIC	4576.989640
BIC	4592.398457
DIC	4537.030468
PDIC	2.163925
log(Z)	-990.106955

         INFO      Range 9-900 translates to channels 5-124                                    SpectrumLike.py:1247

         INFO      Now using 120 bins                                                          SpectrumLike.py:1739

         INFO      Range 9-900 translates to channels 5-123                                    SpectrumLike.py:1247

         INFO      Now using 119 bins                                                          SpectrumLike.py:1739

         INFO      Range 250-30000 translates to channels 1-119                                SpectrumLike.py:1247

         INFO      Now using 119 bins                                                          SpectrumLike.py:1739

         INFO      sampler set to multinest                                                bayesian_analysis.py:202

 *****************************************************
 MultiNest v3.10
 Copyright Farhan Feroz & Mike Hobson
 Release Jul 2015

 no. of live points =  500
 dimensionality =    4
 *****************************************************
  analysing data from chains/fit-.txt ln(ev)=  -1575.1284457690995      +/-  0.19547778320900869
 Total Likelihood Evaluations:        20748
 Sampling finished. Exiting MultiNest

00:24:57 INFO      fit restored to maximum of posterior                                         sampler_base.py:178

         INFO      fit restored to maximum of posterior                                         sampler_base.py:178

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
grb.spectrum.main.Band.K	(2.85 +0.07 +0.35) x 10^-2	1 / (cm2 keV s)
grb.spectrum.main.Band.alpha	(-9.05 +0.29 +1.0) x 10^-1
grb.spectrum.main.Band.xp	(3.97 -0.6 -0.04) x 10^2	keV
grb.spectrum.main.Band.beta	-2.06 -0.23 +0.04

Values of -log(posterior) at the minimum:

	-log(posterior)
b0_interval5	-536.382595
n3_interval5	-523.310250
n4_interval5	-527.894502
total	-1587.587347

Values of statistical measures:

	statistical measures
AIC	3183.288009
BIC	3198.696826
DIC	3137.057165
PDIC	2.849895
log(Z)	-684.069592

         INFO      Range 9-900 translates to channels 5-124                                    SpectrumLike.py:1247

         INFO      Now using 120 bins                                                          SpectrumLike.py:1739

         INFO      Range 9-900 translates to channels 5-123                                    SpectrumLike.py:1247

         INFO      Now using 119 bins                                                          SpectrumLike.py:1739

         INFO      Range 250-30000 translates to channels 1-119                                SpectrumLike.py:1247

         INFO      Now using 119 bins                                                          SpectrumLike.py:1739

         INFO      sampler set to multinest                                                bayesian_analysis.py:202

 *****************************************************
 MultiNest v3.10
 Copyright Farhan Feroz & Mike Hobson
 Release Jul 2015

 no. of live points =  500
 dimensionality =    4
 *****************************************************
  analysing data from chains/fit-.txt ln(ev)=  -1754.6463999276384      +/-  0.18954879192016755
 Total Likelihood Evaluations:        21916
 Sampling finished. Exiting MultiNest

00:25:16 INFO      fit restored to maximum of posterior                                         sampler_base.py:178

         INFO      fit restored to maximum of posterior                                         sampler_base.py:178

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
grb.spectrum.main.Band.K	(1.97 -0.15 +0.09) x 10^-2	1 / (cm2 keV s)
grb.spectrum.main.Band.alpha	-1.004 -0.06 +0.034
grb.spectrum.main.Band.xp	(4.4 -0.4 +1.0) x 10^2	keV
grb.spectrum.main.Band.beta	-2.34 -0.4 +0.15

Values of -log(posterior) at the minimum:

	-log(posterior)
b0_interval6	-609.465978
n3_interval6	-584.286903
n4_interval6	-576.734338
total	-1770.487219

Values of statistical measures:

	statistical measures
AIC	3549.087753
BIC	3564.496570
DIC	3501.676151
PDIC	3.549509
log(Z)	-762.033249

         INFO      Range 9-900 translates to channels 5-124                                    SpectrumLike.py:1247

         INFO      Now using 120 bins                                                          SpectrumLike.py:1739

         INFO      Range 9-900 translates to channels 5-123                                    SpectrumLike.py:1247

         INFO      Now using 119 bins                                                          SpectrumLike.py:1739

         INFO      Range 250-30000 translates to channels 1-119                                SpectrumLike.py:1247

         INFO      Now using 119 bins                                                          SpectrumLike.py:1739

         INFO      sampler set to multinest                                                bayesian_analysis.py:202

 *****************************************************
 MultiNest v3.10
 Copyright Farhan Feroz & Mike Hobson
 Release Jul 2015

 no. of live points =  500
 dimensionality =    4
 *****************************************************
  analysing data from chains/fit-.txt ln(ev)=  -1939.3669472925658      +/-  0.19251446728097965
 Total Likelihood Evaluations:        19293
 Sampling finished. Exiting MultiNest

00:25:30 INFO      fit restored to maximum of posterior                                         sampler_base.py:178

         INFO      fit restored to maximum of posterior                                         sampler_base.py:178

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
grb.spectrum.main.Band.K	(1.65 -0.09 +0.11) x 10^-2	1 / (cm2 keV s)
grb.spectrum.main.Band.alpha	-1.05 +/- 0.05
grb.spectrum.main.Band.xp	(4.4 -0.6 +0.7) x 10^2	keV
grb.spectrum.main.Band.beta	-2.25 -0.4 +0.12

Values of -log(posterior) at the minimum:

	-log(posterior)
b0_interval7	-662.466489
n3_interval7	-640.654608
n4_interval7	-650.705087
total	-1953.826183

Values of statistical measures:

	statistical measures
AIC	3915.765681
BIC	3931.174499
DIC	3868.474468
PDIC	3.175189
log(Z)	-842.256364

         INFO      Range 9-900 translates to channels 5-124                                    SpectrumLike.py:1247

         INFO      Now using 120 bins                                                          SpectrumLike.py:1739

         INFO      Range 9-900 translates to channels 5-123                                    SpectrumLike.py:1247

         INFO      Now using 119 bins                                                          SpectrumLike.py:1739

         INFO      Range 250-30000 translates to channels 1-119                                SpectrumLike.py:1247

         INFO      Now using 119 bins                                                          SpectrumLike.py:1739

         INFO      sampler set to multinest                                                bayesian_analysis.py:202

 *****************************************************
 MultiNest v3.10
 Copyright Farhan Feroz & Mike Hobson
 Release Jul 2015

 no. of live points =  500
 dimensionality =    4
 *****************************************************
  analysing data from chains/fit-.txt ln(ev)=  -2057.8046955957693      +/-  0.19555872563174734
 Total Likelihood Evaluations:        19239
 Sampling finished. Exiting MultiNest

00:25:47 INFO      fit restored to maximum of posterior                                         sampler_base.py:178

         INFO      fit restored to maximum of posterior                                         sampler_base.py:178

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
grb.spectrum.main.Band.K	(1.6626 +0.0021 +0.22) x 10^-2	1 / (cm2 keV s)
grb.spectrum.main.Band.alpha	(-7.98 -0.04 +0.8) x 10^-1
grb.spectrum.main.Band.xp	(3.262 -0.5 -0.005) x 10^2	keV
grb.spectrum.main.Band.beta	-2.102 -0.015 +0.14

Values of -log(posterior) at the minimum:

	-log(posterior)
b0_interval8	-701.913282
n3_interval8	-698.746276
n4_interval8	-665.340467
total	-2066.000025

Values of statistical measures:

	statistical measures
AIC	4140.113365
BIC	4155.522182
DIC	4099.179210
PDIC	2.551884
log(Z)	-893.693224

         INFO      Range 9-900 translates to channels 5-124                                    SpectrumLike.py:1247

         INFO      Now using 120 bins                                                          SpectrumLike.py:1739

         INFO      Range 9-900 translates to channels 5-123                                    SpectrumLike.py:1247

         INFO      Now using 119 bins                                                          SpectrumLike.py:1739

         INFO      Range 250-30000 translates to channels 1-119                                SpectrumLike.py:1247

         INFO      Now using 119 bins                                                          SpectrumLike.py:1739

         INFO      sampler set to multinest                                                bayesian_analysis.py:202

 *****************************************************
 MultiNest v3.10
 Copyright Farhan Feroz & Mike Hobson
 Release Jul 2015

 no. of live points =  500
 dimensionality =    4
 *****************************************************
  analysing data from chains/fit-.txt ln(ev)=  -1878.4244611544707      +/-  0.14435419698869176
 Total Likelihood Evaluations:        13514
 Sampling finished. Exiting MultiNest

00:25:57 INFO      fit restored to maximum of posterior                                         sampler_base.py:178

         INFO      fit restored to maximum of posterior                                         sampler_base.py:178

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
grb.spectrum.main.Band.K	(1.2 -0.5 +0.9) x 10^-2	1 / (cm2 keV s)
grb.spectrum.main.Band.alpha	(-8.2 -2.4 +2.5) x 10^-1
grb.spectrum.main.Band.xp	(1.06 -0.32 +0.6) x 10^2	keV
grb.spectrum.main.Band.beta	-1.82 -0.4 +0.09

Values of -log(posterior) at the minimum:

	-log(posterior)
b0_interval9	-648.509791
n3_interval9	-616.786870
n4_interval9	-616.117301
total	-1881.413962

Values of statistical measures:

	statistical measures
AIC	3770.941239
BIC	3786.350056
DIC	3708.304695
PDIC	-38.542167
log(Z)	-815.789378

00:25:58 INFO      Range 9-900 translates to channels 5-124                                    SpectrumLike.py:1247

         INFO      Now using 120 bins                                                          SpectrumLike.py:1739

         INFO      Range 9-900 translates to channels 5-123                                    SpectrumLike.py:1247

         INFO      Now using 119 bins                                                          SpectrumLike.py:1739

         INFO      Range 250-30000 translates to channels 1-119                                SpectrumLike.py:1247

         INFO      Now using 119 bins                                                          SpectrumLike.py:1739

         INFO      sampler set to multinest                                                bayesian_analysis.py:202

 *****************************************************
 MultiNest v3.10
 Copyright Farhan Feroz & Mike Hobson
 Release Jul 2015

 no. of live points =  500
 dimensionality =    4
 *****************************************************
  analysing data from chains/fit-.txt ln(ev)=  -1322.0428302282203      +/-  0.16636117827400865
 Total Likelihood Evaluations:        16201
 Sampling finished. Exiting MultiNest

00:26:09 INFO      fit restored to maximum of posterior                                         sampler_base.py:178

         INFO      fit restored to maximum of posterior                                         sampler_base.py:178

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
grb.spectrum.main.Band.K	(2.1 +/- 0.5) x 10^-2	1 / (cm2 keV s)
grb.spectrum.main.Band.alpha	(-7.2 -1.5 +1.4) x 10^-1
grb.spectrum.main.Band.xp	(2.1 -0.4 +0.7) x 10^2	keV
grb.spectrum.main.Band.beta	-1.94 -0.5 +0.11

Values of -log(posterior) at the minimum:

	-log(posterior)
b0_interval10	-460.740590
n3_interval10	-437.815150
n4_interval10	-432.948397
total	-1331.504136

Values of statistical measures:

	statistical measures
AIC	2671.121587
BIC	2686.530405
DIC	2633.873210
PDIC	0.029746
log(Z)	-574.155906

         INFO      Range 9-900 translates to channels 5-124                                    SpectrumLike.py:1247

         INFO      Now using 120 bins                                                          SpectrumLike.py:1739

         INFO      Range 9-900 translates to channels 5-123                                    SpectrumLike.py:1247

         INFO      Now using 119 bins                                                          SpectrumLike.py:1739

         INFO      Range 250-30000 translates to channels 1-119                                SpectrumLike.py:1247

         INFO      Now using 119 bins                                                          SpectrumLike.py:1739

         INFO      sampler set to multinest                                                bayesian_analysis.py:202

 *****************************************************
 MultiNest v3.10
 Copyright Farhan Feroz & Mike Hobson
 Release Jul 2015

 no. of live points =  500
 dimensionality =    4
 *****************************************************
  analysing data from chains/fit-.txt ln(ev)=  -813.99277406797080      +/-  0.15434756832240351
 Total Likelihood Evaluations:        12515
 Sampling finished. Exiting MultiNest

00:26:19 INFO      fit restored to maximum of posterior                                         sampler_base.py:178

         INFO      fit restored to maximum of posterior                                         sampler_base.py:178

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
grb.spectrum.main.Band.K	(2.7 -0.6 +4) x 10^-2	1 / (cm2 keV s)
grb.spectrum.main.Band.alpha	(-5.0 -2.1 +4) x 10^-1
grb.spectrum.main.Band.xp	(1.29 -0.5 +0.08) x 10^2	keV
grb.spectrum.main.Band.beta	-2.04 -0.10 +0.28

Values of -log(posterior) at the minimum:

	-log(posterior)
b0_interval11	-292.296133
n3_interval11	-272.386209
n4_interval11	-255.810193
total	-820.492534

Values of statistical measures:

	statistical measures
AIC	1649.098383
BIC	1664.507200
DIC	1610.994972
PDIC	-6.504105
log(Z)	-353.512570

Examine the fits

Now we can look at the fits in count space to make sure they are ok.

[24]:

for a in analysis:
    a.restore_median_fit()
    _ = display_spectrum_model_counts(a, min_rate=[20, 20, 20], step=False)

         INFO      fit restored to median of posterior                                          sampler_base.py:164

00:26:20 INFO      fit restored to median of posterior                                          sampler_base.py:164

00:26:21 INFO      fit restored to median of posterior                                          sampler_base.py:164

         INFO      fit restored to median of posterior                                          sampler_base.py:164

00:26:22 INFO      fit restored to median of posterior                                          sampler_base.py:164

00:26:23 INFO      fit restored to median of posterior                                          sampler_base.py:164

00:26:24 INFO      fit restored to median of posterior                                          sampler_base.py:164

         INFO      fit restored to median of posterior                                          sampler_base.py:164

00:26:25 INFO      fit restored to median of posterior                                          sampler_base.py:164

00:26:26 INFO      fit restored to median of posterior                                          sampler_base.py:164

         INFO      fit restored to median of posterior                                          sampler_base.py:164

00:26:27 INFO      fit restored to median of posterior                                          sampler_base.py:164

../_images/notebooks_grb080916C_41_12.png

../_images/notebooks_grb080916C_41_13.png

../_images/notebooks_grb080916C_41_14.png

../_images/notebooks_grb080916C_41_15.png

../_images/notebooks_grb080916C_41_16.png

../_images/notebooks_grb080916C_41_17.png

../_images/notebooks_grb080916C_41_18.png

../_images/notebooks_grb080916C_41_19.png

../_images/notebooks_grb080916C_41_20.png

../_images/notebooks_grb080916C_41_21.png

../_images/notebooks_grb080916C_41_22.png

../_images/notebooks_grb080916C_41_23.png

Finally, we can plot the models together to see how the spectra evolve with time.

[25]:

fig = plot_spectra(
    *[a.results for a in analysis[::1]],
    flux_unit="erg2/(cm2 s keV)",
    fit_cmap="viridis",
    contour_cmap="viridis",
    contour_style_kwargs=dict(alpha=0.1),
)

../_images/notebooks_grb080916C_43_13.png

This example can serve as a template for performing analysis on GBM data. However, as 3ML provides an abstract interface and modular building blocks, similar analysis pipelines can be built for any time series data.