SpectralFitter

class gdt.core.spectra.fitting.SpectralFitter(pha_list, bkgd_list, rsp_list, statistic, channel_masks=None, method='SLSQP', rng=None)[source]

Bases: object

Base class for jointly-fitting spectra from multiple detectors. This is a (large) wrapper around scipy.optimize.minimize. This class should not be instantiated by the user, but instead a child class for a specific fit statistic should be instantiated.

Child classes should be designed with the following requirements:

  1. An __init__ method that takes in the inputs in the Parent __init__ method. Can have additional arguments/keywords.

  2. A _eval_stat method that accepts the index into the data sets and the source model rate for that data set. It must return a single number that represents the likelihood for a single detector/dataset.

Parameters:

  • pha_list (list of Pha) – The PHA objects containg the count spectrum for each detector

  • bkgd_list (list of BackgroundRates, list of BackgroundSpectrum, or list of Bak) – The background rates object, background spectrum, or Bak object for each detector. If given the background rates object, the times in the corresponding PHA object will be used for the limits of integration.

  • rsp_list (list of Rsp) – The response object for each detector

  • statistic (<function>) – The function to be used as a fit statistic

  • channel_masks (list of np.array(dtype=bool), optional) –

    A boolean mask for each detector that indicates which energy channels are to be fit.

    Note

    The channel_mask will typically be set by the valid_channels attribute in the PHA object. You can set channel_masks to override the PHA valid_channels, but this is usually unnecessary.

  • method (str, optional) –

    The fitting algorithm, which should be one of the options for scipy.optimize.minimize.

    Note

    All solvers, with the exception of ‘dogleg’ and ‘trust-exact’, are supported at this time.

  • rng (Generator, optional) – The RNG object

Attributes Summary

covariance

The covariance matrix of the fitted parameters.

detectors

The detectors used in the fit

dof

The number of degrees of freedom

energy_range

The full energy range of the fit

function_components

The names of the individual model components, if applicable

function_name

The name of the model function being fit

hessian

The Hessian matrix of the fitted parameters, for which the covariance matrix is -inverse(Hessian).

jacobian

The Jacobian vector of the likelihood as a function of the parameters.

message

The message from the fitter on convergence/non-convergence

num_components

Number of function components

num_sets

Number of datasets/detectors being fit

parameters

The parameter values at the maximum likelihood

statistic

The value of the fit statistic

success

True if the fit was successful, False otherwise.

symmetric_errors

The 1-sigma uncertainty on the parameters assuming higher-order terms of the Taylor expansion of the likelihood about the maximum likelihood solution is negligible.

Methods Summary

asymmetric_errors([cl])

Calculate the asymmetric uncertainty on each fitted parameter.

data_count_spectrum([upper_limits_sigma])

The observed source count spectrum, which is (total - background).

fit(function, **kwargs)

Fit the given the function to the data.

load(filename)

Loads a fitter object that has been saved to disk in a compressed numpy file.

model_count_spectrum()

The fitted model count spectrum.

model_variance()

The differential source model variance, accounting for the data variance and background model variance.

residuals([sigma])

The fit residuals for each dataset in units of differential counts or model sigma

sample_flux(erange[, num_samples, num_points])

Produce samples of the energy-integrated photon or energy flux.

sample_parameters(**kwargs)

Produce samples of the fitted parameters.

sample_spectrum(which[, num_samples, ...])

Produce samples of the model photon, energy, or \(\nu F_\nu\) spectrum.

save(filename)

Saves the fitter object to a compressed binary numpy file.

set_rng(rng)

Set/change the generator.

spectrum(which[, num_points, components])

The model photon, energy, or \(\nu F_\nu\) spectrum.

Attributes Documentation

covariance

The covariance matrix of the fitted parameters.

Note

The assumption of negligible higher-order terms in the Taylor expansion may not be valid and therefore the covariance matrix may not be valid.

Type:

(np.array)

detectors

The detectors used in the fit

Type:

(list of str)

dof

The number of degrees of freedom

Type:

(int)

energy_range

The full energy range of the fit

Type:

(float, float)

function_components

The names of the individual model components, if applicable

Type:

(list of str)

function_name

The name of the model function being fit

Type:

(str)

hessian

The Hessian matrix of the fitted parameters, for which the covariance matrix is -inverse(Hessian).

Note

The assumption of negligible higher-order terms in the Taylor expansion may not be valid and therefore the Hessian matrix may not be valid.

Type:

(np.array)

jacobian

The Jacobian vector of the likelihood as a function of the parameters.

Type:

(np.array)

message

The message from the fitter on convergence/non-convergence

Type:

(str)

num_components

Number of function components

Type:

(int)

num_sets

Number of datasets/detectors being fit

Type:

(int)

parameters

The parameter values at the maximum likelihood

Type:

(np.array)

statistic

The value of the fit statistic

Type:

(float)

success

True if the fit was successful, False otherwise.

Type:

(bool)

symmetric_errors

The 1-sigma uncertainty on the parameters assuming higher-order terms of the Taylor expansion of the likelihood about the maximum likelihood solution is negligible. If those terms are not negligible, errors could be very incorrect or event NaN. See attributes covariance and hessian.

Type:

(np.array)

Methods Documentation

asymmetric_errors(cl=0.683)[source]

Calculate the asymmetric uncertainty on each fitted parameter. This function uses scipy.optimize.brentq to find the root of:

crit = -2(LL_max-LL),

where crit is the critical value defined by the number of free parameters and the confidence level, and LL_max is the maximum log-likelihood. The brentq function uses a bracketing range to find the roots, and this function uses the best-fit values of the parameters for one side of the bracket and the min/max values defined in the fitted function for the other side of the bracket.

Parameters:

cl (float, optional) – The confidence level of the uncertainties. Default is 0.683.

Returns:

(np.array) –

A (nparams, 2) array where each row reprensents the

(negative, positive) uncertainty on each parameter.

data_count_spectrum(upper_limits_sigma=2.0)[source]

The observed source count spectrum, which is (total - background). Data bins with source counts less than the model variance are converted to upper limits.

Parameters:

upper_limits_sigma (float, optional) – The upper limits sigma. Default is 2.

Returns:

(tuple) – A 5-tuple containing, for each dataset:

  • list of np.array: Energy centroids

  • list of np.array: Energy channel half widths (in log-space)

  • list of np.array: differential count spectrum

  • list of np.array: 1-sigma count spectrum errors

  • list of np.array: Upper limits boolean mask

The upper limits mask is a boolean mask where True indicates the element is an upper limit.

fit(function, **kwargs)[source]

Fit the given the function to the data.

Parameters:

  • function (Function) – The function object to use

  • **kwargs – Keyword arguments passed to the fitter

classmethod load(filename)[source]

Loads a fitter object that has been saved to disk in a compressed numpy file. Thanks to the cool properties of the numpy format, this method will load the complete state of the fitter when it was saved.

Parameters:

filename (str) – The complete filename of the saved file

model_count_spectrum()[source]

The fitted model count spectrum.

Returns:

(list of EnergyBins)

model_variance()[source]

The differential source model variance, accounting for the data variance and background model variance.

Returns:

(list of np.array)

residuals(sigma=True)[source]

The fit residuals for each dataset in units of differential counts or model sigma

Parameters:

sigma (bool, optional) – If True, return in units of model sigma, otherwise return in units of differential counts. Default is True.

Returns:

(tuple) – A 4-tuple containing:

  • list of np.array: Energy centroids

  • list of np.array: Energy channel half-widths (in log-space)

  • list of np.array: Fit residuals

  • list of np.array: Residual uncertainties

sample_flux(erange, num_samples=1000, num_points=1000, **kwargs)[source]

Produce samples of the energy-integrated photon or energy flux. This uses the covariance matrix and assumes a multivariate normal distribuion for the parameters. Caveat Emptor.

Parameters:

  • erange (float, float) – The energy range over which to calculate the flux

  • num_samples (int) – The number of spectra to sample. Default is 1000.

  • num_points (int) – The number of grid points in energy. Default is 1000.

  • **kwargs – Options to be passed to integrate().

Returns:

(np.array)

sample_parameters(**kwargs)[source]

Produce samples of the fitted parameters. This uses the covariance matrix and assumes a multivariate normal distribuion. Caveat Emptor.

Parameters:

**kwargs – Options (like size) to be passed to numpy.random.multivariate_normal

Returns:

(np.array)

sample_spectrum(which, num_samples=1000, num_points=1000, components=False)[source]

Produce samples of the model photon, energy, or \(\nu F_\nu\) spectrum. This uses the covariance matrix and assumes a multivariate normal distribuion for the parameters. Caveat Emptor.

Parameters:

  • which (str) – Which spectrum to return. Either ‘photon’, ‘energy’, or ‘nufnu’

  • num_samples (int, optional) – The number of spectrum to sample. Default is 1000.

  • num_points (int, optional) – The number of grid points in energy. Default is 1000.

  • components (bool, optional) – If True, calculate the spectrum for each individual model components (if there are multiple components). Default is False.

Returns:

(tuple) – A 2-tuple containing:

  • np.array: The energy grid at which the spectrum is calculated

  • np.array: Array of shape (num_samples, num_points) if there is only one component or (num_samples, num_components, num_points) if there are multiple components.

save(filename)[source]

Saves the fitter object to a compressed binary numpy file. The full state of the fitter is serialized and saved, and it can be restored using the load() method.

Parameters:

filename (str) – The complete filename to save the fitter object to

set_rng(rng)[source]

Set/change the generator.

Parameters:

rng (numpy.random.Generator) – random number generator

spectrum(which, num_points=1000, components=False)[source]

The model photon, energy, or \(\nu F_\nu\) spectrum.

Parameters:

  • which (str) – Which spectrum to return. Either ‘photon’, ‘energy’, or ‘nufnu’

  • num_points (int, optional) – The number of grid points in energy. Default is 1000.

  • components (bool, optional) – If True, calculate the spectrum for each individual model components (if there are multiple components). Default is False.

Returns:

(tuple) – A 2-tuple containing:

  • np.array: The energy grid at which the spectrum is calculated

  • np.array or list of np.array: model spectrum values, or spectrum values for each component