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")
19:53:15 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
nameradectrigger_timet90
objectfloat64float64float64float64
GRB080916009119.800-56.60054725.008861362.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)
19:54:15 INFO      Auto-determined polynomial order: 0                                binned_spectrum_series.py:389
19:54:31 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
19:54:32 INFO      Interval set to 1.28-64.257 for n3                                    time_series_builder.py:290
         INFO      Auto-probed noise models:                                                    SpectrumLike.py:469
         INFO      - observation: poisson                                                       SpectrumLike.py:470
         INFO      - background: gaussian                                                       SpectrumLike.py:471
         INFO      Range 9-900 translates to channels 5-124                                    SpectrumLike.py:1210
19:54:35 INFO      Now using 120 bins                                                          SpectrumLike.py:1673
19:54:37 INFO      Auto-determined polynomial order: 1                                binned_spectrum_series.py:389
19:54:54 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
19:54:55 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:469
         INFO      - observation: poisson                                                       SpectrumLike.py:470
         INFO      - background: gaussian                                                       SpectrumLike.py:471
         INFO      Range 9-900 translates to channels 5-123                                    SpectrumLike.py:1210
         INFO      Now using 119 bins                                                          SpectrumLike.py:1673
19:54:57 INFO      Auto-determined polynomial order: 1                                binned_spectrum_series.py:389
19:55: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
19:55: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
         INFO      Auto-probed noise models:                                                    SpectrumLike.py:469
         INFO      - observation: poisson                                                       SpectrumLike.py:470
         INFO      - background: gaussian                                                       SpectrumLike.py:471
         INFO      Range 250-30000 translates to channels 1-119                                SpectrumLike.py:1210
         INFO      Now using 119 bins                                                          SpectrumLike.py:1673
../_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)
19:55:16 INFO      sampler set to multinest                                                bayesian_analysis.py:197

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
 *****************************************************
 ln(ev)=  -3106.2403997395095      +/-  0.23382643502267439
 Total Likelihood Evaluations:        22056
  analysing data from chains/fit-.txt
 Sampling finished. Exiting MultiNest
19:55: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.448 +0.007 +0.04) x 10^-2 1 / (cm2 keV s)
GRB080916009...alpha -1.090 -0.010 +0.013
GRB080916009...break_energy (1.97 -0.19 -0.06) x 10^2 keV
GRB080916009...break_scale (0.01 +0.06 +1.2) x 10^-1
GRB080916009...beta -2.005 +0.022 +0.08
Values of -log(posterior) at the minimum:

-log(posterior)
b0 -1051.881177
n3 -1021.846394
n4 -1014.173979
total -3087.901550
Values of statistical measures:

statistical measures
AIC 6185.973555
BIC 6205.205766
DIC 6178.363946
PDIC 1.624398
log(Z) -1349.023065

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)
19:57:23 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")
19:57:24 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
19:57:25 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
19:57:26 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 = (0, 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:1210
         INFO      Now using 120 bins                                                          SpectrumLike.py:1673
         INFO      Range 9-900 translates to channels 5-123                                    SpectrumLike.py:1210
         INFO      Now using 119 bins                                                          SpectrumLike.py:1673
         INFO      Range 250-30000 translates to channels 1-119                                SpectrumLike.py:1210
         INFO      Now using 107 bins                                                          SpectrumLike.py:1673
         INFO      sampler set to multinest                                                bayesian_analysis.py:197
 *****************************************************
 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.09135677499012      +/-  0.17426220190747665
 Total Likelihood Evaluations:        18422
 Sampling finished. Exiting MultiNest

19:57:42 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.6 -1.2 +1.4) x 10^-1
grb.spectrum.main.Band.xp (3.1 -0.4 +0.9) x 10^2 keV
grb.spectrum.main.Band.beta -2.01 -0.5 +0.07
Values of -log(posterior) at the minimum:

-log(posterior)
b0_interval0 -285.673029
n3_interval0 -250.106253
n4_interval0 -267.907306
total -803.686588
Values of statistical measures:

statistical measures
AIC 1615.486491
BIC 1630.895309
DIC 1571.034073
PDIC 2.638142
log(Z) -342.263727
         INFO      Range 9-900 translates to channels 5-124                                    SpectrumLike.py:1210
         INFO      Now using 120 bins                                                          SpectrumLike.py:1673
         INFO      Range 9-900 translates to channels 5-123                                    SpectrumLike.py:1210
         INFO      Now using 119 bins                                                          SpectrumLike.py:1673
         INFO      Range 250-30000 translates to channels 1-119                                SpectrumLike.py:1210
         INFO      Now using 119 bins                                                          SpectrumLike.py:1673
         INFO      sampler set to multinest                                                bayesian_analysis.py:197
 *****************************************************
 MultiNest v3.10
 Copyright Farhan Feroz & Mike Hobson
 Release Jul 2015

 no. of live points =  500
 dimensionality =    4
 *****************************************************
 ln(ev)=  -1944.5379816653137      +/-  0.21914724208112216
  analysing data from chains/fit-.txt
 Total Likelihood Evaluations:        22203
 Sampling finished. Exiting MultiNest
19:58:04 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.07 -0.13 +0.09) x 10^-2 1 / (cm2 keV s)
grb.spectrum.main.Band.alpha (-8.66 -0.28 +0.17) x 10^-1
grb.spectrum.main.Band.xp (6.43 -0.30 +0.7) x 10^2 keV
grb.spectrum.main.Band.beta -2.183 -0.15 -0.011
Values of -log(posterior) at the minimum:

-log(posterior)
b0_interval1 -673.690462
n3_interval1 -642.149841
n4_interval1 -645.854999
total -1961.695302
Values of statistical measures:

statistical measures
AIC 3931.503918
BIC 3946.912736
DIC 3872.522899
PDIC 2.495586
log(Z) -844.502115
         INFO      Range 9-900 translates to channels 5-124                                    SpectrumLike.py:1210
         INFO      Now using 120 bins                                                          SpectrumLike.py:1673
         INFO      Range 9-900 translates to channels 5-123                                    SpectrumLike.py:1210
         INFO      Now using 119 bins                                                          SpectrumLike.py:1673
         INFO      Range 250-30000 translates to channels 1-119                                SpectrumLike.py:1210
         INFO      Now using 115 bins                                                          SpectrumLike.py:1673
         INFO      sampler set to multinest                                                bayesian_analysis.py:197
 *****************************************************
 MultiNest v3.10
 Copyright Farhan Feroz & Mike Hobson
 Release Jul 2015

 no. of live points =  500
 dimensionality =    4
 *****************************************************
 ln(ev)=  -907.65839243289236      +/-  0.19080663463205622
 Total Likelihood Evaluations:        21407
 Sampling finished. Exiting MultiNest
  analysing data from chains/fit-.txt
19:58:24 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.580 -0.014 +0.5) x 10^-2 1 / (cm2 keV s)
grb.spectrum.main.Band.alpha -1.04 -0.05 +0.13
grb.spectrum.main.Band.xp (5.61 -2.1 +0.04) x 10^2 keV
grb.spectrum.main.Band.beta -1.87 -0.06 +0.14
Values of -log(posterior) at the minimum:

-log(posterior)
b0_interval2 -324.320940
n3_interval2 -288.922222
n4_interval2 -312.309663
total -925.552825
Values of statistical measures:

statistical measures
AIC 1859.218964
BIC 1874.627781
DIC 1804.920868
PDIC 2.068446
log(Z) -394.191031
         INFO      Range 9-900 translates to channels 5-124                                    SpectrumLike.py:1210
         INFO      Now using 120 bins                                                          SpectrumLike.py:1673
         INFO      Range 9-900 translates to channels 5-123                                    SpectrumLike.py:1210
         INFO      Now using 119 bins                                                          SpectrumLike.py:1673
         INFO      Range 250-30000 translates to channels 1-119                                SpectrumLike.py:1210
         INFO      Now using 109 bins                                                          SpectrumLike.py:1673
         INFO      sampler set to multinest                                                bayesian_analysis.py:197
 *****************************************************
 MultiNest v3.10
 Copyright Farhan Feroz & Mike Hobson
 Release Jul 2015

 no. of live points =  500
 dimensionality =    4
 *****************************************************
 ln(ev)=  -788.85051640242341      +/-  0.18182012474359940
 Total Likelihood Evaluations:        15796
 Sampling finished. Exiting MultiNest
  analysing data from chains/fit-.txt
19:58:40 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.94 -0.4 +0.20) x 10^-2 1 / (cm2 keV s)
grb.spectrum.main.Band.alpha (-9.1 -1.3 +0.5) x 10^-1
grb.spectrum.main.Band.xp (3.35 -0.33 +1.0) x 10^2 keV
grb.spectrum.main.Band.beta -2.26 -0.6 -0.05
Values of -log(posterior) at the minimum:

-log(posterior)
b0_interval3 -298.318496
n3_interval3 -242.451951
n4_interval3 -262.645793
total -803.416239
Values of statistical measures:

statistical measures
AIC 1614.945793
BIC 1630.354610
DIC 1571.105556
PDIC 3.054340
log(Z) -342.593426
         INFO      Range 9-900 translates to channels 5-124                                    SpectrumLike.py:1210
         INFO      Now using 120 bins                                                          SpectrumLike.py:1673
         INFO      Range 9-900 translates to channels 5-123                                    SpectrumLike.py:1210
         INFO      Now using 119 bins                                                          SpectrumLike.py:1673
         INFO      Range 250-30000 translates to channels 1-119                                SpectrumLike.py:1210
         INFO      Now using 119 bins                                                          SpectrumLike.py:1673
         INFO      sampler set to multinest                                                bayesian_analysis.py:197
 *****************************************************
 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)=  -2271.2239357411418      +/-  0.19814931923289439
 Total Likelihood Evaluations:        21286
 Sampling finished. Exiting MultiNest

19:58:59 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.02 -0.04 +0.17) x 10^-2 1 / (cm2 keV s)
grb.spectrum.main.Band.alpha (-9.91 -0.18 +0.6) x 10^-1
grb.spectrum.main.Band.xp (4.1 -0.6 +0.4) x 10^2 keV
grb.spectrum.main.Band.beta -1.94 -0.22 +0.06
Values of -log(posterior) at the minimum:

-log(posterior)
b0_interval4 -778.629754
n3_interval4 -756.872764
n4_interval4 -747.087004
total -2282.589522
Values of statistical measures:

statistical measures
AIC 4573.292357
BIC 4588.701175
DIC 4528.835737
PDIC 3.696591
log(Z) -986.380022
         INFO      Range 9-900 translates to channels 5-124                                    SpectrumLike.py:1210
         INFO      Now using 120 bins                                                          SpectrumLike.py:1673
         INFO      Range 9-900 translates to channels 5-123                                    SpectrumLike.py:1210
         INFO      Now using 119 bins                                                          SpectrumLike.py:1673
         INFO      Range 250-30000 translates to channels 1-119                                SpectrumLike.py:1210
         INFO      Now using 119 bins                                                          SpectrumLike.py:1673
         INFO      sampler set to multinest                                                bayesian_analysis.py:197
 *****************************************************
 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)=  -1573.6747124769859      +/-  0.19549391648304051
 Total Likelihood Evaluations:        18383
 Sampling finished. Exiting MultiNest

19:59:17 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.79 -0.23 +0.13) x 10^-2 1 / (cm2 keV s)
grb.spectrum.main.Band.alpha (-9.1 -0.6 +0.4) x 10^-1
grb.spectrum.main.Band.xp (4.22 -0.35 +0.8) x 10^2 keV
grb.spectrum.main.Band.beta -2.10 -0.4 +0.05
Values of -log(posterior) at the minimum:

-log(posterior)
b0_interval5 -536.958539
n3_interval5 -523.525117
n4_interval5 -527.553004
total -1588.036660
Values of statistical measures:

statistical measures
AIC 3184.186635
BIC 3199.595452
DIC 3136.455522
PDIC 3.175818
log(Z) -683.438244
         INFO      Range 9-900 translates to channels 5-124                                    SpectrumLike.py:1210
         INFO      Now using 120 bins                                                          SpectrumLike.py:1673
         INFO      Range 9-900 translates to channels 5-123                                    SpectrumLike.py:1210
         INFO      Now using 119 bins                                                          SpectrumLike.py:1673
         INFO      Range 250-30000 translates to channels 1-119                                SpectrumLike.py:1210
         INFO      Now using 119 bins                                                          SpectrumLike.py:1673
         INFO      sampler set to multinest                                                bayesian_analysis.py:197
 *****************************************************
 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)=  -1756.1459273624052      +/-  0.19546968038357976
 Total Likelihood Evaluations:        19565
 Sampling finished. Exiting MultiNest

19:59: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 (1.95 -0.04 +0.15) x 10^-2 1 / (cm2 keV s)
grb.spectrum.main.Band.alpha -1.012 -0.015 +0.05
grb.spectrum.main.Band.xp (4.45 -0.6 +0.20) x 10^2 keV
grb.spectrum.main.Band.beta -2.29 -0.21 +0.17
Values of -log(posterior) at the minimum:

-log(posterior)
b0_interval6 -609.530873
n3_interval6 -584.191725
n4_interval6 -576.735082
total -1770.457679
Values of statistical measures:

statistical measures
AIC 3549.028673
BIC 3564.437490
DIC 3499.844409
PDIC 2.735885
log(Z) -762.684486
         INFO      Range 9-900 translates to channels 5-124                                    SpectrumLike.py:1210
         INFO      Now using 120 bins                                                          SpectrumLike.py:1673
         INFO      Range 9-900 translates to channels 5-123                                    SpectrumLike.py:1210
         INFO      Now using 119 bins                                                          SpectrumLike.py:1673
         INFO      Range 250-30000 translates to channels 1-119                                SpectrumLike.py:1210
         INFO      Now using 119 bins                                                          SpectrumLike.py:1673
         INFO      sampler set to multinest                                                bayesian_analysis.py:197
 *****************************************************
 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)=  -1940.8268798664146      +/-  0.19535807852060383
 Total Likelihood Evaluations:        20109
 Sampling finished. Exiting MultiNest

19:59: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 (1.73 -0.05 +0.14) x 10^-2 1 / (cm2 keV s)
grb.spectrum.main.Band.alpha -1.02 -0.04 +0.06
grb.spectrum.main.Band.xp (4.2 -0.7 +0.4) x 10^2 keV
grb.spectrum.main.Band.beta -2.52 -0.17 +0.14
Values of -log(posterior) at the minimum:

-log(posterior)
b0_interval7 -662.253823
n3_interval7 -642.286476
n4_interval7 -650.391592
total -1954.931892
Values of statistical measures:

statistical measures
AIC 3917.977098
BIC 3933.385915
DIC 3869.347472
PDIC 2.302357
log(Z) -842.890404
         INFO      Range 9-900 translates to channels 5-124                                    SpectrumLike.py:1210
         INFO      Now using 120 bins                                                          SpectrumLike.py:1673
         INFO      Range 9-900 translates to channels 5-123                                    SpectrumLike.py:1210
         INFO      Now using 119 bins                                                          SpectrumLike.py:1673
         INFO      Range 250-30000 translates to channels 1-119                                SpectrumLike.py:1210
         INFO      Now using 119 bins                                                          SpectrumLike.py:1673
         INFO      sampler set to multinest                                                bayesian_analysis.py:197
 *****************************************************
 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)=  -2056.0178913742088      +/-  0.19024138082103711
 Total Likelihood Evaluations:        19504
 Sampling finished. Exiting MultiNest

20:00: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 (1.59 +0.04 +0.23) x 10^-2 1 / (cm2 keV s)
grb.spectrum.main.Band.alpha (-8.22 +0.15 +1.0) x 10^-1
grb.spectrum.main.Band.xp (3.50 -0.6 -0.14) x 10^2 keV
grb.spectrum.main.Band.beta -2.13 -0.22 +0.10
Values of -log(posterior) at the minimum:

-log(posterior)
b0_interval8 -702.382627
n3_interval8 -698.406012
n4_interval8 -665.520648
total -2066.309287
Values of statistical measures:

statistical measures
AIC 4140.731889
BIC 4156.140706
DIC 4098.611711
PDIC 2.954964
log(Z) -892.917225
         INFO      Range 9-900 translates to channels 5-124                                    SpectrumLike.py:1210
         INFO      Now using 120 bins                                                          SpectrumLike.py:1673
         INFO      Range 9-900 translates to channels 5-123                                    SpectrumLike.py:1210
         INFO      Now using 119 bins                                                          SpectrumLike.py:1673
         INFO      Range 250-30000 translates to channels 1-119                                SpectrumLike.py:1210
         INFO      Now using 119 bins                                                          SpectrumLike.py:1673
         INFO      sampler set to multinest                                                bayesian_analysis.py:197
 *****************************************************
 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)=  -1879.1886894987244      +/-  0.14913583245844605
 Total Likelihood Evaluations:        12511
 Sampling finished. Exiting MultiNest

20:00:24 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.31 +0.7) x 10^-2 1 / (cm2 keV s)
grb.spectrum.main.Band.alpha (-8.8 -1.7 +2.3) x 10^-1
grb.spectrum.main.Band.xp (1.17 -0.35 +0.5) x 10^2 keV
grb.spectrum.main.Band.beta -1.88 -0.34 +0.10
Values of -log(posterior) at the minimum:

-log(posterior)
b0_interval9 -648.461109
n3_interval9 -617.109986
n4_interval9 -616.239811
total -1881.810905
Values of statistical measures:

statistical measures
AIC 3771.735125
BIC 3787.143942
DIC 3724.385773
PDIC -22.501877
log(Z) -816.121278
         INFO      Range 9-900 translates to channels 5-124                                    SpectrumLike.py:1210
         INFO      Now using 120 bins                                                          SpectrumLike.py:1673
         INFO      Range 9-900 translates to channels 5-123                                    SpectrumLike.py:1210
         INFO      Now using 119 bins                                                          SpectrumLike.py:1673
         INFO      Range 250-30000 translates to channels 1-119                                SpectrumLike.py:1210
         INFO      Now using 119 bins                                                          SpectrumLike.py:1673
         INFO      sampler set to multinest                                                bayesian_analysis.py:197
 *****************************************************
 MultiNest v3.10
 Copyright Farhan Feroz & Mike Hobson
 Release Jul 2015

 no. of live points =  500
 dimensionality =    4
 *****************************************************
 ln(ev)=  -1322.3359560867518      +/-  0.16742218527554686
 Total Likelihood Evaluations:        15089
 Sampling finished. Exiting MultiNest
  analysing data from chains/fit-.txt
20:00: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
grb.spectrum.main.Band.K (2.1 -0.4 +0.6) x 10^-2 1 / (cm2 keV s)
grb.spectrum.main.Band.alpha (-7.2 -1.4 +1.6) x 10^-1
grb.spectrum.main.Band.xp (2.1 -0.4 +0.5) x 10^2 keV
grb.spectrum.main.Band.beta -1.91 -0.34 +0.07
Values of -log(posterior) at the minimum:

-log(posterior)
b0_interval10 -460.936580
n3_interval10 -437.554291
n4_interval10 -433.094760
total -1331.585631
Values of statistical measures:

statistical measures
AIC 2671.284577
BIC 2686.693394
DIC 2634.042059
PDIC 0.606781
log(Z) -574.283209
         INFO      Range 9-900 translates to channels 5-124                                    SpectrumLike.py:1210
         INFO      Now using 120 bins                                                          SpectrumLike.py:1673
         INFO      Range 9-900 translates to channels 5-123                                    SpectrumLike.py:1210
         INFO      Now using 119 bins                                                          SpectrumLike.py:1673
         INFO      Range 250-30000 translates to channels 1-119                                SpectrumLike.py:1210
         INFO      Now using 119 bins                                                          SpectrumLike.py:1673
         INFO      sampler set to multinest                                                bayesian_analysis.py:197
 *****************************************************
 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.43776871464650      +/-  0.14334945303256250
 Total Likelihood Evaluations:        12828
 Sampling finished. Exiting MultiNest

20:00:51 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.7 +2.1) x 10^-2 1 / (cm2 keV s)
grb.spectrum.main.Band.alpha (-5.0 -1.8 +3.5) x 10^-1
grb.spectrum.main.Band.xp (1.31 -0.31 +0.28) x 10^2 keV
grb.spectrum.main.Band.beta -2.20 -0.4 +0.29
Values of -log(posterior) at the minimum:

-log(posterior)
b0_interval11 -292.463715
n3_interval11 -272.512822
n4_interval11 -255.820360
total -820.796897
Values of statistical measures:

statistical measures
AIC 1649.707109
BIC 1665.115927
DIC 1617.415968
PDIC 0.016019
log(Z) -352.402945

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
20:00:52 INFO      fit restored to median of posterior                                          sampler_base.py:164
20:00:53 INFO      fit restored to median of posterior                                          sampler_base.py:164
20:00:54 INFO      fit restored to median of posterior                                          sampler_base.py:164
         INFO      fit restored to median of posterior                                          sampler_base.py:164
20:00:55 INFO      fit restored to median of posterior                                          sampler_base.py:164
20:00:56 INFO      fit restored to median of posterior                                          sampler_base.py:164
20:00:57 INFO      fit restored to median of posterior                                          sampler_base.py:164
20:00:58 INFO      fit restored to median of posterior                                          sampler_base.py:164
20:00:59 INFO      fit restored to median of posterior                                          sampler_base.py:164
20:01:00 INFO      fit restored to median of posterior                                          sampler_base.py:164
         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.