As an example of a chi^2 fit, we use the flux points from the Crab nebula as measured by H.E.S.S. in 2006: http://adsabs.harvard.edu/abs/2006A%26A...457..899A
In [1]:
%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
plt.style.use('ggplot')
We start by loading the flux points from a text file.
It is of course possible to load this data using just Python or Numpy, but we'll use the pandas.read_table function because it's very flexible, i.e. by setting a few arguments you'll be able to load most ascii tables.
In [2]:
data = pd.read_table('spectrum_crab_hess_2006.txt',
comment='#', sep='\s*', engine='python')
data
Out[2]:
In the paper they fit a power-law with an exponential cutoff and find the following parameters (see row "all" in table 6):
The flux1 is the differential flux at 1 TeV.
Let's code up that model ...
TODO: extend this tutorial to also consider a power-law model and compare the two models via chi2 / ndf.
In [3]:
def flux_ecpl(energy, flux1, gamma, energy_cut):
return flux1 * energy ** (-gamma) * np.exp(-energy / energy_cut)
In [4]:
energy = np.logspace(-0.5, 1.6, 100)
flux = flux_ecpl(energy, flux1=3.76e-11, gamma=2.39, energy_cut=14.3)
plt.plot(energy, flux)
plt.errorbar(x=data['energy'],
y=data['flux'],
yerr=data['flux_err'],
fmt='.'
)
plt.loglog();
In this case we'll use a chi^2 likelihood function to fit the model to the data.
Note that the likelihood function combines the data and model, and just depends on the model parameters that
shall be estimated (whereas the model function flux_ecpl has an extra parameter energy).
Also note that we're accessing data and model flux_ecpl from the global scope instead of passing them in explicitly as parameters. Modeling and fitting frameworks like e.g. Sherpa have more elaborate ways to combine data and models and likelihood functions, but for simple, small code examples like we do here, using the global scope to tie things together works just fine.
In [5]:
def chi2(flux1, gamma, energy_cut):
energy = data['energy']
flux_model = flux_ecpl(energy, flux1, gamma, energy_cut)
chi = (data['flux'] - flux_model) / data['flux_err']
return np.sum(chi ** 2)
In [6]:
# TODO: visualise the likelihood as a 1D profile or
# 2D contour to check that the implementation is OK
# before fitting. E.g. reproduce Fig 19 from the paper?
# Maybe talk about how chi2 differences relate to
# confidence levels here?
In [7]:
from iminuit import Minuit
# TODO
In [8]:
import emcee
# TODO
Often the fit parameters are not the quantities quoted when reporting the analysis results. Instead some other quantity that is a function of the fit parameters is of interest.
In the Crab paper, the integral flux above 1 TeV is the main flux quantity that's reported int the abstract.
Let's see how we can compute it and propagate the error on that quantity.
TODO:
In [ ]: