Tutorial

The first step in fitting a model to an observed spectrum is to read the spectrum into the appropriate format. See Data format for an explanation of the format and an example, and Units system for a brief explanation of the unit system used in Naima. We load the spectral data with astropy.io.ascii:

from astropy.io import ascii

data = ascii.read("RXJ1713_HESS_2007.dat")

Building the model and prior functions

The model function is the function that will be called to compare with the observed spectrum. It must take two parameters: an array of the free parameters of the model, and the data table.

Naima includes several models in the naima.models module that make it easier to fit common functional forms for spectra (PowerLaw, ExponentialCutoffPowerLaw, BrokenPowerLaw, and LogParabola), as well as several radiative models (InverseCompton, Synchrotron, Bremsstrahlung, and PionDecay; see Radiative Models for a detailed explanation of these). Once initialized with the relevant parameters, all model instances can be called with an energy array to obtain the flux of the model at the values of the energy array. If they are called with a data table as argument, the energy values from the energy column will be used.

Building the model function from one of the radiative models is easy. In the following example, the three model parameters in the pars array are the amplitude, the spectral index, and the logarith or the cutoff energy. We first add the necessary units for the radiative model and then compute and return the model flux for the energies contained in the data table:

from naima.models import ExponentialCutoffPowerLaw, InverseCompton
import astropy.units as u


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

    ECPL = ExponentialCutoffPowerLaw(amplitude, 10 * u.TeV, alpha, e_cutoff)
    IC = InverseCompton(
        ECPL,
        seed_photon_fields=[
            "CMB",
            ["FIR", 26.5 * u.K, 0.415 * u.eV / u.cm ** 3],
        ],
    )

    return IC.flux(data, distance=1.0 * u.kpc)

In addition, we must build a function to return the prior function, i.e., a function that encodes any previous knowledge you have about the parameters, such as previous measurements or physically acceptable ranges. Two simple priors functions are included with Naima: normal_prior, and uniform_prior, and loguniform_prior. uniform_prior can be used to set parameter limits. Following the example above, we might want to limit the amplitude to be positive, and the spectral index to be between -1 and 5:

from naima import uniform_prior


def lnprior(pars):
    lnprior = uniform_prior(pars[0], 0.0, np.inf) + uniform_prior(pars[1], -1, 5)
    return lnprior

Selecting a starting point for the sampling

Before starting the MCMC run, we must provide the procedure with initial estimates of the parameters and their names:

p0 = np.array((1e33, 3.0, np.log10(30)))
labels = ["norm", "index", "log10(cutoff)"]

This example is relatively simple, but for more complicated models (in particular those with more than one emission channel, such as Synchrotron and Inverse Compton) it may be difficult to provide an adequate initial parameter vector without comparing it visually with the spectra. For these models, an initial parameter vector far from the maximum likelihood vector may mean that during the minimization or sampling process the algorithm gets stuck in a local maximum.

In order to make an adequate estimation easier, Naima provides a tool to interactively see the output of the model and compare it with the observed data while changing the parameter values. This tool can be accessed in two ways: the first is setting interactive=True in the options of get_sampler or run_sampler (see below in sampling for details on these functions). This will launch an interactive matplotlib window in which the parameters can be changed and the resulting model compared to the observed spectra. With each change, the log probability of the model given the observed spectra is computed. In addition, a Nelder-Mead fit can be launched from a button in the window to find the maximum likelihood parameter vector. Once you are happy that the current model is a good approximation to the observed spectrum, closing the window (whether through the window manager or the Close window button) will use the current parameter vector as a starting point for the sampling procedure.

The alternative way of accessing the interactive fitter is to access it directly through the class naima.InteractiveModelFitter from an interactive Python interpreter. The parameter vector shown in the interactive window can be accessed through the imf.pars attribute, and then copied to a new p0 variable to be used in the sampling:

>> imf = InteractiveModelFitter(model, p0, data=data, labels=labels)
>> # interactive fitting done
>> p0 = imf.pars

Note that the data argument can be omitted and an energy range specified through the e_range argument to inspect the behaviour of the model independently of the data:

>> imf = InteractiveModelFitter(model, p0, e_range=[1 * u.GeV, 100 * u.TeV])

Sampling the posterior distribution function

All the objects above can then be provided to run_sampler, the main fitting function in Naima:

sampler, pos = naima.run_sampler(
    data_table=data,
    p0=p0,
    label=labels,
    model=model_function,
    prior=lnprior,
    nwalkers=128,
    nburn=50,
    nrun=10,
    threads=4,
)

The nwalkers parameter specifies how many walkers will be used in the sampling procedure, nburn specifies how many steps to be run as burn-in, and nrun specifies how many steps to run after the burn-in and save these samples in the sampler object. For details on these parameters, see the documentation of the emcee package.

Inspecting and analysing results of the run

The results stored in the sampler object can be analysed through the plotting procedures of Naima: plot_chain, plot_fit, and plot_data. In addition, two convenience functions can be used to generate a collection of plots that illustrate the results and the stability of the fitting procedure. These are save_diagnostic_plots and save_results_table:

naima.save_diagnostic_plots(
    "RXJ1713_IC",
    sampler,
    blob_labels=[
        "Spectrum",
        "Electron energy distribution",
        "$W_e (E_e>1$ TeV)",
    ],
)

naima.save_results_table("RXJ1713_naima_fit", sampler)

The saved table will include information in the metadata about the run such as the number of walkers n_walkers and steps n_run sampled, the initial parameter vectori p0, the parameter vector with the maximum likelihood ML_pars and the maximum value of the negative log-likelihood MaxLogLikelihood. The table itself shows the median and upper and lower uncertainties (50th, 84th, and 16th percentiles of the posterior distribution) for the parameters sampled:

n_walkers         : 128
n_run             : 100
p0                : [1.3884430349732936e+32, 2.5717536865466712, 1.6821145924203185]
ML_pars           : [1.3697204402514948e+32, 2.5839150825284958, 1.7002798209378214]
MaxLogLikelihood  : -17.98655803890747

------------- ------- ------- --------
        label   median   unc_lo  unc_hi
------------- -------- -------- -------
         norm 1.38e+32 8.57e+30 1.1e+31
        index     2.57    0.115   0.103
log10(cutoff)     1.69    0.109   0.115
       cutoff     48.9     10.9    14.8

The table is saved by default in ECSV format which can be easily accesed with the astropy.io.ascii module.

Plotting functions: chains

The function plot_chain will show information about the MCMC chain for a given parameter. It shows the parameter value with respect to the step number of the chain, which can be used to assess the stability of the chain, a plot of the posterior distribution, and several statistics of the posterior distribution. One of these is the median and 16th and 84th percentiles of the distribution, which can be reported as the inferred marginalised value of the parameter and associated \(1\sigma\) uncertainty. For the electron index (parameter 1), plot_chain shows:

_images/RXJ1713_IC_chain_index.png

For parameters that have been sampled in logarithmic space and their parameter label includes log10 or log, plot_chain will also compute the value and percentiles in linear space:

_images/RXJ1713_IC_chain_cutoff.png

The relationship between the samples of the different parameters can be seen though a corner plot with plot_corner which is a wrapper around corner.corner. The maximum likelihood parameter vector can be indicated with cross:

_images/RXJ1713_IC_corner.png

Plotting functions: fit

The plot function plot_fit allows for several ways to represent the results of the MCMC fitting. By default, it will show the Maximum Likelihood model with a black line, and 100 samples from the posterior distribution in gray:

_images/RXJ1713_IC_model_samples.png

The 100 samples are taken from the blobs stored in the sampler, so they only contain the model values at the observed flux points. If you want to show the samples and ML model for a wider energy range (or between energy bands like X-ray and gamma-ray) you can use the e_range parameter. Note that the model will be recomputed n_samples times (by default 100) when plot_fit is called, so this may significantly slow the plot speed if the model function calls are expensive. Setting e_range=[100*u.GeV, 100*u.TeV], we obtain the following plot:

_images/RXJ1713_IC_model_samples_erange.png

The spread of the parameters in the posterior distribution can also be visualized as confidence bands. Using the confs parameter of plot_fit, a confidence band will be computed for each of the confidence levels (in sigma) given in confs. Setting confs=[3,1], the confidence bands at \(1\sigma\) and \(3\sigma\) are plotted. Note that no samples are shown if the confs parameter is set:

_images/RXJ1713_IC_model_confs.png

As for the plot showing the samples, the energy range for the confidence bands can be set through the e_range parameter. The number of samples needed will be computed so that the highest confidence level given can be constrained. This results in 740 samples for a \(3\sigma\) confidence level:

_images/RXJ1713_IC_model_confs_erange.png

Saving and retrieving the results of the sampling run

The parameter chain and metadata blobs resulting from the sampling procedure can be saved with naima.save_run. This function will save the results of the run to an HDF5 file that can be archived and analysed after the fact. The sampler properties saved to the HDF5 file are:

  • parameter vector chain: sampler.get_chain()

  • log-probability for all the parameter vectors in the chain: sampler.get_log_prob()

  • metadata blobs: sampler.get_blobs()

  • parameter labels: sampler.labels

  • data table: sampler.data

Note that only metadata blobs that can be converted into a numpy array will be stored. Blobs consisting of other classes, or with different lengths within the same blob, will raise a warning and may be discarded on the naima.save_run call.

The saved sampler can be retrieved with the naima.read_run function, which will return an EnsembleSampler-like object. However, the model function cannot be saved in the HDF5 file, so a model function has to be provided at read time. This function will be used to compute the model output from the parameter vectors in the chain. Without a model function, the sampler read from the HDF5 file can be passed onto plot_chain, plot_corner, plot_fit, and plot_data for analysis. If a model function is provided to naima.read_run, the resulting sampler can also be used with plot_fit when setting the e_range parameter, which requires the model function to be accesible in the modelfn attribute of the sampler object provided.

Saving additional information — Metadata blobs

If we wish to save additional information at each of the model computations, extra information can be returned from the model call. This extra information (known as metadata blobs; see details in the emcee documentation) is stored in the sampler object returned from the fitting procedure and can be accessed later. There are three formats for the data stored as a metadata blob that will be understood by the plotting routines of Naima:

  • A Quantity scalar. A histogram and distribution properties (median, 16th and 84th percentiles, etc.) will be plotted.

  • A Quantity array with the same length as the observed spectrum energy array. If it has a physical type of flux or luminosity, it will be interpreted as a photon spectrum and plotted against the observed spectrum energy array.

  • A pair (tuple or list) of Quantity arrays of equal length. They will be plotted against each other.

When fitting a radiative output to a spectrum, information on the particle distribution (e.g., the actual particle distribution, or the total energy in relativistic particles) can be saved as a metadata blob. Below is an example that does precisely this with an Inverse Compton emission model:

from naima.models import ExponentialCutoffPowerLaw, InverseCompton
import astropy.units as u
import numpy as np


def model_function(pars, data):
    amplitude = pars[0] * (1 / u.eV)
    alpha = pars[1]
    e_cutoff = (10 ** pars[2]) * u.TeV
    e_0 = 10 * u.TeV

    ECPL = ExponentialCutoffPowerLaw(amplitude, e_0, alpha, e_cutoff)
    IC = InverseCompton(
        ECPL,
        seed_photon_fields=[
            "CMB",
            ["FIR", 26.5 * u.K, 0.415 * u.eV / u.cm ** 3],
        ],
    )

    # The total enegy in electrons of model IC can be accessed through the
    # attribute We or obtained for a given range with compute_We
    We = IC.compute_We(Eemin=1 * u.TeV)

    # We can also save the particle distribution between 100 MeV and 100 TeV
    electron_e = np.logspace(11, 15, 100) * u.eV
    electron_dist = ECPL(electron_e)

    # The first object returned must be the model photon spectrum, and
    # subsequent objects will be stored as metadata blobs
    return IC(data), (electron_e, electron_dist), We

The additional quantities we have stored can the be accesed via the sampler.get_blobs method. The function plot_blob allows to plot them and extract distribution properties. For the blobs that are a tuple or have the same length as data['energy'], they will be plotted as spectra:

_images/RXJ1713_IC_pdist.png

and for the ones that are a scalar value, such as the total energy in electrons that we returned as the third object, a histogram and distribution properties will be plotted:

_images/RXJ1713_IC_We.png