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")

18:44:06 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)

18:45:42 INFO      Auto-determined polynomial order: 0                                binned_spectrum_series.py:389

18:45:51 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

18:45:52 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

18:45:54 INFO      Now using 120 bins                                                          SpectrumLike.py:1739

18:45:56 INFO      Auto-determined polynomial order: 1                                binned_spectrum_series.py:389

18:46:05 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

18:46:06 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

18:46:07 INFO      Auto-determined polynomial order: 1                                binned_spectrum_series.py:389

18:46:18 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

         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

18:46:19 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)

18:46:20 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)=  -3103.0077213195996      +/-  0.23614490739484614
 Total Likelihood Evaluations:        23357
 Sampling finished. Exiting MultiNest

18:46:39 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.474 -0.010 +0.018) x 10^-2	1 / (cm2 keV s)
GRB080916009...alpha	-1.078 -0.008 +0.022
GRB080916009...break_energy	(2.14 -0.16 +0.11) x 10^2	keV
GRB080916009...break_scale	(2.3 -0.4 +0.6) x 10^-1
GRB080916009...beta	-2.05 -0.04 +0.06

Values of -log(posterior) at the minimum:

	-log(posterior)
b0	-1055.400631
n3	-1025.586097
n4	-1017.443989
total	-3098.430717

Values of statistical measures:

	statistical measures
AIC	6207.031889
BIC	6226.264099
DIC	6177.672528
PDIC	3.278476
log(Z)	-1347.619131

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)

18:48:55 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")

18:48:56 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

         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

18:48:57 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

18:48:58 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)=  -788.14561911054818      +/-  0.17479123791490644
 Total Likelihood Evaluations:        16835
 Sampling finished. Exiting MultiNest

18:49:13 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.6 +/- 0.6) x 10^-2	1 / (cm2 keV s)
grb.spectrum.main.Band.alpha	(-5.5 -1.5 +1.2) x 10^-1
grb.spectrum.main.Band.xp	(3.2 -0.4 +0.9) x 10^2	keV
grb.spectrum.main.Band.beta	-2.05 -0.5 +0.10

Values of -log(posterior) at the minimum:

	-log(posterior)
b0_interval0	-285.656120
n3_interval0	-250.235678
n4_interval0	-268.033190
total	-803.924987

Values of statistical measures:

	statistical measures
AIC	1615.963289
BIC	1631.372106
DIC	1570.813562
PDIC	2.287559
log(Z)	-342.287293

         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)=  -1943.4773096929068      +/-  0.21388511217313472
 Total Likelihood Evaluations:        24334
 Sampling finished. Exiting MultiNest

18:49:32 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	(4.11 -0.09 +0.17) x 10^-2	1 / (cm2 keV s)
grb.spectrum.main.Band.alpha	(-8.59 -0.20 +0.29) x 10^-1
grb.spectrum.main.Band.xp	(6.29 -0.6 +0.33) x 10^2	keV
grb.spectrum.main.Band.beta	-2.14 -0.07 +0.08

Values of -log(posterior) at the minimum:

	-log(posterior)
b0_interval1	-673.933346
n3_interval1	-641.728738
n4_interval1	-645.630544
total	-1961.292628

Values of statistical measures:

	statistical measures
AIC	3930.698571
BIC	3946.107389
DIC	3872.732801
PDIC	3.401111
log(Z)	-844.041471

         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 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.09441746626987      +/-  0.18777394544406945
 Total Likelihood Evaluations:        19059
 Sampling finished. Exiting MultiNest

18:49: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	(3.04 +0.08 +0.6) x 10^-2	1 / (cm2 keV s)
grb.spectrum.main.Band.alpha	(-9.29 -0.17 +1.0) x 10^-1
grb.spectrum.main.Band.xp	(3.48 -1.0 -0.16) x 10^2	keV
grb.spectrum.main.Band.beta	-1.734 -0.012 +0.09

Values of -log(posterior) at the minimum:

	-log(posterior)
b0_interval2	-323.453556
n3_interval2	-288.101920
n4_interval2	-310.858319
total	-922.413795

Values of statistical measures:

	statistical measures
AIC	1852.940905
BIC	1868.349722
DIC	1807.206940
PDIC	1.940962
log(Z)	-395.248984

         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.52166806866194      +/-  0.17842334049526717
 Total Likelihood Evaluations:        16874
 Sampling finished. Exiting MultiNest

18:50:02 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.81 -0.25 +0.5) x 10^-2	1 / (cm2 keV s)
grb.spectrum.main.Band.alpha	(-9.4 -0.8 +0.9) x 10^-1
grb.spectrum.main.Band.xp	(3.6 -0.8 +0.6) x 10^2	keV
grb.spectrum.main.Band.beta	-2.24 -0.4 +0.20

Values of -log(posterior) at the minimum:

	-log(posterior)
b0_interval3	-298.565905
n3_interval3	-242.676066
n4_interval3	-262.709417
total	-803.951389

Values of statistical measures:

	statistical measures
AIC	1616.016092
BIC	1631.424909
DIC	1570.313743
PDIC	2.878773
log(Z)	-342.450609

         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)=  -2270.8416346074600      +/-  0.19884112875948648
 Total Likelihood Evaluations:        21257
 Sampling finished. Exiting MultiNest

18:50: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	(2.04 -0.08 +0.13) x 10^-2	1 / (cm2 keV s)
grb.spectrum.main.Band.alpha	(-9.86 -0.26 +0.5) x 10^-1
grb.spectrum.main.Band.xp	(4.0 +/- 0.4) x 10^2	keV
grb.spectrum.main.Band.beta	-1.94 -0.11 +0.05

Values of -log(posterior) at the minimum:

	-log(posterior)
b0_interval4	-778.503923
n3_interval4	-756.958823
n4_interval4	-746.878415
total	-2282.341161

Values of statistical measures:

	statistical measures
AIC	4572.795637
BIC	4588.204455
DIC	4527.490815
PDIC	3.347241
log(Z)	-986.213991

18:50:19 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.0539725563331      +/-  0.19651863373558909
 Total Likelihood Evaluations:        20477
 Sampling finished. Exiting MultiNest

18:50:35 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.05 +0.33) x 10^-2	1 / (cm2 keV s)
grb.spectrum.main.Band.alpha	(-8.94 +0.11 +0.8) x 10^-1
grb.spectrum.main.Band.xp	(4.01 -0.6 -0.08) x 10^2	keV
grb.spectrum.main.Band.beta	-2.04 -0.17 +0.07

Values of -log(posterior) at the minimum:

	-log(posterior)
b0_interval5	-536.733857
n3_interval5	-523.670734
n4_interval5	-527.283745
total	-1587.688336

Values of statistical measures:

	statistical measures
AIC	3183.489986
BIC	3198.898803
DIC	3136.412396
PDIC	2.864210
log(Z)	-684.037249

         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)=  -1755.8637754415133      +/-  0.19682960866125088
 Total Likelihood Evaluations:        19534
 Sampling finished. Exiting MultiNest

18:50:49 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.09 +0.11) x 10^-2	1 / (cm2 keV s)
grb.spectrum.main.Band.alpha	-1.00 +/- 0.04
grb.spectrum.main.Band.xp	(4.4 -0.4 +0.6) x 10^2	keV
grb.spectrum.main.Band.beta	-2.303 -0.5 +0.018

Values of -log(posterior) at the minimum:

	-log(posterior)
b0_interval6	-609.402825
n3_interval6	-584.270380
n4_interval6	-576.649163
total	-1770.322368

Values of statistical measures:

	statistical measures
AIC	3548.758050
BIC	3564.166868
DIC	3500.167931
PDIC	2.644541
log(Z)	-762.561949

         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.3916734459458      +/-  0.19210514055637976
 Total Likelihood Evaluations:        19090
 Sampling finished. Exiting MultiNest

18:51:03 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.68 -0.12 +0.09) x 10^-2	1 / (cm2 keV s)
grb.spectrum.main.Band.alpha	-1.033 -0.06 +0.033
grb.spectrum.main.Band.xp	(4.3 -0.4 +0.8) x 10^2	keV
grb.spectrum.main.Band.beta	-2.28 -0.4 +0.14

Values of -log(posterior) at the minimum:

	-log(posterior)
b0_interval7	-662.247107
n3_interval7	-640.869581
n4_interval7	-650.507074
total	-1953.623762

Values of statistical measures:

	statistical measures
AIC	3915.360839
BIC	3930.769656
DIC	3868.873066
PDIC	3.351158
log(Z)	-842.267102

         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)=  -2063.7856212196361      +/-  0.20552485216321686
 Total Likelihood Evaluations:        19870
 Sampling finished. Exiting MultiNest

18:51:22 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.779 -0.007 +0.18) x 10^-2	1 / (cm2 keV s)
grb.spectrum.main.Band.alpha	(-7.56 -0.12 +0.4) x 10^-1
grb.spectrum.main.Band.xp	(2.816 -0.4 +0.029) x 10^2	keV
grb.spectrum.main.Band.beta	-1.875 +0.006 +0.034

Values of -log(posterior) at the minimum:

	-log(posterior)
b0_interval8	-703.827443
n3_interval8	-699.024630
n4_interval8	-665.717469
total	-2068.569542

Values of statistical measures:

	statistical measures
AIC	4145.252398
BIC	4160.661215
DIC	4105.303346
PDIC	1.673619
log(Z)	-896.290707

         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.9120382875662      +/-  0.14437639354503232
 Total Likelihood Evaluations:        13273
 Sampling finished. Exiting MultiNest

18:51:33 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.07 -0.15 +2.4) x 10^-2	1 / (cm2 keV s)
grb.spectrum.main.Band.alpha	(-8.8 -0.9 +5) x 10^-1
grb.spectrum.main.Band.xp	(1.17 -0.6 +0.13) x 10^2	keV
grb.spectrum.main.Band.beta	-1.89 -0.13 +0.19

Values of -log(posterior) at the minimum:

	-log(posterior)
b0_interval9	-648.404089
n3_interval9	-617.111523
n4_interval9	-616.286261
total	-1881.801873

Values of statistical measures:

	statistical measures
AIC	3771.717061
BIC	3787.125878
DIC	3691.916058
PDIC	-55.134489
log(Z)	-816.001130

         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.5230233053551      +/-  0.16835516882689766
 Total Likelihood Evaluations:        15984
 Sampling finished. Exiting MultiNest

18:51:46 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.4 +0.7) x 10^-2	1 / (cm2 keV s)
grb.spectrum.main.Band.alpha	(-7.2 +/- 1.6) x 10^-1
grb.spectrum.main.Band.xp	(2.1 -0.5 +0.6) x 10^2	keV
grb.spectrum.main.Band.beta	-1.94 -0.33 +0.12

Values of -log(posterior) at the minimum:

	-log(posterior)
b0_interval10	-460.821341
n3_interval10	-437.616896
n4_interval10	-433.191219
total	-1331.629456

Values of statistical measures:

	statistical measures
AIC	2671.372226
BIC	2686.781044
DIC	2632.086245
PDIC	-1.531309
log(Z)	-574.364451

         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)=  -811.99074192251965      +/-  0.14731152285519392
 Total Likelihood Evaluations:        11880
 Sampling finished. Exiting MultiNest

18:51:54 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.8 -0.8 +2.1) x 10^-2	1 / (cm2 keV s)
grb.spectrum.main.Band.alpha	(-4.8 -2.0 +3.1) x 10^-1
grb.spectrum.main.Band.xp	(1.29 -0.30 +0.28) x 10^2	keV
grb.spectrum.main.Band.beta	-2.04 -0.5 +0.11

Values of -log(posterior) at the minimum:

	-log(posterior)
b0_interval11	-292.358135
n3_interval11	-272.446421
n4_interval11	-255.838918
total	-820.643474

Values of statistical measures:

	statistical measures
AIC	1649.400262
BIC	1664.809080
DIC	1617.142696
PDIC	0.071654
log(Z)	-352.643099

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

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

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

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

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

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

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

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

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

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

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

18:52:04 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.