Introduction

There are two operations performed in XSPEC that require statistics. The first is parameter estimation, which comprises finding the parameters for a given model that provide the best fit to the data and then estimating uncertainties on these parameters. The second operation is testing whether the model and its best-fit parameters actually match the data. This is usually referred to as determining the goodness-of-fit.

Which statistics should be used for these two operations depends on the probability distributions underlying the data. Almost all astronomical data are drawn from one of two distributions: Gaussian (or normal) and Poisson. The Poisson distribution is the familiar case of counting statistics and is valid whenever the only source of experimental noise is due to the number of events arriving at the detector. This is a good approximation for modern CCD instruments. If some other sort of noise is dominant then it is usually described by the Gaussian distribution. A common example of this is detectors that require background to be modeled in some way, rather than directly measured. The uncertainty in the background modeling is assumed to be Gaussian.

In the limit of large numbers of counts the Poisson distribution can be well approximated by a Gaussian so the latter is often used for detectors with high counting rates. In most cases this will cause no errors and does simplify the handling of background uncertainties however care should be exercised that no systematic offsets are introduced.

A fuller discussion of many of the issues discussed in this appendix can be found in Siemiginowska (2011).

Parameter Estimation

The standard statistic used in parameter estimation is the maximum likelihood. This is based on the intuitive idea that the best values of the parameters are those that maximize the probability of the observed data given the model. The likelihood is defined as the total probability of observing the data given the model and current parameters. In practice, the statistic used is twice the negative log likelihood.

**For Gaussian data (chi)**

The likelihood for Gaussian data is

_{}

where *y _{i}* are the
observed data rates,

_{}

commonly referred to as χ^{2}
and used for the **statistic chi** option.

**For Gaussian data with
background (chi)**

The previous section assumed that
the only contribution to the observed data was from the model. In practice,
there is usually background. This can either be included in the model or taken
from another spectrum file (read in using the **back** command). In the
latter case the *y _{i }*become observed data rates from the source
spectrum subtracted by the background spectrum and the

For Poisson data (cstat)

The likelihood for Poisson distributed data is:

_{}

where *S _{i}* are the
observed counts,

_{}

Castor (priv. comm) has pointed out that modifying this by a quantity that depends only on the data (and hence makes no difference to the best-fit parameters) to give :

_{}

provides a statistic which
asymptotes to *S ^{2 }*in the limit of large number of counts. This
is what is used for the

For Poisson data with Poisson background (cstat)

This case is more difficult than that of Gaussian data because the difference between two Poisson variables is not another Poisson variable so the background data cannot be subtracted from the source and used within the C statistic. The combined likelihood for the source and background observations can be written as:

_{}

where *t _{s}*
and

As a simple illustration suppose the source spectrum is source.pha and the background spectrum back.pha. The source model is an absorbed apec and the background model is a power-law. Further suppose that the background model requires a different response matrix to the source, backmod.rsp say. The fit is set up by:

XSPEC12> data 1:1 source.pha 2:2 background.pha

XSPEC12> resp 2:1 backmod.rsp 2:2 backmod.rsp

XSPEC12> model 2:backmodel pow

where the normalization of the apec model is fixed to zero for the second data group (i.e. the background spectrum) and the parameters of the background model are linked between the data groups.

If there
is no appropriate model for the background it is still possible to proceed.
Suppose that each bin in the background spectrum is given its own parameter so
that the background model is *b _{i}* =

_{}

where

_{}

and

_{}

The sign
of *d _{i}* in

_{}

If *B _{i}*

_{}

otherwise:

_{}

This W
statistic is used for **statistic cstat** if a background spectrum with
Poisson statistics has been read in. In practice, it works well for many cases
but for weak sources can generate an obviously wrong best fit. It is not clear
why this happens although binning to ensure that every bin contains at least
one count often seems to fix the problem.

In the
limit of large numbers of counts per spectrum bin a second-order Taylor
expansion shows that *W* tends to :

_{}

which is
distributed as χ^{2 }with *N – M* degrees of freedom, where
the model *m _{i}* has M parameters (include the normalization).

For Poisson data with Gaussian background (pgstat)

Another possible background option is if the background spectrum is not Poisson. For instance, it may have been generated by some model based on correlations between the background counts and spacecraft orbital position. In this case there may be an uncertainty associated with the background which is assumed to be Gaussian. In this case the same technique as above can be used to derive a profile likelihood statistic :

_{}

where

_{}

and

_{}

There is
a special case for any bin with *S _{i}*

_{}

This is
what is used for the **statistic pgstat** option.

Bayesian analysis of Poisson data with Poisson background (lstat)

An alternative approach to fitting Poisson data with background is to use Bayesian methods. In this case instead of solving for the background rate parameters we marginalize over them writing the joint probability distribution of the source parameters as :

_{}

where {θ_{j}}
are the source parameter, {b_{k}} the background rate parameters and I
any prior information. Using Bayes theorem, that the {θ_{j}} and
independent of the {b_{k}}, that the {b_{k}} are individually
independent and that the observed counts are Poisson gives :

_{}

where :

_{}

To
calculate J_k we need to make an assumption about the prior background
probability distribution, p(b_{k}|I). We follow Loredo (1992) and
assume a uniform prior between 0 and b_{i}^{max}. Expanding the
binomial gives :

_{}

where :

_{}

Again,
follow Loredo we assume that _{}and using the approximation
_{}when
_{}gives
:

_{}

Note that
for m_{k} = 0 only the j = 0 term in the summation is non-zero. Now, we
define lstat by calculating -2 ln P and ignoring all additive terms which are
independent of the model parameters :

_{}

**For
power spectra from time series data (whittle)**

XSPEC has been used by a number of researchers to fit models to power spectra from time series data. In this case the x-axis is frequency (in Hz) and not keV so plots have to be modified appropriately. The correct fit statistic is that due to Whittle as discussed in Vaughan (2010) and Barret & Vaughan (2012) :

_{}

** **

Parameter confidence regions

Fisher Matrix

XSPEC provides several different methods to estimate the precision with which parameters are determined. The simplest, and least reliable, is based on the inverse of the second derivative of the statistic with respect to the parameter at the best fit. The first derivative must be zero by construction and the second derivative provides a measure of how rapidly the statistic increases away from the best-fit. The faster the statistic increases, i.e. the larger the second derivative, the more precisely the parameter is determined. The matrix of second derivatives is often referred to as the Fisher information. Its inverse is the covariance matrix, written out at the end of an XSPEC fit.

The +/- numbers provided for each parameter in the standard fit output are estimates of the one-sigma uncertainty, calculated as the square root of the diagonal elements of the covariance matrix. As such, these ignore any correlations between parameters. Whether correlations are important can be seen by comparing with the off-diagonal elements of the covariance matrix. In general, these estimates should be considered lower limits to the true uncertainty.

Correlation information is also given in the table of variances and principal axes which also appears at the end of a fit. Each row in this table is an eigenvalue and associated eigenvector of the Fisher matrix. If the parameters are independent then each eigenvector will have a contribution from only one parameter. For instance, if there are three independent parameters then the eigenvectors will be (1,0,0), (0,1,0), and (0,0,1). If the parameters are not independent then each eigenvector will show contributions from more than one parameter.

Delta Statistic

The next most reliable method for deriving parameter confidence regions is to find surfaces of constant delta statistic from the best-fit value, i.e. where :

_{}

This is the method used by the **error** command, which
searches for the parameter value where the statistic differs from that at the
best fit by a value (Δ) specified in the command. For each value of the
parameter being tested all other free parameters are allowed to vary. The
results of the **error** command can be checked using **steppar**,** **which
can also be used to find simultaneous confidence regions of multiple
parameters. The specific values of Δ which generate particular confidence
regions are calculated by assuming that_{}is distributed as χ^{2
}with the number of degrees of freedom equal to the number of parameters
being tested (e.g. when using the **error** command there is one degree of
freedom, when using **steppar** for two parameters followed by **plot
contour** there are two degrees of freedom). This assumption is correct for
the *S ^{2}* statistic and is asymptotically correct for other
statistic choices.

Monte Carlo

The best but most computationally expensive methods for estimating parameter confidence regions are using two different Monte Carlo techniques. The first technique is to start with the best fit model and parameters and simulate datasets with identical properties (responses, exposure times, etc.) to those observed. For each simulation, perform a fit and record the best-fit parameters. The sets of best-fit parameters now map out the multi-dimensional probability distribution for the parameters assuming that the original best-fit parameters are the true ones. While this is unlikely to be true, the relative distribution should still be accurate so can be used to estimate confidence regions. There is no explicit command in XSPEC to use this technique however it is easy to construct scripts to perform the simulations and store the results.

The second technique is Markov Chain Monte Carlo (MCMC) and is
of much wider applicability. In MCMC a chain of sets of parameter values is
generated which describe the parameter probability distribution. This
determines both the best-fit (the mode) and the confidence regions. The **chain**
command runs MCMC chains which can be converted to probability distributions
using **margin **(which takes the same arguments as **steppar**). The
results can be plotted in 1- or 2-D using **plot margin **however this is
not quite as useful as it might be because what is plotted is the probability,
not the probability within some region. If MCMC chains are in use then the **error**
command will use them to estimate the parameter uncertainty.

** **

Goodness-of-fit

Parameter
values and confidence regions only mean anything if the model actual fits the
data. The standard way of assessing this is to perform a test to reject the
null hypothesis that the observed data are drawn from the model. Thus we
calculate some statistic T and if T_{obs} > T_{critical}
then we reject the model at the confidence level corresponding to T_{critical}.
Ideally, T_{critical} is independent of the model so all that is
required to evaluate the test is a table giving T_{critical }values for
different confidence levels. This is the case for χ^{2} which is
one of the reasons why it is used so widely. However, for other test statistics
this may not be true and the distribution of T must be estimated for the model
in use then the observed value compared to that distribution. This is done in
XSPEC using the **goodness** command. The model is simulated many times and
a value of T calculated for each fake dataset. These are then ordered and a
distribution constructed. This distribution can be plotted using **plot
goodness**. Now suppose that T_{obs} exceeds 90% of the simulated T
values, we can reject the model at 90% confidence.

It is worth emphasizing that goodness-of-fit testing only allows us to reject a model with a certain level of confidence, it never provides us with a probability that this is the correct model.

**Chi-square (chi)**

The
standard goodness-of-fit test for Gaussian data is χ^{2 }(as
defined above). At the end of a **fit**, XSPEC writes out the reduced χ^{2
}(= χ^{2}/ν , where ν is the number of degrees of
freedom = number of data bins – number of free parameters). A rough rule of
thumb is that the reduced χ^{2} should be approximately one. If
the reduced χ^{2 }is much greater than one then the observed data
are likely not drawn from the model. If the reduced χ^{2 }is much
less than one then the Gaussian sigma associated with the data are likely
over-estimated. XSPEC also writes out the null hypothesis probability, which is
the probability of the observed data being drawn from the model given the value
of χ^{2 }and^{ }the number of degrees of freedom.

**Pearson
chi-square (pchi)**

Pearson's original (1900) chi-square test was not for Gaussian data but for the case of dividing counts up between cells. This corresponds to the case of Poisson data with no background.

_{}

**Kolmogorov-Smirnov
(ks)**

There are a number of test statistics based on the empirical distribution function (EDF). The EDF is the cumulative spectrum :

_{}and _{}

The EDF
can be plotted using **plot icounts**. The best known of these tests is
Kolmogorov-Smirnov whose statistic is simply the largest difference between the
observed and model EDFs :

_{}

The XSPEC
**statistic test ks** option returns log D. The significance of the ks value
can be determined using the **goodness** command.

**Cramer-von
Mises (cvm)**

The Cramer-von Mises statistic is the sum of the squared differences of the EDFs :

_{}

The XSPEC
**statistic test cvm** option returns log w^{2} and its significance
should be determined using the **goodness** command.

**Anderson-Darling
(ad)**

Anderson-Darling is a modification of Cramer-von Mises which places more weight on the tails of distribution :

_{}

**Runs
(runs)**

The Runs
(or Wald-Wolfowitz) test checks that residuals are randomly distributed above
and below zero and do not cluster. Suppose N_{p }is the number of
channels with +ve residuals, N_{n} the number of channels with negative
residuals, and R the number of runs then the Runs statistic is :

_{}

where :

_{}and _{}

As for the EDF tests, XSPEC returns log Runs.

** **

References

Barret, D. & Vaughan, S., 2012. “Maximum likelihood fitting of X-ray power density spectra: application to high-frequency quasi-periodic oscillation from the neutron star X-ray binary 4U1608-522”, ApJ 746, 131.

Cash, W., 1979. “Parameter estimation in astronomy through application of the likelihood ratio”, ApJ 228, 939.

Loredo, T., 1992. In “Statistical Challenges in Modern Astronomy”, eds. Feigelson, E.D. and Babu, G.J., pp 275-297.

Siemiginowska, A., 2011. In “Handbook of X-ray Astronomy” eds. Arnaud, K.A., Smith, R.K. and Siemiginowska, A., Cambridge University Press, Cambridge.

Vaughan, S. 2010, “A Bayesian test for periodic signals in red noise”, MNRAS 402, 307.