Spectral Fitting

In the following you will see how to fit spectral data in OGIP format. The format is described at 1D counts spectra. An example dataset is available in the $GAMMAPY_DATA repo. For a description of the available fit statstics see Fit statistics.

Getting Started

The following example shows how to fit a power law simultaneously to two simulated crab runs using the Fit class.

from gammapy.spectrum import SpectrumDatasetOnOff
from gammapy.modeling import Fit
from gammapy.modeling.models import PowerLawSpectralModel
import matplotlib.pyplot as plt

path = "$GAMMAPY_DATA/joint-crab/spectra/hess/"
obs_1 = SpectrumDatasetOnOff.from_ogip_files(path + "pha_obs23523.fits")
obs_2 = SpectrumDatasetOnOff.from_ogip_files(path + "pha_obs23592.fits")

model = PowerLawSpectralModel(
    index=2,
    amplitude='1e-12  cm-2 s-1 TeV-1',
    reference='1 TeV',
)

obs_1.model = model
obs_2.model = model

fit = Fit([obs_1, obs_2])
result = fit.run()

model.parameters.covariance = result.parameters.covariance You can check the fit results by looking at the result and model object:

>>> print(result)

    OptimizeResult

    backend    : minuit
    method     : minuit
    success    : True
    nfev       : 115
    total stat : 65.36
    message    : Optimization terminated successfully.


>>> print(model)

    PowerLawSpectralModel

    Parameters:

           name     value     error        unit      min max frozen
        --------- --------- --------- -------------- --- --- ------
            index 2.781e+00 1.120e-01                nan nan  False
        amplitude 5.201e-11 4.965e-12 cm-2 s-1 TeV-1 nan nan  False
        reference 1.000e+00 0.000e+00            TeV nan nan   True

    Covariance:

           name     index   amplitude reference
        --------- --------- --------- ---------
            index 1.255e-02 3.578e-13 0.000e+00
        amplitude 3.578e-13 2.465e-23 0.000e+00
        reference 0.000e+00 0.000e+00 0.000e+00

Interactive Sherpa Fit

If you want to do something specific you can always fit the PHA data directly using Sherpa. The following example illustrates how to do this with the example dataset used above. It makes use of the Sherpa datastack module.

from pathlib import Path
import os
from sherpa.astro import datastack
from sherpa.models import PowLaw1D

pha1 = str(Path(os.environ["GAMMAPY_DATA"]) / "joint-crab/spectra/hess/pha_obs23592.fits")
pha2 = str(Path(os.environ["GAMMAPY_DATA"]) / "joint-crab/spectra/hess/pha_obs23523.fits")
phalist = ','.join([pha1, pha2])

ds = datastack.DataStack()
ds.load_pha(phalist)

model = PowLaw1D('powlaw1d.default')
model.ampl = 1
model.ref = 1e9
model.gamma = 2

ds.set_source(model*1e-20)

for i in range(1, len(ds.datasets) + 1):
    datastack.ignore_bad(i)
    datastack.ignore_bad(i, 1)

datastack.set_stat('wstat')
ds.fit()
datastack.covar()

This should give the following output

Datasets              = 1, 2
Method                = levmar
Statistic             = wstat
Initial fit statistic = 253.552
Final fit statistic   = 65.361 at function evaluation 25
Data points           = 82
Degrees of freedom    = 80
Probability [Q-value] = 0.88159
Reduced statistic     = 0.817012
Change in statistic   = 188.191
   powlaw1d.default.gamma   2.78053      +/- 0.121423
   powlaw1d.default.ampl   5.20034      +/- 0.510299
Datasets              = 1, 2
Confidence Method     = covariance
Iterative Fit Method  = None
Fitting Method        = levmar
Statistic             = wstat
covariance 1-sigma (68.2689%) bounds:
   Param            Best-Fit  Lower Bound  Upper Bound
   -----            --------  -----------  -----------
   powlaw1d.default.gamma      2.78053    -0.112025     0.112025
   powlaw1d.default.ampl      5.20034    -0.496564     0.496564