Spectra Module

The classes and methods in gbm.spectra are for spectral fitting.

Spectral Fitting

These classes and functions define the likelihood functions and maximum-likelihood fitters associated with those functions.

SpectralFitter

class gbm.spectra.fitting.SpectralFitter(pha_list, bkgd_list, rsp_list, statistic, channel_masks=None, method='SLSQP')

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 _fold_model method that accepts the model function and the list of parameters. 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 or list of BackgroundSpectrum) – The background rates object for each detector or background spectrum 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 that ‘TNC’ and ‘SLSQP’ are the two algorithms that provide support for fitting with bounds and constraints. Other algorithms will ignore bounds and constraints.

Attributes

• covariance (np.array) – 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.

• detectors (list of str) – The detectors used in the fit

• dof (int) – The number of degrees of freedom

• energy_range (float, float) – The full energy range of the fit

• function_components (list of str) – The names of the individual model components, if appropriate

• function_name (str) – The name of the model function being fit

• hessian (np.array) – 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.

• jacobian (np.array) – The Jacobian vector of the likelihood as a function of the parameters.

• message (str) – The message from the fitter on convergence/non-convergence

• num_components (int) – Number of function components

• num_sets (int) – Number of datasets/detectors being fit

• parameters (np.array) – The parameter values at the maximum likelihood

• statistic (float) – The fit statistic (likelihood)

• success (bool) – True if the fit was successful, False otherwise.

• symmetric_errors (np.array) – 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.

asymmetric_errors(cl=0.683)

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 chi-square 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 (1 sigma).

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)

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:

• list of np.array: Energy centroids

• list of np.array: Energy channel widths

• list of np.array: differential count spectrum

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

• list of np.array: Upper limits boolean mask

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

fit(function, **kwargs)

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)

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()

The fitted model count spectrum

Returns

list of EnergyBins – An EnergyBins object for each dataset

model_variance()

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

Returns

list of np.array – The model variance for each dataset.

residuals(sigma=True)

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 widths

• list of np.array: Fit residuals

• list of np.array: Residual uncertainties

sample_flux(erange, num_samples=1000, numpoints=1000, **kwargs)

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) – Which spectrum to return. Either ‘photon’, ‘energy’, or ‘nufnu’

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

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

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

Returns

np.array – The sampled flux

sample_parameters(**kwargs)

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 – The random parameter samples

sample_spectrum(which, num_samples=1000, numpoints=1000, components=False)

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.

• numpoints (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, numpoints) if there is only one component or (num_samples, num_components, numpoints) if there are multiple components.

save(filename)

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

spectrum(which, numpoints=1000, components=False)

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

Parameters

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

• numpoints (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

---

SpectralFitterChisq

class gbm.spectra.fitting.SpectralFitterChisq(pha_list, bkgd_list, rsp_list, **kwargs)

Bases: gbm.spectra.fitting.SpectralFitter

Jointly-fit datasets using Chi-Squared Likelihood.

Parameters

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

• bkgd_list (list of BackgroundRates) – The background rates object for each detector

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

• 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 that ‘TNC’ and ‘SLSQP’ are the two algorithms that provide support for fitting with bounds and constraints. Other algorithms will ignore bounds and constraints.

Attributes

• covariance (np.array) – 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.

• detectors (list of str) – The detectors used in the fit

• dof (int) – The number of degrees of freedom

• energy_range (float, float) – The full energy range of the fit

• function_components (list of str) – The names of the individual model components, if appropriate

• function_name (str) – The name of the model function being fit

• hessian (np.array) – 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.

• jacobian (np.array) – The Jacobian vector of the likelihood as a function of the parameters.

• message (str) – The message from the fitter on convergence/non-convergence

• num_components (int) – Number of function components

• num_sets (int) – Number of datasets/detectors being fit

• parameters (np.array) – The parameter values at the maximum likelihood

• statistic (float) – The fit statistic (likelihood)

• success (bool) – True if the fit was successful, False otherwise.

• symmetric_errors (np.array) – 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.

---

chisq

gbm.spectra.fitting.chisq(obs_counts, back_rates, back_var, mod_rates, exposure)

Chi-Square Likelihood

Parameters

• obs_counts (np.array) – The total observed counts

• back_rates (np.array) – The background model count rates

• back_var (np.array) – The background model count rate variance

• mod_rates (np.array) – The model source rates

• exposure (float) – The source exposure

Returns

float – The chi-squared statistic

SpectralFitterCstat

class gbm.spectra.fitting.SpectralFitterCstat(pha_list, bkgd_list, rsp_list, back_exp=None, **kwargs)

Bases: gbm.spectra.fitting.SpectralFitter

Jointly-fit datasets using C-Stat (Poisson source with Poisson background)

Parameters

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

• bkgd_list (list of BackgroundRates) – The background rates object for each detector

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

• 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.

• back_exp (list of float, optional) – Set if different than the source exposure

• method (str, optional) – The fitting algorithm, which should be one of the options for scipy.optimize.minimize. Note that ‘TNC’ and ‘SLSQP’ are the two algorithms that provide support for fitting with bounds and constraints. Other algorithms will ignore bounds and constraints.

Attributes

• covariance (np.array) – 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.

• detectors (list of str) – The detectors used in the fit

• dof (int) – The number of degrees of freedom

• energy_range (float, float) – The full energy range of the fit

• function_components (list of str) – The names of the individual model components, if appropriate

• function_name (str) – The name of the model function being fit

• hessian (np.array) – 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.

• jacobian (np.array) – The Jacobian vector of the likelihood as a function of the parameters.

• message (str) – The message from the fitter on convergence/non-convergence

• num_components (int) – Number of function components

• num_sets (int) – Number of datasets/detectors being fit

• parameters (np.array) – The parameter values at the maximum likelihood

• statistic (float) – The fit statistic (likelihood)

• success (bool) – True if the fit was successful, False otherwise.

• symmetric_errors (np.array) – 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.

---

cstat

gbm.spectra.fitting.cstat(obs_counts, mod_rates, exposure, back_counts, back_exp)

C-Statistic Likelihood for Poisson source + Poisson background. The “W” statistic from the XSPEC Statistics Appendix.

Parameters

• obs_counts (np.array) – The total observed counts

• mod_rates (np.array) – The model source rates

• exposure (float) – The source exposure

• back_counts (np.array) – The background model counts

• back_exp (float) – The background exposure

Returns

float – The C-statistic

SpectralFitterPgstat

class gbm.spectra.fitting.SpectralFitterPgstat(pha_list, bkgd_list, rsp_list, back_exp=None, **kwargs)

Bases: gbm.spectra.fitting.SpectralFitter

Jointly-fit datasets using Profile Gaussian likelihood (PGSTAT).

Parameters

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

• bkgd_list (list of BackgroundRates) – The background rates object for each detector

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

• 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.

• back_exp (list of float, optional) – Set if different than the source exposure

• method (str, optional) – The fitting algorithm, which should be one of the options for scipy.optimize.minimize. Note that ‘TNC’ and ‘SLSQP’ are the two algorithms that provide support for fitting with bounds and constraints. Other algorithms will ignore bounds and constraints.

Attributes

• covariance (np.array) – 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.

• detectors (list of str) – The detectors used in the fit

• dof (int) – The number of degrees of freedom

• energy_range (float, float) – The full energy range of the fit

• function_components (list of str) – The names of the individual model components, if appropriate

• function_name (str) – The name of the model function being fit

• hessian (np.array) – 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.

• jacobian (np.array) – The Jacobian vector of the likelihood as a function of the parameters.

• message (str) – The message from the fitter on convergence/non-convergence

• num_components (int) – Number of function components

• num_sets (int) – Number of datasets/detectors being fit

• parameters (np.array) – The parameter values at the maximum likelihood

• statistic (float) – The fit statistic (likelihood)

• success (bool) – True if the fit was successful, False otherwise.

• symmetric_errors (np.array) – 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.

---

pgstat

gbm.spectra.fitting.pgstat(obs_counts, mod_rates, exposure, back_counts, back_var, back_exp)

Profile Gaussian Likelihood. From the XSPEC Statistics Appendix.

Parameters

• obs_counts (np.array) – The total observed counts

• mod_rates (np.array) – The model source rates

• exposure (float) – The source exposure

• back_counts (np.array) – The background model counts

• back_var (np.array) – The background model count variance

• back_exp (float) – The background exposure

Returns

float – The PG-stat

---

SpectralFitterPstat

class gbm.spectra.fitting.SpectralFitterPstat(pha_list, bkgd_list, rsp_list, **kwargs)

Bases: gbm.spectra.fitting.SpectralFitter

Jointly-fit datasets using P-Stat (Poisson source with known background). This statistic assumes non-zero counts, therefore any channels with zero counts will be masked out and not used in fitting.

Parameters

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

• bkgd_list (list of BackgroundRates) – The background rates object for each detector

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

• 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.

• back_exp (list of float, optional) – Set if different than the source exposure

• method (str, optional) – The fitting algorithm, which should be one of the options for scipy.optimize.minimize. Note that ‘TNC’ and ‘SLSQP’ are the two algorithms that provide support for fitting with bounds and constraints. Other algorithms will ignore bounds and constraints.

Attributes

• covariance (np.array) – 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.

• detectors (list of str) – The detectors used in the fit

• dof (int) – The number of degrees of freedom

• energy_range (float, float) – The full energy range of the fit

• function_components (list of str) – The names of the individual model components, if appropriate

• function_name (str) – The name of the model function being fit

• hessian (np.array) – 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.

• jacobian (np.array) – The Jacobian vector of the likelihood as a function of the parameters.

• message (str) – The message from the fitter on convergence/non-convergence

• num_components (int) – Number of function components

• num_sets (int) – Number of datasets/detectors being fit

• parameters (np.array) – The parameter values at the maximum likelihood

• statistic (float) – The fit statistic (likelihood)

• success (bool) – True if the fit was successful, False otherwise.

• symmetric_errors (np.array) – 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.

---

pstat

gbm.spectra.fitting.pstat(obs_counts, mod_rates, exposure, back_rates)

Likelihood for Poisson source + known background. The “pstat” statistic from the XSPEC Statistics Appendix.

Note

Elements with zero counts are masked out and not figured in the statistic.

Parameters

• obs_counts (np.array) – The total observed counts

• mod_rates (np.array) – The model source rates

• exposure (float) – The source exposure

• back_rates (np.array) – The background model count rates

Returns

float – The C-statistic

---

Spectral Function Base Classes

The base classes represent photon model functions with metadata that allow these functions to be fit to data. Specifically, Function is the base class for all functions and SuperFunction inherits from Function and is the base class for either a sum or product of functions.

Function

class gbm.spectra.functions.Function

Bases: object

Abstract class for a 1-D function (e.g. photon model).

Sub-classes should be designed with the following requirements:

1. Attribute nparams containing the total number of parameters of the function (even parameters that are fixed to a particular value during fitting)

2. Method eval() that accepts two arguments: an array of parameter values of length nparams; and the array of abscissa values.

Sub-classes can optionally add the following for complete and enhanced functionality:

1. Attribute param_list that is a list of parameter names in the order they are accepted in the eval() method. Each entry in param_list is a tuple of the form (<param_name>, <units>, <description>). If no units or description are needed, those fields should be a null string ‘’. If omitted, these will be filled to default values (‘param0’, ‘’, ‘’), etc.

2. Attribute default_values that is a list of default parameter values. This is primarily used for parameter estimation. Default is 1.0 for each parameter

3. Attribute delta_abs that is a list of max absolute fitting step for each parameter. This is used for parameter estimation. Default is 0.1 for each parameter.

4. Attribute delta_rel that is a list of max relative fitting step for each parameter. This is used for paramter estimation. Default is 0.01 for each parameter.

5. Attribute min_values that is a list of minimum values each parameter can take. This is used primarily for parameter estimation. Default is -np.inf for each parameter.

6. Attribute max_values that is a list of maximum values each parameter can take. This is used primarily for parameter estimation. Default is np.inf for each parameter.

7. Attribute free that is a boolean list indicating if each parameter is allowed to be free for fitting. False indicates that the parameter will be fixed at its default value. This is used for parameter estimation. Default is True for each parameter.

Attributes

• default_values (list) – The default values for each parameter. If not defined, each defaults to 1.0

• delta_abs (list) – The maximum absolute fitting step size for each parameter. If not defined, each default to 0.1

• delta_rel (list) – The maximum relative fitting step size for each parameter. If not defined, each default to 0.01

• free (list) – For each parameter, mark True if left free for fitting, otherwise mark False to fix at the default value. If not defined, each default to True.

• max_values (list) – The maximum possible value for each parameter. If not defined, each default to np.inf

• min_values (list) – The minimum possible value for each parameter. If not defined, each default to -np.inf

• name (str) – The name of the function

• nparams (int) – Number of parameters in the function.

• num_components (int) – Number of function components. For a normal function this is one, for a SuperFunction, this can be > 1.

• param_list (list of 3-tuples) – A list of tuples corresponding to the metadata for the parameters: (<param_name>, <units>, <description>)

classmethod add(func1, func2)

Add two Functions. Can also use the + operator.

Parameters

• func1 (Function) – A function

• func2 (Function) – A function to add to func1.

Returns

SuperFunction – A function containing the sum of the two functions

eval(params, x)

Evaluate the function

Parameters

• params (iterable) – The parameter vector of the free parameters

• x (np.array) – The values at which to evaluate the function

Returns

np.array – The function evaluated at the x points

fit_eval(params, x, **kwargs)

Evaluate the function for a fit. This differs from the eval() because the params vector should only contain the free fit parameters. The parameters that are set to fixed will automatically be passed to eval() at their fixed value. This enables a fit to be performed withfixed parameters and the associated jacobian and covariance matrices will have the appropriate shape.

Parameters

• params (iterable) – The parameter vector of the free parameters

• x (np.array) – The values at which to evaluate the function

Returns

np.array – The function evaluated at the x points

integrate(params, erange, num_points=1000, log=True, energy=False)

Integrate the function over energy to produce a flux.

Parameters

• params (iterable) – The parameter vector, which can be of length nparams or equal to the length of the non-fixed (free) parameters.

• erange (float, float) – The energy range over which to integrate

• num_points (int, optional) – The number of points used for integration. Default is 1000.

• log (bool, optional) – If True, the integration grid is made in log-space, False it is linear. Default is True.

• energy (bool, optional) – If True, perform integration of E*F(E) to produce an energy flux, otherwise integration is over F(E) to produce a photon flux. Default is False.

Returns

float – The photon or energy flux

classmethod multiply(func1, func2)

Multiply two Functions. Can also use the * operator.

Parameters

• func1 (Function) – A function

• func2 (Function) – A function to multiply to func1.

Returns

SuperFunction – A function containing the product of the two functions

parameter_bounds(apply_state=True)

The valid parameter bounds

Parameters

apply_state (bool, optional) – If True, will only return the bounds for free parameters. Default is True.

Returns

[(float, float), …] – The list of bounds [(min, max), …]

---

SuperFunction

class gbm.spectra.functions.SuperFunction

Bases: gbm.spectra.functions.Function

Abstract class for a sum or product of functions.

Attributes

• default_values (list) – The default values for each parameter. If not defined, each defaults to 1.0

• delta_abs (list) – The maximum absolute fitting step size for each parameter. If not defined, each default to 0.1

• delta_rel (list) – The maximum relative fitting step size for each parameter. If not defined, each default to 0.01

• free (list) – For each parameter, mark True if left free for fitting, otherwise mark False to fix at the default value. If not defined, each default to True.

• max_values (list) – The maximum possible value for each parameter. If not defined, each default to np.inf

• min_values (list) – The minimum possible value for each parameter. If not defined, each default to -np.inf

• name (str) – The name of the function

• nparams (int) – Number of parameters in the function.

• num_components (int) – Number of function components. For a normal function this is one, for a SuperFunction, this can be > 1.

• param_list (list of 3-tuples) – A list of tuples corresponding to the metadata for the parameters: (<param_name>, <units>, <description>)

classmethod add(func1, func2)

Add two Functions. Can also use the + operator.

Parameters

• func1 (Function) – A function

• func2 (Function) – A function to add to func1.

Returns

SuperFunction – A function containing the sum of the two functions

eval(params, x, components=False)

Evaluate the function

Parameters

• params (iterable) – The parameter vector of the free parameters

• x (np.array) – The values at which to evaluate the function

• components (bool, optional) – If True, evaluate for each component. Default is False.

Returns

np.array or list of np.array – The function evaluated at x points or the function evaluated at x points for each component.

fit_eval(params, x, **kwargs)

Evaluate the function for a fit. This differs from the eval() because the params vector should only contain the free fit parameters. The parameters that are set to fixed will automatically be passed to eval() at their fixed value. This enables a fit to be performed withfixed parameters and the associated jacobian and covariance matrices will have the appropriate shape.

Parameters

• params (iterable) – The parameter vector of the free parameters

• x (np.array) – The values at which to evaluate the function

Returns

np.array – The function evaluated at the x points

integrate(params, erange, num_points=1000, log=True, energy=False)

Integrate the function over energy to produce a flux.

Parameters

• params (iterable) – The parameter vector, which can be of length nparams or equal to the length of the non-fixed (free) parameters.

• erange (float, float) – The energy range over which to integrate

• num_points (int, optional) – The number of points used for integration. Default is 1000.

• log (bool, optional) – If True, the integration grid is made in log-space, False it is linear. Default is True.

• energy (bool, optional) – If True, perform integration of E*F(E) to produce an energy flux, otherwise integration is over F(E) to produce a photon flux. Default is False.

Returns

float – The photon or energy flux

classmethod multiply(func1, func2)

Multiply two Functions. Can also use the * operator.

Parameters

• func1 (Function) – A function

• func2 (Function) – A function to multiply to func1.

Returns

SuperFunction – A function containing the product of the two functions

parameter_bounds(apply_state=True)

The valid parameter bounds

Parameters

apply_state (bool, optional) – If True, will only return the bounds for free parameters. Default is True.

Returns

[(float, float), …] – The list of bounds [(min, max), …]

---

Spectral Functions

The functions listed below inherit from Function and have the optional metadata described in that class documentation.

PowerLaw

class gbm.spectra.functions.PowerLaw

Bases: gbm.spectra.functions.Function

A power law function:

\(F(E) = A \bigl(\frac{E}{E_{piv}}\bigr)^{\lambda}\),

where

• \(A\) is the amplitude in ph/s/cm^2/keV

• \(\lambda\) is the power-law index

• \(E_{piv}\) is the pivot energy in keV. Nominally fixed to 100 keV.

---

Comptonized

class gbm.spectra.functions.Comptonized

Bases: gbm.spectra.functions.Function

An exponentially cut-off power law, parametrized by the peak of the \(\nu F_{\nu}\) spectrum:

\(F(E) = A e^{-(2+\lambda)E/E_{peak}} \bigl(\frac{E}{E_{piv}}\bigr)^\lambda\)

where

• \(A\) is the amplitude in ph/s/cm^2/keV

• \(E_{peak}\) is the \(\nu F_{\nu}\) peak in keV

• \(\lambda\) is the power-law index

• \(E_{piv}\) is the pivot energy in keV. Nominally fixed to 100 keV.

---

Band

class gbm.spectra.functions.Band

Bases: gbm.spectra.functions.Function

The Band GRB function (a type of smoothly broken power law), parametrized by the peak of the \(\nu F_{\nu}\) spectrum:

\(F(E) = A \begin{cases} \bigl(\frac{E}{E_{piv}}\bigr)^\alpha e^{-(2+\alpha)E/E_{peak}}, \text{ if } {E < \frac{(\alpha-\beta)E_{peak}}{2+\alpha}}\\ \bigl(\frac{(\alpha-\beta)E_{peak}}{(2+\alpha)E_{piv}}\bigr)^{\alpha-\beta} \exp{\beta-\alpha}\bigl(\frac{E}{E_{piv}}\bigr)^\beta, \text{ otherwise } \end{cases}\)

where

• \(A\) is the amplitude in ph/s/cm^2/keV

• \(E_{peak}\) is the \(\nu F_{\nu}\) peak in keV

• \(\alpha\) is the low-energy power-law index

• \(\beta\) is the high-energy power-law index

• \(E_{piv}\) is the pivot energy in keV. Nominally fixed to 100 keV.

---

BandOld

class gbm.spectra.functions.BandOld

Bases: gbm.spectra.functions.Function

The Band GRB function (a type of smoothly broken power law), using the original parameterization:

\(F(E) = A \begin{cases} \bigl(\frac{E}{E_{piv}}\bigr)^\alpha e^{-\frac{E}{E_0}}, \text{ if } {E < (\alpha-\beta) E_0}\\ \bigl(\frac{(\alpha-\beta) E_0}{E_{piv}}\bigr)^{(\alpha-\beta)} e^{(\beta-\alpha)} \bigl(\frac{E}{E_{piv}}\bigr)^\beta, \text{ otherwise } \end{cases}\)

where

• \(A\) is the amplitude in ph/s/cm^2/keV

• \(E_0\) is the break energy in keV

• \(\alpha\) is the low-energy power-law index

• \(\beta\) is the high-energy power-law index

• \(E_{piv}\) is the pivot energy in keV. Nominally fixed to 100 keV.

References

Band, D. et al. 1993, ApJ, 413, 281

---

BrokenPowerLaw

class gbm.spectra.functions.BrokenPowerLaw

Bases: gbm.spectra.functions.Function

A sharply-broken power law:

\(F(E) = A \begin{cases} \bigl(\frac{E}{E_{piv}}\bigr)^\alpha, & \text{ if } {E \leq E_b}\\ \bigl(\frac{E_b}{E_{piv}}\bigr)^\alpha \bigl(\frac{E}{E_b}\bigr)^\beta, & \text{ otherwise } \end{cases}\)

where

• \(A\) is the amplitude in ph/s/cm^2/keV

• \(E_{b}\) is the break energy in keV

• \(\alpha\) is the low-energy power-law index

• \(\beta\) is the high-energy power-law index

• \(E_{piv}\) is the pivot energy in keV. Nominally fixed to 100 keV.

---

DoubleBrokenPowerLaw

class gbm.spectra.functions.DoubleBrokenPowerLaw

Bases: gbm.spectra.functions.Function

A sharply-broken power law with two breaks:

\(F(E) = A \begin{cases} \bigl(\frac{E}{E_{piv}}\bigr)^\alpha, \text{ if } {E \leq E_{b1}}\\ \bigl(\frac{E_{b1}}{E_{piv}}\bigr)^\alpha \bigl(\frac{E}{E_b}\bigr)^\beta, \text{ if } {E > E_{b1} \text{ and } E \leq E_{b2}}\\ \bigl(\frac{E_{b1}}{E_{piv}}\bigr)^\alpha \bigl(\frac{E_{b1}}{E_{b2}}\bigr)^\beta \bigl(\frac{E}{E_{b2}}\bigr)^\gamma, \text{ otherwise} \end{cases}\)

where

• \(A\) is the amplitude in ph/s/cm^2/keV

• \(E_{b1}\) is the first break energy in keV

• \(E_{b2}\) is the second break energy in keV

• \(\alpha\) is the low-energy power-law index

• \(\beta\) is the mid-energy power-law index

• \(\gamma\) is the high-energy power-law index

• \(E_{piv}\) is the pivot energy in keV. Nominally fixed to 100 keV.

---

SmoothlyBrokenPowerLaw

class gbm.spectra.functions.SmoothlyBrokenPowerLaw

Bases: gbm.spectra.functions.Function

A smoothly-broken power law with a break scale:

\(F(E) = A \bigl(\frac{E}{E_{piv}}\bigr)^\beta 10^{(\beta-\beta_{piv})}, \\ \text{where }\\ \beta = m \Delta \ln\bigl(\frac{e^\alpha + e^{-\alpha}}{2}\bigr); \text{ } \beta_{piv} = m \Delta \ln\bigl(\frac{e^{\alpha_{piv}} + e^{-\alpha_{piv}}}{2}\bigr); \\ \alpha = \frac{\log(E/E_b)}{\Delta}; \text{ } \alpha_{piv} = \frac{\log(E_{piv}/E_b)}{\Delta}; \\ m = \frac{\lambda_2 - \lambda_1}{2}; \text{ } b = \frac{\lambda_1 + \lambda_2}{2};\)

and where

• \(A\) is the amplitude in ph/s/cm^2/keV

• \(E_{b}\) is the break energy in keV

• \(\Delta\) is the break scale in decades of energy

• \(\lambda_1\) is the low-energy power-law index

• \(\lambda_2\) is the high-energy power-law index

• \(E_{piv}\) is the pivot energy in keV. Nominally fixed to 100 keV.

---

LogNormal

class gbm.spectra.functions.LogNormal

Bases: gbm.spectra.functions.Function

A Log-Normal function:

\(F(E) = \frac{A}{\sqrt{2\pi} \sigma E} e^{-\bigl(\frac{\ln E - \mu}{2\sigma}\bigr)^2}\)

where

• \(A\) is the amplitude in ph/s/cm^2/keV

• \(\mu\) is the natural log of energy

• \(\sigma\) is the logarithmic standard deviation

---

GaussianLog

class gbm.spectra.functions.GaussianLog

Bases: gbm.spectra.functions.Function

A Gaussian with log_10(E):

\(F(E) = \frac{A}{\sqrt{2\pi} \sigma} e^{-\bigl(\frac{\log E-\log E_{cen} }{2\sigma}\bigr)^2}, \\ \text{ where } \\ \sigma = \frac{W}{2.35482}\)

and where

• \(A\) is the amplitude in ph/s/cm^2/keV

• \(E_{cen}\) is the centroid energy in keV

• \(W\) is the \(\log_{10}\) FWHM at \(E_{cen}\) in decades of energy

---

GaussianLogVaryingFWHM

class gbm.spectra.functions.GaussianLogVaryingFWHM

Bases: gbm.spectra.functions.Function

A Gaussian with log_10(E) and linearly-varying FWHM:

\(F(E) = \frac{A}{\sqrt{2\pi} \sigma} e^{-\bigl(\frac{\log E-\log E_{cen} }{2\sigma}\bigr)^2}, \\ \text{ where } \\ \sigma = \frac{W + s (\log E - \log E_{cen})}{2.35482}\)

and where

• \(A\) is the amplitude in ph/s/cm^2/keV

• \(E_{cen}\) is the centroid energy in keV

• \(W\) is the \(\log_{10}\) FWHM at \(E_{cen}\) in decades of energy

• \(s\) is the slope of \(W\) with respect to \(\log_{10}(E)\) in decades per \(\log_{10}(E)\)

---

SunyaevTitarchuk

class gbm.spectra.functions.SunyaevTitarchuk

Bases: gbm.spectra.functions.Function

Sunyaev-Titarchuk Comptonization spectra:

\(F(E) = A E^2 e^{-E} \int_0^\infty t^{n-1} e^{-t \bigl(1+\frac{t}{E}\bigr)^{n+3}} \text{d}t, \\ \text{ where } \\ n = -\frac{3}{2} + \sqrt{\gamma + \frac{9}{4}}; \\ \gamma = \frac{511 \pi^2}{G kT \bigl(\tau + \frac{2}{3} \bigr)^2}\)

and where

• \(A\) is the amplitude in ph/s/cm^2/keV

• \(kT\) is the electron energy in keV

• \(\tau\) is the optical depth

• \(G\) is the geometry factor, fixed at 3 for a spherical cloud and at 12 for a disk of electrons

References

Sunyaev, R. A. and Titarchuk, L. G. 1980, A&A, 86, 121

---

OTTB

class gbm.spectra.functions.OTTB

Bases: gbm.spectra.functions.Function

Optically-Thin Thermal Bremsstrahlung:

\(F(E) = A e^{-\frac{E}{kT}} e^{\frac{E_{piv}}{kT}} \biggl(\frac{E_{piv}}{E} \biggr)\)

where

• \(A\) is the amplitude in ph/s/cm^2/keV

• \(kT\) is the electron energy in keV

• \(E_{piv}\) is the pivot energy in keV

---

OTTS

class gbm.spectra.functions.OTTS

Bases: gbm.spectra.functions.Function

Optically-Thin Thermal Synchrotron:

\(F(E) = A e^{(E/E_C)^{1/3}}\)

where

• \(A\) is the amplitude in ph/s/cm^2/keV

• \(E_C\) is the energy scale in keV

References

Liang, E. P., et al., 1983, ApJ, 271, 776

---

BlackBody

class gbm.spectra.functions.BlackBody

Bases: gbm.spectra.functions.Function

Black Body spectrum:

\(F(E) = A \frac{E^2}{e^{E/kT}-1}\)

where

• \(A\) is the amplitude in ph/s/cm^2/keV

• \(kT\) is the temperature in keV

---

GaussLine

class gbm.spectra.functions.GaussLine

Bases: gbm.spectra.functions.Function

A Gaussian Line:

\(F(E) = A \frac{\text{erfc}(a_l)-\text{erfc}(a_r)}{2(E_r-E_l)}, \\ \text{ where } \\ a_l = \frac{E_l - E_{cen}}{\sqrt{2}\sigma}; \text{ } a_r = \frac{E_r - E_{cen}}{\sqrt{2}\sigma}; \\ \sigma = W/2.35482\)

and where

• \(A\) is the amplitude in ph/s/cm^2/keV

• \(E_{cen}\) is the centroid energy in keV

• \(W\) is the FHWM in keV

---

YangSoongPulsar

class gbm.spectra.functions.YangSoongPulsar

Bases: gbm.spectra.functions.Function

Yang Soong’s pulsar spectral form

\(F(E) = A C (1 - E_W G), \\ \text{ where} \\ C = \begin{cases} \bigl(\frac{E}{E_{piv}}\bigr)^{-\alpha}, & \text{ if } {E \leq E_b} \\ \frac{M}{E} e^{-E/E_{piv}}, & \text{ otherwise } \end{cases}; \\ M = e^{E_b/E_F} E_b \bigl(\frac{E_b}{E_{piv}}\bigr)^{-\alpha}; \\ G = \frac{0.94}{W} e^{-2.76 (\frac{E-E_{cen}}{W})^2};\)

and where

• \(A\) is the amplitude in ph/s/cm^2/keV

• \(\alpha\) is the power-law index

• \(E_b\) is the break energy in keV

• \(E_F\) is the e-folding energy in keV

• \(E_W\) is the equivalent line width in keV

• \(E_{cen}\) is the line centroid in keV

• \(W\) is the FWHM of the line, in keV

• \(E_{piv}\) is the pivot energy in keV

References

Soong, Y., et al. 1990, ApJ, 348, 641

---

TanakaPulsar

class gbm.spectra.functions.TanakaPulsar

Bases: gbm.spectra.functions.Function

Tanaka’s pulsar model

\(F(E) = A \frac{C}{C_{piv}} e^{L_{piv}-L}, \\ \text{ where} \\ C = A E^{-\alpha} e^{-E/kT}; \text{ } C_{piv} E_{piv}^{-\alpha} e^{-E_{piv}/kT}; \\ L = \sum_{i=1}^{N_L} \frac{\tau_i (i W)^2 [E/(i E_L)]^2}{(E-i E_L)^2 + (i W)^2}; \\ L_{piv} = \sum_{i=1}^{N_L} \frac{\tau_i (i W)^2 [E_{piv}/(i E_L)]^2} {(E_{piv}-i E_L)^2 + (i W)^2}; \\\)

and where

• \(A\) is the amplitude in ph/s/cm^2/keV

• \(\alpha\) is the power-law index

• \(kT\) is the temperature in keV

• \(\tau_1\) is the optical depth of the first line

• \(\tau_2\) is the optical depth of the second line

• \(N_L\) is the number of lines; fixed at either 0, 1, or 2

• \(E_L\) is the line centroid in keV of the first line

• \(W\) is the line FWHM in keV

• \(E_{piv}\) is the pivot energy in keV

References

Grove, J. E., et al. 1995, ApJ, 438, L25

---

LowEnergyCutoff

class gbm.spectra.functions.LowEnergyCutoff

Bases: gbm.spectra.functions.Function

A Low-Energy cutoff (Multiplicative Component):

\(F(E) = \begin{cases} (E/E_{cut})^{E_{cut}/E_F} e^{(E_{cut}-E)/E_F} & \text{ if } {E \leq E_{cut}} \\ 1 & \text{ otherwise } \end{cases}\)

where

• \(E_{cut}\) is the cutoff energy in keV

• \(E_{cen}\) is the folding energy in keV

---

HighEnergyCutoff

class gbm.spectra.functions.HighEnergyCutoff

Bases: gbm.spectra.functions.Function

A High-Energy cutoff (Multiplicative Component):

\(F(E) = \begin{cases} (E/E_{cut})^{E_{cut}/E_F} e^{(E_{cut}-E)/E_F} & \text{ if } {E > E_{cut}} \\ 1 & \text{ otherwise } \end{cases}\)

where

• \(E_{cut}\) is the cutoff energy in keV

• \(E_{cen}\) is the folding energy in keV

---

PowerLawMult

class gbm.spectra.functions.PowerLawMult

Bases: gbm.spectra.functions.Function

A Power Law (Multiplicative Component):

\(F(E) = \bigl(\frac{E}{E_{piv}} \bigr)^\lambda\)

where

• \(\lambda\) is the index

• \(E_{piv}\) is the pivot energy in keV

---

GaussLineMult

class gbm.spectra.functions.GaussLineMult

Bases: gbm.spectra.functions.Function

A Gaussian Line (Multiplicative Component):

\(F(E) = e^{I G}, \\ \text{ where } \\ G = \frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{(E-E_{cen})^2}{2\sigma^2}}; \\ \sigma = W/2.35482\)

and where

• \(I\) is the intensity

• \(E_{cen}\) is the centroid energy in keV

• \(W\) is the FWHM in keV

---

LorentzLineMult

class gbm.spectra.functions.LorentzLineMult

Bases: gbm.spectra.functions.Function

A Lorentzian Line (Multiplicative Component):

\(F(E) = e^{\frac{I}{W^2 + (E-E_{cen})^2}}\)

where

• \(I\) is the intensity

• \(E_{cen}\) is the centroid energy in keV

• \(W\) is the FWHM in keV

---