SED fitting with naima

In this notebook we will carry out a fit of an IC model to the HESS spectrum of RX J1713.7-3946 with the naima wrapper around emcee. This tutorial will follow loosely the tutorial found on the naima documentation.

The first step is to load the data, which we can find in the same directory as this notebook. The data format required by naima for the data files can be found in the documentation.



In [ ]:

    
import naima
import numpy as np
from astropy.io import ascii
import astropy.units as u
%matplotlib inline
import matplotlib.pyplot as plt

hess_spectrum = ascii.read('RXJ1713_HESS_2007.dat', format='ipac')
fig = naima.plot_data(hess_spectrum)

The we define the model to be fit. The model function must take a tuple of free parameters as first argument and a data table as second. It must return the model flux at the energies given by data['energy'] in first place, and any extra objects will be saved with the MCMC chain.

emcee does not accept astropy Quantities as parameters, so we have to give them units before setting the attributes of the particle distribution function.

Here we define an IC model with an Exponential Cutoff Power-Law with the amplitude, index, and cutoff energy as free parameters. Because the amplitude and cutoff energy may be considered to have a uniform prior in log-space, we sample their decimal logarithms (we could also use a log-uniform prior). We also place a uniform prior on the particle index with limits between -1 and 5.



In [ ]:

    
from naima.models import ExponentialCutoffPowerLaw, InverseCompton
from naima import uniform_prior

ECPL = ExponentialCutoffPowerLaw(1e36/u.eV, 5*u.TeV, 2.7, 50*u.TeV)
IC = InverseCompton(ECPL, seed_photon_fields=['CMB', ['FIR', 30*u.K, 0.4*u.eV/u.cm**3]])

# define labels and initial vector for the parameters
labels = ['log10(norm)', 'index', 'log10(cutoff)']
p0 = np.array((34, 2.7, np.log10(30)))

# define the model function
def model(pars, data):
    ECPL.amplitude = (10**pars[0]) / u.eV
    ECPL.alpha = pars[1]
    ECPL.e_cutoff = (10**pars[2]) * u.TeV

    return IC.flux(data['energy'], distance=2.0*u.kpc), IC.compute_We(Eemin=1*u.TeV)

from naima import uniform_prior

def lnprior(pars):
    lnprior = uniform_prior(pars[1], -1, 5)
    return lnprior

We take the data, model, prior, parameter vector, and labels and call the main fitting procedure: naima.run_sampler. This function is a wrapper around emcee, and the details of the MCMC run can be configured through its arguments:

nwalkers: number of emcee walkers.
nburn: number of steps to take for the burn-in period. These steps will be discarded in the final results.
nrun: number of steps to take and save to the sampler chain.
prefit: whether to do a Nelder-Mead fit before starting the MCMC run (reduces the burn-in steps required).
interactive: whether to launch an interactive model fitter before starting the run to set the initial vector. This will only work in matplotlib is using a GUI backend (qt4, qt5, gtkagg, tkagg, etc.). The final parameters when you close the window will be used as starting point for the run.
threads: How many different threads (CPU cores) to use when computing the likelihood.



In [ ]:

    
sampler, pos = naima.run_sampler(data_table=hess_spectrum, model=model, prior=lnprior, p0=p0, labels=labels,
                                nwalkers=32, nburn=50, nrun=100, prefit=True, threads=4)



In [ ]:

    
# inspect the chains stored in the sampler for the three free parameters
f = naima.plot_chain(sampler, 0)
f = naima.plot_chain(sampler, 1)
f = naima.plot_chain(sampler, 2)



In [ ]:

    
# make a corner plot of the parameters to show covariances
f = naima.plot_corner(sampler)



In [ ]:

    
# Show the fit
f = naima.plot_fit(sampler)
f.axes[0].set_ylim(bottom=1e-13)



In [ ]:

    
# Inspect the metadata blob saved
f = naima.plot_blob(sampler,1, label='$W_e (E_e>1$ TeV)')



In [ ]:

    
# There is also a convenience function that will plot all the above files to pngs or a single pdf
naima.save_diagnostic_plots('RXJ1713_naima_fit', sampler, blob_labels=['Spectrum','$W_e (E_e>1$ TeV)'])

Simultaneous fitting of two radiative components: Synchrotron and IC.

Use the Suzaku XIS spectrum of RX J1713 to do a simultaneous fit of the synchrotron and inverse Compton spectra and derive an estimate of the magnetic field strength under the assumption of a leptonic scenario.



In [ ]:

    
suzaku_spectrum = ascii.read('RXJ1713_Suzaku-XIS.dat')
f=naima.plot_data(suzaku_spectrum)

Note that in all naima functions (including run_sampler) you can provide a list of spectra, so you can consider both the HESS and Suzaku spectra:



In [ ]:

    
f=naima.plot_data([suzaku_spectrum, hess_spectrum], sed=True)

Below is the model, labels, parameters and prior defined above for the IC-only fit. Modify it as needed and feed it to naima.run_sampler to obtain an estimate of the magnetic field strength.



In [ ]:

    
#from naima.models import ExponentialCutoffPowerLaw, InverseCompton
#from naima import uniform_prior

#ECPL = ExponentialCutoffPowerLaw(1e36/u.eV, 10*u.TeV, 2.7, 50*u.TeV)
#IC = InverseCompton(ECPL, seed_photon_fields=['CMB', ['FIR', 30*u.K, 0.4*u.eV/u.cm**3]])

## define labels and initial vector for the parameters
#labels = ['log10(norm)', 'index', 'log10(cutoff)']
#p0 = np.array((34, 2.7, np.log10(30)))

## define the model function
#def model(pars, data):
#    ECPL.amplitude = (10**pars[0]) / u.eV
#    ECPL.alpha = pars[1]
#    ECPL.e_cutoff = (10**pars[2]) * u.TeV

#    return IC.flux(data['energy'], distance=2.0*u.kpc), IC.compute_We(Eemin=1*u.TeV)

#from naima import uniform_prior

#def lnprior(pars):
#    lnprior = uniform_prior(pars[1], -1, 5)
#    return lnprior