Quickstart¶

A short walk-through of fitting a gamma-ray spectrum with BaySpec. It covers the three stages of a typical session:

Data — joint spectra from Fermi/GBM’s NaI and BGO detectors.
Model — a cutoff power-law (cpl).
Inference — Bayesian posterior sampling via multinest (emcee works as a drop-in replacement).

import numpy as np
from bayspec.model.local import *
from bayspec import DataUnit, Data, BayesInfer, Plot, MaxLikeFit

savepath = 'quickstart'

Load spectra and response data.

nai = DataUnit(
    src='./ME/me.src',
    bkg='./ME/me.bkg',
    rsp='./ME/me.rsp',
    notc=[8, 900],
    stat='pgstat',
    rebn={'min_sigma': 2, 'max_bin': 10})

bgo = DataUnit(
    src='./HE/he.src',
    bkg='./HE/he.bkg',
    rsp='./HE/he.rsp',
    notc=[300, 38000],
    stat='pgstat',
    rebn={'min_sigma': 2, 'max_bin': 10})

data = Data([('nai', nai),
             ('bgo', bgo)])

data.save(savepath)
data

Data

Name	Noticing	Statistic	Grouping	Time
nai	[[8, 900]]	pgstat	None	None
bgo	[[300, 38000]]	pgstat	None	None

Data Parameters

par#	Component	Parameter	Value	Prior	Frozen
1	nai	sf	1	None	True
2	nai	bf	1	None	True
3	nai	rf	1	None	True
4	nai	ra	0	None	True
5	nai	dec	0	None	True
6	bgo	sf	1	None	True
7	bgo	bf	1	None	True
8	bgo	rf	1	None	True
9	bgo	ra	0	None	True
10	bgo	dec	0	None	True

Define the spectral model.

model = cpl()
model.save(savepath)
model

Model

cpl [add]

power-law model with high-energy cutoff

Model Configurations

cfg#	Component	Parameter	Value
1	cpl	redshift	0.000
2	cpl	pivot_energy	1.000
3	cpl	vfv_peak	True

Model Parameters

par#	Component	Parameter	Value	Prior
1	cpl	$\alpha$	-1	unif(-2, 2)
2	cpl	log$E_p$	2	unif(0, 4)
3	cpl	log$A$	0	unif(-10, 10)

Run Bayesian inference.

infer = BayesInfer([(data, model)])
infer.save(savepath)
infer

Bayesian Inference

Configurations

cfg#	Class	Expression	Component	Parameter	Value
1	model	cpl	cpl	redshift	0.000
2	model	cpl	cpl	pivot_energy	1.000
3	model	cpl	cpl	vfv_peak	True

Parameters

par#	Class	Expression	Component	Parameter	Value	Prior
1*	model	cpl	cpl	$\alpha$	-1	unif(-2, 2)
2*	model	cpl	cpl	log$E_p$	2	unif(0, 4)
3*	model	cpl	cpl	log$A$	0	unif(-10, 10)

post = infer.multinest(nlive=400, resume=True, verbose=False, savepath=savepath)
post.save(savepath)
post

Posterior Results

Parameters

par#	Class	Expression	Component	Parameter	Mean	Median	Best	1sigma Best	1sigma CI
1	model	cpl	cpl	$\alpha$	-1.562	-1.562	-1.561	-1.561	[-1.572, -1.551]
2	model	cpl	cpl	log$E_p$	2.331	2.331	2.330	2.330	[2.321, 2.341]
3	model	cpl	cpl	log$A$	2.353	2.354	2.352	2.352	[2.338, 2.368]

Statistics

Data	Model	Statistic	Value	Bins
nai	cpl	pgstat	405.171	119
bgo	cpl	pgstat	149.560	117
Total	Total	stat/dof	554.731/233	236

Information Criterias

AIC	AICc	BIC	lnZ
560.731	560.834	571.123	-295.844

fig = Plot.infer(post, style='CE')
fig.save(f'{savepath}/ctsspec')

fig = Plot.post_corner(post, ploter='plotly')
fig.save(f'{savepath}/corner')

earr = np.logspace(1, 3, 100)

modelplot = Plot.model(ploter='plotly', style='vFv', post=True)
modelplot.add_model(model, E=earr)
fig = modelplot.get_fig()
fig.save(f'{savepath}/model')

ergflux = model.best_ergflux(emin=10, emax=1000, ngrid=1000)
ergflux_sample = model.ergflux_sample(emin=10, emax=1000, ngrid=1000)

print(ergflux)
print(ergflux_sample)

8.367943905909297e-06
{'mean': 8.369501130327865e-06, 'median': 8.369420051840691e-06, 'Isigma': array([8.31271099e-06, 8.42362544e-06]), 'IIsigma': array([8.25891689e-06, 8.48377819e-06]), 'IIIsigma': array([8.20528996e-06, 8.54385474e-06])}