Calculating errors on flux from model components

How xspec calculates flux errors

Flux is not in general a parameter in a model so we cannot estimate a confidence interval in the usual fashion using the error command, which works by varying the parameter value and looking for a change in the fit statistic by a set amount (eg 2.706 for 90% confidence). Suppose instead we consider the 90% confidence region of parameter space for all interesting parameters. For instance, if there are two interesting parameters this would be the region inside the 90% confidence contours produced after a steppar and plot contour. An estimate of the 90% confidence interval on the flux would be to calculate the flux for all sets of parameters which lie within this 90% confidence region. However, this is not practical in general because with a large number of parameters we would have to systematically step through a high-dimension grid, which is likely to be computationally expensive.

So, instead we use a method of selecting sets of parameters which will draw them out of a probability distribution which should approximate the confidence contours. We can then calculate the flux for each set of parameters, order the fluxes and determine the range into which 90% of the values fall. This is then an estimate for the 90% confidence region. So, it just remains to find a way of selecting appropriate sets of parameters. We use the covariance matrix. At the best-fit position the first derivatives of the fit statistic with respect to the parameters are by definition zero. The matrix of second derivatives is called the covariance matrix (and its inverse is often referred to as the Fisher information matrix). We draw sets of parameters from a multi-dimensional Normal distribution whose mean is given by the best-fit and whose sigma is the covariance matrix. It is important to note that this is an approximation and if the fit statistic depends in a complicated fashion on the parameter values it can be very wrong. We are developing a better algorithm based on the use of Markov Chain Monte Carlo.

Calculating the error on the flux from part of the model

A common extension to estimating the flux and its confidence region for the whole model is to do the same for some part of the model. For instance, we might want to know the flux and error without the effects of Galactic absorption. Unfortunately, this is not as easy. We can set the Galactic absorption column to zero and calculate the flux however the covariance matrix we stored from the best-fit is now not appropriate. If we continue as before then the sets of parameters generated will have non-zero absorption columns and hence the fluxes are not correct. (Note that xspec v11 will actually appear to run in this case and return an error but the value is meaningless). What we need to do is to generate sets of parameters from the best-fit values and the covariance matrix and then within each set of parameters set the absorption column to zero and calculate the flux. These fluxes can then be ordered and the 90% confidence region found. Unfortunately, we cannot do this within xspec at the moment. However, there is a way of writing a Tcl script to do this calculation external to xspec. There is a tclout option, simpar, which returns a vector of parameters based on the best-fit values and the covariance matrix.

As an example we have supplied a script which can be used to calculate the 90% confidence interval on the unabsorbed flux ie the flux for the model without an initial Galactic absorption component. To use this script download it then within xspec :

XSPEC12>source fluxerror.tcl
XSPEC12>fluxerror 500 {1} 90

This runs 500 realizations of the parameters, sets the first parameter to zero (ie switch off the Galactic absorption), then calculate the flux. The {} is a Tcl list of parameters to set to zero so other components can be switched off by zeroing their normalizations.

Xspec Home Page

Last modified: Monday, 07-May-2007 12:09:38 EDT