Background Modeling

When fitting a spectrum with a background, it is invalid to simply subtract off the background if the background is part of the data’s generative model van Dyk et al. (2001). Therefore, we are often left with the task of modeling the statistical process of the background along with our source.

In typical spectral modeling, we find a few common cases when background is involved. If we have total counts ($$S_i$$) in $$i^{\rm th}$$ on $$N$$ bins observed for an exposure of $$t_{\rm s}$$ and also a measurement of $$B_i$$ background counts from looking off source for $$t_{\rm b}$$ seconds, we can then suppose a model for the source rate ($$m_i$$) and background rate ($$b_i$$).

Poisson source with Poisson background

This is described by a likelihood of the following form:

$L = \prod^N_{i=1} \frac{(t_{\rm s}(m_i+b_i))^{S_i} e^{-t_{\rm s}(m_i+b_i)}}{S_i!} \times \frac{(t_{\rm b} b_i)^{B_i} e^{-t_{\rm b}b_i}}{B_i!}$

which is a Poisson likelihood for the total model ($$m_i +b_i$$) conditional on the Poisson distributed background observation. This is the typical case for e.g. aperture x-ray instruments that observe a source region and then a background region. Both observations are Poisson distributed.

Poisson source with Gaussian background

This likelihood is similar, but the conditonal background distribution is described by Gaussian:

$L = \prod^N_{i=1} \frac{(t_{\rm s}(m_i+b_i))^{S_i} e^{-t_{\rm s}(m_i+b_i)}}{S_i!} \times \frac{1}{\sigma_{b,i}\sqrt{2 \pi}} \exp \left[ \frac{({B_i} - t_{\rm b} b_i)^2} {2 \sigma_{b,i}^2} \right]$

where the $$\sigma_{b,i}$$ are the measured errors on $$B_i$$. This situation occurs e.g. when the background counts are estimated from a fitted model such as time-domain instruments that estimate the background counts from temporal fits to the lightcurve.

In 3ML, we can fit a background model along with the the source model which allows for arbitrarily low background counts (in fact zero) in channels. The alternative is to use profile likelihoods where we first differentiate the likelihood with respect to the background model

$\frac{ \partial L}{{\partial b_i}} = 0$

and solve for the $$b_i$$ that maximize the likelihood. Both the Poisson and Gaussian background profile likelihoods are described in the XSPEC statistics guide. This implicitly yields $$N$$ parameters to the model thus requiring at least one background count per channel. These profile likelihoods are the default Poisson likelihoods in 3ML when a background model is not used with a SpectrumLike (and its children, DispersionSpectrumLike and OGIPLike) plugin.

Let’s examine how to handle both cases.

[1]:

import warnings

warnings.simplefilter("ignore")
import numpy as np

np.seterr(all="ignore")

[1]:

{'divide': 'warn', 'over': 'warn', 'under': 'ignore', 'invalid': 'warn'}

[2]:

%%capture
from threeML import *

[3]:

from jupyterthemes import jtplot

%matplotlib inline
jtplot.style(context="talk", fscale=1, ticks=True, grid=False)
set_threeML_style()
silence_warnings()
import astropy.units as u


First we will create an observation where we have a simulated broken power law source spectrum along with an observed background spectrum. The background is a powerl law continuum with a Gaussian line.

[4]:

# create the simulated observation

energies = np.logspace(1, 4, 151)

low_edge = energies[:-1]
high_edge = energies[1:]

# get a BPL source function
source_function = Broken_powerlaw(K=2, xb=300, piv=300, alpha=0.0, beta=-3.0)

# power law background function
background_function = Powerlaw(K=0.5, index=-1.5, piv=100.0) + Gaussian(
F=50, mu=511, sigma=20
)

spectrum_generator = SpectrumLike.from_function(
"fake",
source_function=source_function,
background_function=background_function,
energy_min=low_edge,
energy_max=high_edge,
)

spectrum_generator.view_count_spectrum()

08:30:38 INFO      Auto-probed noise models:                                                    SpectrumLike.py:491

         INFO      - observation: poisson                                                       SpectrumLike.py:492

         INFO      - background: None                                                           SpectrumLike.py:493

08:30:40 INFO      Auto-probed noise models:                                                    SpectrumLike.py:491

         INFO      - observation: poisson                                                       SpectrumLike.py:492

         INFO      - background: None                                                           SpectrumLike.py:493

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

         INFO      - observation: poisson                                                       SpectrumLike.py:492

         INFO      - background: poisson                                                        SpectrumLike.py:493

08:30:41 INFO      Auto-probed noise models:                                                    SpectrumLike.py:491

         INFO      - observation: poisson                                                       SpectrumLike.py:492

         INFO      - background: poisson                                                        SpectrumLike.py:493

findfont: Font family ['sans-serif'] not found. Falling back to DejaVu Sans.
findfont: Generic family 'sans-serif' not found because none of the following families were found: Helvetica
findfont: Generic family 'sans-serif' not found because none of the following families were found: Helvetica
findfont: Generic family 'sans-serif' not found because none of the following families were found: Helvetica

[4]:


Using a profile likelihood

We have very few counts counts in some channels (in fact sometimes zero), but let’s assume we do not know the model for the background. In this case, we will use the profile Poisson likelihood.

[5]:

# instance our source spectrum
bpl = Broken_powerlaw(piv=300, xb=500)

# instance a point source
ra, dec = 0, 0
ps_src = PointSource("source", ra, dec, spectral_shape=bpl)

# instance the likelihood model
src_model = Model(ps_src)

# pass everything to a joint likelihood object
jl_profile = JointLikelihood(src_model, DataList(spectrum_generator))

# fit the model
_ = jl_profile.fit()

# plot the fit in count space
_ = spectrum_generator.display_model(step=False)

08:30:43 INFO      set the minimizer to minuit                                             joint_likelihood.py:1042

Best fit values:


result unit
parameter
source.spectrum.main.Broken_powerlaw.K 1.85 +/- 0.13 1 / (cm2 keV s)
source.spectrum.main.Broken_powerlaw.xb (3.08 +/- 0.13) x 10^2 keV
source.spectrum.main.Broken_powerlaw.alpha (-8 +/- 7) x 10^-2
source.spectrum.main.Broken_powerlaw.beta -3.12 +/- 0.21

Correlation matrix:


 1 -0.59 0.76 0.02 -0.59 1 -0.45 -0.6 0.76 -0.45 1 0.01 0.02 -0.6 0.01 1

Values of -log(likelihood) at the minimum:


-log(likelihood)
fake 423.44375
total 423.44375

Values of statistical measures:


statistical measures
AIC 855.163362
BIC 866.930041
findfont: Font family ['sans-serif'] not found. Falling back to DejaVu Sans.
findfont: Generic family 'sans-serif' not found because none of the following families were found: Helvetica


Our fit recovers the simulated parameters. However, we should have binned the spectrum up such that there is at least one background count per spectral bin for the profile to be valid.

[6]:

spectrum_generator.rebin_on_background(1)

spectrum_generator.view_count_spectrum()

_ = jl_profile.fit()

_ = spectrum_generator.display_model(step=False)

08:30:47 INFO      Now using 87 bins                                                           SpectrumLike.py:1735

Best fit values:


result unit
parameter
source.spectrum.main.Broken_powerlaw.K 1.93 +/- 0.15 1 / (cm2 keV s)
source.spectrum.main.Broken_powerlaw.xb (2.98 +/- 0.15) x 10^2 keV
source.spectrum.main.Broken_powerlaw.alpha (-1 +/- 8) x 10^-2
source.spectrum.main.Broken_powerlaw.beta -3.03 +/- 0.20

Correlation matrix:


 1 -0.66 0.77 0.13 -0.66 1 -0.5 -0.67 0.77 -0.5 1 0.09 0.13 -0.67 0.09 1

Values of -log(likelihood) at the minimum:


-log(likelihood)
fake 305.473654
total 305.473654

Values of statistical measures:


statistical measures
AIC 619.223171
BIC 630.989850

Modeling the background

Now let’s try to model the background assuming we know that the background is a power law with a Gaussian line. We can extract a background plugin from the data by passing the original plugin to a classmethod of spectrum like.

[7]:

# extract the background from the spectrum plugin.
# This works for OGIPLike plugins as well, though we could easily also just read
# in a bakcground PHA
background_plugin = SpectrumLike.from_background("bkg", spectrum_generator)

08:30:50 INFO      Auto-probed noise models:                                                    SpectrumLike.py:491

         INFO      - observation: poisson                                                       SpectrumLike.py:492

         INFO      - background: None                                                           SpectrumLike.py:493


This constructs a new plugin with only the observed background so that we can first model it.

[8]:

background_plugin.view_count_spectrum()

[8]:


We now construct our background model and fit it to the data. Let’s assume we know that the line occurs at 511 keV, but we are unsure of its strength an width. We do not need to bin the data up because we are using a simple Poisson likelihood which is valid even when we have zero counts Cash (1979).

[9]:

# instance the spectrum setting the line's location to 511
bkg_spectrum = Powerlaw(piv=100) + Gaussian(F=50, mu=511)

# setup model parameters
# fix the line's location
bkg_spectrum.mu_2.fix = True

# nice parameter bounds
bkg_spectrum.K_1.bounds = (1e-4, 10)
bkg_spectrum.F_2.bounds = (0.0, 1000)
bkg_spectrum.sigma_2.bounds = (2, 30)

ps_bkg = PointSource("bkg", 0, 0, spectral_shape=bkg_spectrum)

bkg_model = Model(ps_bkg)

jl_bkg = JointLikelihood(bkg_model, DataList(background_plugin))

_ = jl_bkg.fit()

_ = background_plugin.display_model(
step=False, data_color="#1A68F0", model_color="#FF9700"
)

08:30:51 INFO      set the minimizer to minuit                                             joint_likelihood.py:1042

Best fit values:


result unit
parameter
bkg.spectrum.main.composite.K_1 (3.20 -0.23 +0.25) x 10^-1 1 / (cm2 keV s)
bkg.spectrum.main.composite.index_1 -1.35 +/- 0.04
bkg.spectrum.main.composite.F_2 (2.5 +/- 0.5) x 10 1 / (cm2 s)
bkg.spectrum.main.composite.sigma_2 (1.98 +/- 0.33) x 10 keV

Correlation matrix:


 1 0.08 -0.05 -0.01 0.08 1 -0.05 -0.01 -0.05 -0.05 1 0.03 -0.01 -0.01 0.03 1

Values of -log(likelihood) at the minimum:


-log(likelihood)
bkg 223.157239
total 223.157239

Values of statistical measures:


statistical measures
AIC 454.590340
BIC 466.357019

We now have a model and estimate for the background which we can use when fitting with the source spectrum. We now create a new plugin with just the total observation and pass our background plugin as the background argument.

[10]:

modeled_background_plugin = SpectrumLike(
"full",
# here we use the original observation
observation=spectrum_generator.observed_spectrum,
# we pass the background plugin as the background!
background=background_plugin,
)

08:30:52 INFO      Background modeled from plugin: bkg                                          SpectrumLike.py:482

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

         INFO      - observation: poisson                                                       SpectrumLike.py:492

         INFO      - background: poisson                                                        SpectrumLike.py:493


When we look at out count spectrum now, we will see the predicted background, rather than the measured one:

[11]:

modeled_background_plugin.view_count_spectrum()

[11]:


Now we simply fit the spectrum as we did in the profiled case. The background plugin’s parameters are stored in our new plugin as nuissance parameters:

[12]:

modeled_background_plugin.nuisance_parameters

[12]:

OrderedDict([('cons_full',
Parameter cons_full = 1.0 []
(min_value = 0.8, max_value = 1.2, delta = 0.05, free = False)),
('bkg_bkg_position_ra_full',
Parameter ra = 0.0 [deg]
(min_value = 0.0, max_value = 360.0, delta = 0.1, free = False)),
('bkg_bkg_position_dec_full',
Parameter dec = 0.0 [deg]
(min_value = -90.0, max_value = 90.0, delta = 0.1, free = False)),
('bkg_bkg_spectrum_main_composite_K_1_full',
Parameter K_1 = 0.32001053275945907 [1 / (cm2 keV s)]
(min_value = 0.0001, max_value = 10.0, delta = 0.1, free = True)),
('bkg_bkg_spectrum_main_composite_piv_1_full',
Parameter piv_1 = 100.0 [keV]
(min_value = None, max_value = None, delta = 0.1, free = False)),
('bkg_bkg_spectrum_main_composite_index_1_full',
Parameter index_1 = -1.3461600262984823 []
(min_value = -10.0, max_value = 10.0, delta = 0.20099999999999998, free = True)),
('bkg_bkg_spectrum_main_composite_F_2_full',
Parameter F_2 = 24.570007209793534 [1 / (cm2 s)]
(min_value = 0.0, max_value = 1000.0, delta = 0.1, free = True)),
('bkg_bkg_spectrum_main_composite_mu_2_full',
Parameter mu_2 = 511.0 [keV]
(min_value = None, max_value = None, delta = 0.1, free = False)),
('bkg_bkg_spectrum_main_composite_sigma_2_full',
Parameter sigma_2 = 19.777423534922427 [keV]
(min_value = 2.0, max_value = 30.0, delta = 0.1, free = True)),
('bkg_cons_bkg_full',
Parameter cons_bkg = 1.0 []
(min_value = 0.8, max_value = 1.2, delta = 0.05, free = False))])


and the fitting engine will use them in the fit. The parameters will still be connected to the background plugin and its model and thus we can free/fix them there as well as set priors on them.

[13]:

# instance the source model... the background plugin has it's model already specified
bpl = Broken_powerlaw(piv=300, xb=500)

bpl.K.bounds = (1e-5, 1e1)
bpl.xb.bounds = (1e1, 1e4)

ps_src = PointSource("source", 0, 0, bpl)

src_model = Model(ps_src)

jl_src = JointLikelihood(src_model, DataList(modeled_background_plugin))

_ = jl_src.fit()

08:30:54 INFO      set the minimizer to minuit                                             joint_likelihood.py:1042

Best fit values:


result unit
parameter
source.spectrum.main.Broken_powerlaw.K 1.77 +/- 0.13 1 / (cm2 keV s)
source.spectrum.main.Broken_powerlaw.xb (3.05 +/- 0.14) x 10^2 keV
source.spectrum.main.Broken_powerlaw.alpha (-8 +/- 8) x 10^-2
source.spectrum.main.Broken_powerlaw.beta -2.95 +/- 0.19
K_1 (3.40 +/- 0.31) x 10^-1 1 / (cm2 keV s)
index_1 -1.31 +/- 0.04
F_2 (2.5 +/- 0.5) x 10 1 / (cm2 s)
sigma_2 (2.14 +/- 0.31) x 10 keV

Correlation matrix:


 1 -0.62 0.75 0.07 0.09 -0.19 0 -0 -0.62 1 -0.47 -0.63 0.02 0.27 -0.07 -0.02 0.75 -0.47 1 0.05 0.33 -0.28 -0 -0 0.07 -0.63 0.05 1 -0.23 -0.39 -0.03 -0.02 0.09 0.02 0.33 -0.23 1 0.04 -0.02 -0 -0.19 0.27 -0.28 -0.39 0.04 1 -0 0 0 -0.07 -0 -0.03 -0.02 -0 1 0.14 -0 -0.02 -0 -0.02 -0 0 0.14 1

Values of -log(likelihood) at the minimum:


-log(likelihood)
full 544.959925
total 544.959925

Values of statistical measures:


statistical measures
AIC 1106.941126
BIC 1130.004931
[14]:

# over plot the joint background and source fits
fig = modeled_background_plugin.display_model(step=False)

_ = background_plugin.display_model(
data_color="#1A68F0", model_color="#FF9700", model_subplot=fig.axes, step=False
)

[ ]: