Analyzing GRB 080916C
(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")
22:57:14 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:
Free parameters (5):
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)
| 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)
22:58:11 INFO Auto-determined polynomial order: 0 binned_spectrum_series.py:389
22:58:20 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
22:58:21 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
22:58:23 INFO Now using 120 bins SpectrumLike.py:1739
22:58:25 INFO Auto-determined polynomial order: 1 binned_spectrum_series.py:389
22:58:35 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
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
22:58:36 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
22:58:38 INFO Auto-determined polynomial order: 1 binned_spectrum_series.py:389
22:58:49 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
22:58:50 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: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
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)
22:58:51 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)
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.7881721717995 +/- 0.22134783228260471
Total Likelihood Evaluations: 27495
Sampling finished. Exiting MultiNest
22:59:11 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.467 -0.018 +0.021) x 10^-2 | 1 / (keV s cm2) |
| GRB080916009...alpha | -1.0921 -0.0008 +0.04 | |
| GRB080916009...break_energy | (1.88 +0.15 +0.8) x 10^2 | keV |
| GRB080916009...break_scale | (0.0 +1.8 +3.5) x 10^-1 | |
| GRB080916009...beta | -1.96 -0.26 -0.05 |
Values of -log(posterior) at the minimum:
| -log(posterior) | |
|---|---|
| n3 | -1019.832436 |
| n4 | -1010.904188 |
| b0 | -1050.827196 |
| total | -3081.563820 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 6173.298095 |
| BIC | 6192.530306 |
| DIC | 6180.668774 |
| PDIC | 4.612109 |
| log(Z) | -1346.655193 |
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
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")
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)
23:01:17 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)
[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")
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)
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
23:01:18 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
23:01:19 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.65905072157466 +/- 0.17841336460984905
Total Likelihood Evaluations: 15787
Sampling finished. Exiting MultiNest
23:01:31 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.5 -0.5 +0.7) x 10^-2 | 1 / (keV s cm2) |
| grb.spectrum.main.Band.alpha | (-5.8 -1.0 +1.5) x 10^-1 | |
| grb.spectrum.main.Band.xp | (3.2 -0.5 +0.7) x 10^2 | keV |
| grb.spectrum.main.Band.beta | -2.04 -0.4 +0.11 |
Values of -log(posterior) at the minimum:
| -log(posterior) | |
|---|---|
| n3_interval0 | -250.261638 |
| n4_interval0 | -268.005932 |
| b0_interval0 | -285.703224 |
| total | -803.970794 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 1616.054903 |
| BIC | 1631.463720 |
| DIC | 1570.380534 |
| PDIC | 2.502667 |
| log(Z) | -342.510274 |
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)= -1974.3529309397813 +/- 0.23619326152835496
Total Likelihood Evaluations: 22922
Sampling finished. Exiting MultiNest
23:01:48 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.091 -0.012 +0.06) x 10^-2 | 1 / (keV s cm2) |
| grb.spectrum.main.Band.alpha | (-7.4500 -0.06 +0.0031) x 10^-1 | |
| grb.spectrum.main.Band.xp | (3.8052 -0.12 -0.0013) x 10^2 | keV |
| grb.spectrum.main.Band.beta | -1.7923 +0.0020 +0.016 |
Values of -log(posterior) at the minimum:
| -log(posterior) | |
|---|---|
| n3_interval1 | -645.514550 |
| n4_interval1 | -643.489494 |
| b0_interval1 | -690.408206 |
| total | -1979.412250 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 3966.937814 |
| BIC | 3982.346632 |
| DIC | 3917.516285 |
| PDIC | 0.989985 |
| log(Z) | -857.450583 |
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)= -908.95526290612497 +/- 0.19918528323311785
Total Likelihood Evaluations: 20413
Sampling finished. Exiting MultiNest
23:02: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.84 -0.4 +0.14) x 10^-2 | 1 / (keV s cm2) |
| grb.spectrum.main.Band.alpha | (-9.92 -1.1 +0.31) x 10^-1 | |
| grb.spectrum.main.Band.xp | (4.2 -0.5 +2.8) x 10^2 | keV |
| grb.spectrum.main.Band.beta | -1.81 -0.5 +0.05 |
Values of -log(posterior) at the minimum:
| -log(posterior) | |
|---|---|
| n3_interval2 | -288.465710 |
| n4_interval2 | -311.286092 |
| b0_interval2 | -323.699338 |
| total | -923.451140 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 1855.015594 |
| BIC | 1870.424412 |
| DIC | 1805.484287 |
| PDIC | 1.650166 |
| log(Z) | -394.754255 |
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)= -792.48946267470023 +/- 0.18510546944558171
Total Likelihood Evaluations: 17692
Sampling finished. Exiting MultiNest
23:02: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 | (2.95 -0.17 +0.6) x 10^-2 | 1 / (keV s cm2) |
| grb.spectrum.main.Band.alpha | (-9.4 -0.4 +1.0) x 10^-1 | |
| grb.spectrum.main.Band.xp | (2.96 -0.7 +0.34) x 10^2 | keV |
| grb.spectrum.main.Band.beta | -1.91 +0.04 +0.14 |
Values of -log(posterior) at the minimum:
| -log(posterior) | |
|---|---|
| n3_interval3 | -242.254457 |
| n4_interval3 | -261.851411 |
| b0_interval3 | -298.847099 |
| total | -802.952967 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 1614.019248 |
| BIC | 1629.428066 |
| DIC | 1571.725700 |
| PDIC | 1.654590 |
| log(Z) | -344.173801 |
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
23:02:14 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)= -2271.4013348635499 +/- 0.20239646827124147
Total Likelihood Evaluations: 19977
Sampling finished. Exiting MultiNest
23:02:26 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.07 -0.08 +0.13) x 10^-2 | 1 / (keV s cm2) |
| grb.spectrum.main.Band.alpha | (-9.77 -0.24 +0.5) x 10^-1 | |
| grb.spectrum.main.Band.xp | (3.9 +/- 0.4) x 10^2 | keV |
| grb.spectrum.main.Band.beta | -1.96 +/- 0.06 |
Values of -log(posterior) at the minimum:
| -log(posterior) | |
|---|---|
| n3_interval4 | -757.198020 |
| n4_interval4 | -746.558302 |
| b0_interval4 | -778.429692 |
| total | -2282.186014 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 4572.485343 |
| BIC | 4587.894161 |
| DIC | 4526.456549 |
| PDIC | 2.793861 |
| log(Z) | -986.457066 |
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)= -1572.8379336512810 +/- 0.19009096217207097
Total Likelihood Evaluations: 20872
Sampling finished. Exiting MultiNest
23:02: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.82 -0.24 +0.12) x 10^-2 | 1 / (keV s cm2) |
| grb.spectrum.main.Band.alpha | (-9.1 -0.6 +0.4) x 10^-1 | |
| grb.spectrum.main.Band.xp | (4.16 -0.34 +0.9) x 10^2 | keV |
| grb.spectrum.main.Band.beta | -2.13 -0.31 +0.09 |
Values of -log(posterior) at the minimum:
| -log(posterior) | |
|---|---|
| n3_interval5 | -523.552781 |
| n4_interval5 | -527.678010 |
| b0_interval5 | -536.696935 |
| total | -1587.927726 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 3183.968767 |
| BIC | 3199.377584 |
| DIC | 3137.012038 |
| PDIC | 3.433953 |
| log(Z) | -683.074836 |
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)= -1779.7947536540012 +/- 0.20592974147790152
Total Likelihood Evaluations: 19592
Sampling finished. Exiting MultiNest
23:02:52 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.05 +1.7) x 10^-2 | 1 / (keV s cm2) |
| grb.spectrum.main.Band.alpha | -1.12 -0.05 +0.5 | |
| grb.spectrum.main.Band.xp | (3.4 -1.7 +0.7) x 10^2 | keV |
| grb.spectrum.main.Band.beta | -1.89 -0.06 +0.18 |
Values of -log(posterior) at the minimum:
| -log(posterior) | |
|---|---|
| n3_interval6 | -593.150150 |
| n4_interval6 | -581.642232 |
| b0_interval6 | -610.517788 |
| total | -1785.310171 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 3578.733656 |
| BIC | 3594.142473 |
| DIC | 3495.722983 |
| PDIC | -40.777663 |
| log(Z) | -772.955040 |
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)= -1938.9720369688273 +/- 0.18982528363775455
Total Likelihood Evaluations: 21037
Sampling finished. Exiting MultiNest
23:03: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 | (1.66 -0.11 +0.10) x 10^-2 | 1 / (keV s cm2) |
| grb.spectrum.main.Band.alpha | -1.04 -0.05 +0.04 | |
| grb.spectrum.main.Band.xp | (4.3 -0.5 +0.9) x 10^2 | keV |
| grb.spectrum.main.Band.beta | -2.27 -0.4 +0.14 |
Values of -log(posterior) at the minimum:
| -log(posterior) | |
|---|---|
| n3_interval7 | -640.867944 |
| n4_interval7 | -650.481960 |
| b0_interval7 | -662.309470 |
| total | -1953.659375 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 3915.432064 |
| BIC | 3930.840881 |
| DIC | 3869.346733 |
| PDIC | 3.533885 |
| log(Z) | -842.084856 |
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)= -2054.3764041470190 +/- 0.18884837770541682
Total Likelihood Evaluations: 19282
Sampling finished. Exiting MultiNest
23:03: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 | (1.52 -0.08 +0.16) x 10^-2 | 1 / (keV s cm2) |
| grb.spectrum.main.Band.alpha | (-8.5 -0.5 +0.8) x 10^-1 | |
| grb.spectrum.main.Band.xp | (3.8 -0.6 +0.4) x 10^2 | keV |
| grb.spectrum.main.Band.beta | -2.32 -0.18 +0.21 |
Values of -log(posterior) at the minimum:
| -log(posterior) | |
|---|---|
| n3_interval8 | -698.434599 |
| n4_interval8 | -666.207579 |
| b0_interval8 | -702.375280 |
| total | -2067.017458 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 4142.148231 |
| BIC | 4157.557048 |
| DIC | 4097.442549 |
| PDIC | 3.045416 |
| log(Z) | -892.204336 |
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.4023371054529 +/- 0.14302829647889262
Total Likelihood Evaluations: 13388
Sampling finished. Exiting MultiNest
23:03:25 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.29 +1.4) x 10^-2 | 1 / (keV s cm2) |
| grb.spectrum.main.Band.alpha | (-8.7 -1.5 +4) x 10^-1 | |
| grb.spectrum.main.Band.xp | (1.2 -0.5 +0.4) x 10^2 | keV |
| grb.spectrum.main.Band.beta | -1.87 -0.34 +0.15 |
Values of -log(posterior) at the minimum:
| -log(posterior) | |
|---|---|
| n3_interval9 | -617.095666 |
| n4_interval9 | -616.260893 |
| b0_interval9 | -648.547760 |
| total | -1881.904318 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 3771.921951 |
| BIC | 3787.330768 |
| DIC | 3692.886255 |
| PDIC | -54.219346 |
| log(Z) | -815.779770 |
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.3251398048114 +/- 0.16899122660723453
Total Likelihood Evaluations: 15304
Sampling finished. Exiting MultiNest
23:03: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.00 -0.35 +0.5) x 10^-2 | 1 / (keV s cm2) |
| grb.spectrum.main.Band.alpha | (-7.5 +/- 1.3) x 10^-1 | |
| grb.spectrum.main.Band.xp | (2.2 -0.4 +0.6) x 10^2 | keV |
| grb.spectrum.main.Band.beta | -1.95 -0.25 +0.07 |
Values of -log(posterior) at the minimum:
| -log(posterior) | |
|---|---|
| n3_interval10 | -437.558558 |
| n4_interval10 | -433.380356 |
| b0_interval10 | -460.928265 |
| total | -1331.867179 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 2671.847672 |
| BIC | 2687.256489 |
| DIC | 2633.508541 |
| PDIC | 0.291402 |
| log(Z) | -574.278511 |
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)= -812.01821106508737 +/- 0.14761952222093477
Total Likelihood Evaluations: 12362
Sampling finished. Exiting MultiNest
23:03:45 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.8 +1.4) x 10^-2 | 1 / (keV s cm2) |
| grb.spectrum.main.Band.alpha | (-5.0 -2.0 +2.5) x 10^-1 | |
| grb.spectrum.main.Band.xp | (1.30 -0.23 +0.33) x 10^2 | keV |
| grb.spectrum.main.Band.beta | -2.12 -0.5 +0.13 |
Values of -log(posterior) at the minimum:
| -log(posterior) | |
|---|---|
| n3_interval11 | -272.322554 |
| n4_interval11 | -255.869941 |
| b0_interval11 | -292.337228 |
| total | -820.529723 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 1649.172759 |
| BIC | 1664.581577 |
| DIC | 1618.299121 |
| PDIC | 1.221754 |
| log(Z) | -352.655028 |
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
23:03:46 INFO fit restored to median of posterior sampler_base.py:164
23:03:47 INFO fit restored to median of posterior sampler_base.py:164
23:03:48 INFO fit restored to median of posterior sampler_base.py:164
23:03:49 INFO fit restored to median of posterior sampler_base.py:164
INFO fit restored to median of posterior sampler_base.py:164
23:03:50 INFO fit restored to median of posterior sampler_base.py:164
23:03:52 INFO fit restored to median of posterior sampler_base.py:164
23:03:53 INFO fit restored to median of posterior sampler_base.py:164
23:03:54 INFO fit restored to median of posterior sampler_base.py:164
23:03:55 INFO fit restored to median of posterior sampler_base.py:164
23:03:57 INFO fit restored to median of posterior sampler_base.py:164
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),
)
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.