On the use of Cstat in testing models for Xray spectra
^{1} SRON Netherlands Institute for Space Research, Sorbonnelaan 2, 3584 CA Utrecht, The Netherlands
email: J.S.Kaastra@sron.nl
^{2} Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands
^{3} Department of Physics and Astronomy, Universiteit Utrecht, PO Box 80000, 3508 TA Utrecht, The Netherlands
Received: 15 July 2016
Accepted: 22 July 2017
Context. It has been shown that for the analysis of Xray spectra the Cstatistic, contrary to the χ^{2}statistic, provides unbiased estimates of the model parameters and their uncertainty ranges.
Aims. However, it is often stated that the Cstatistic cannot be used to carry out statistical tests on the goodness of fit of the model, and therefore several investigations are still based on χ^{2}statistics.
Methods. Here we show that it is straightforward to calculate the expected value and variance of the Cstatistic so that it can be used in tests.
Results. We provide formulae and simple numerical approximations to evaluate these expected values and variances. We also give examples indicating that tests based on only the expected value and variance of the Cstatistic are reliable for spectra even with only ~30 counts.
Conclusions. The Cstatistic can be used for statistical tests such as assessing the goodness of fit of a spectral model.
Key words: instrumentation: spectrographs / methods: data analysis / methods: statistical / Xrays: general
© ESO, 2017
1. Introduction
Xray spectra of astrophysical sources are often characterised by relatively low numbers of counts per spectral bin. In the early days, spectral models were often tested using χ^{2}statistics. The goodness of fit is expressed as (1)where the summation over i is over all n bins of the spectrum, N_{i} is the observed number of counts, s_{i} is the expected number of counts for the tested model, and for Poissonian statistics. There are three important remarks to make here.
First, both the model s_{i} and the observed spectrum N_{i} should include the source plus background counts to properly use Poissonian statistics.
Secondly, minimisation of χ^{2} to obtain the bestfit parameters of the model is easier when is approximated by N_{i}, which is a reasonable approximation when N_{i} is large and the Poissonian distribution approaches a normal distribution, but it fails for small N_{i}, which can be easily seen by putting σ_{i} = 0 in (1). This method leads to biased results, even for higher count rates (e.g. Nousek & Shue 1989; Mighell 1999). There are methods to compensate for this, however all of these have some drawbacks. For example the Gehrels (1986) approximation for N_{i} in the denominator of (1) causes problems when spectra are rebinned. Using (e.g. Wheaton et al. 1995) properly requires multiple passes through the minimisation procedure, updating at each step and sometimes leading to oscillatory behaviour.
Finally, it is often common practice in case of small N_{i} to rebin the spectra to at least 25 counts per bin or a similar sized number. This has the risk of washing out spectral details.
It has been pointed out by Cash (1979) that (2)is a much better statistic and can be applied to bins with a small number of counts without any bias in the derived parameters. Also, this statistic can be used to derive uncertainty ranges on the parameters of the model. It is often stated that the Cashstatistic cannot be used to measure the goodness of the fit using a quantity corresponding for instance to the reduced where p is the number of free parameters. For example, Nousek & Shue (1989) noted, “The principal disadvantage of the C statistic is that there is no value corresponding to the reduced χ^{2} value with which we can measure the goodness of the fit”, and Humphrey et al. (2009) wrote, “Since the absolute value of the Cstatistic cannot be directly interpreted as a goodnessoffit indicator observers typically prefer instead to minimize the betterknown χ^{2}fit statistic”. Maybe because of this many Xray astronomers keep using χ^{2}statistics, even in cases where the approximation breaks down or leads to biased parameters.
In this paper, I present calculations that can be used to evaluate the expected value for C and its rootmeansquared (rms) deviation to assess the goodness of fit of a spectral model derived from a best fit to the observed spectrum.
I use a modification of the original Cashstatistic that has been attributed to Castor and is implemented in current fitting packages such as XSPEC (Arnaud 1996), SHERPA (Freeman et al. 2001), and SPEX (Kaastra et al. 1996). This modified Cstatistic, designated here as cstat, is defined as (3)and has similar properties as the original Cash statistic, but in addition it can be used to assign a goodnessoffit measure to the fit. For a spectrum with many counts per bin C → χ^{2}, but where the predicted number of counts per bin is small, the expected value for C can be substantially smaller than the number of bins n.
2. Expected value of C and its variance
Fig. 1 Expected value of the contribution per bin to C, and its rms uncertainty as a function of the mean expected number of counts μ. 

Open with DEXTER 
The expected contribution C_{e,i} to the total C from any individual bin i and its variance C_{v,i} are given by with P_{k}(μ) the Poisson distribution (7)and μ the expected number of counts in the relevant bin. We show C_{e,i} and in Fig. 1.
The above equations hold for a single bin. For the full spectrum, the expected values C_{e,i} and variances C_{v,i} can be simply added over all bins i, yielding the expected value and variance for the full spectrum.
While the above procedure yields the exact expected mean value and variance for C, in general this does not mean that C has a Gaussian distribution with that mean and variance. A Gaussian distribution occurs when each bin, or most bins, i, have enough counts such that the distribution of C_{i} becomes asymptotically Gaussian, or for lower numbers of counts when the number of spectral bins is large enough. In the latter case we can use the central limit theorem, which states that when independent random variables are added, their sum tends towards a normal distribution even if the original variables themselves are not normally distributed.
In the above case, acceptable spectral models typically have ΣC_{e,i}(μ_{i})−f[ΣC_{v,i}(μ_{i})]^{0.5}<C< ΣC_{e,i}(μ_{i}) + f[ΣC_{v,i}(μ_{i})]^{0.5} with f a factor of order unity corresponding to the required significance level, for instance f = 1 for 68% confidence.
When these above conditions are not met, i.e. for low number of counts or low number of spectral bins, the distribution of C is not Gaussian. In that case the higher order moments of the distribution of C can be calculated analogous to (4) and (6). These can be used in principle to build the distribution of C. This can be a cumbersome task, however, and alternatively, using the bestfit model, one may test it simply by running multiple simulations of the spectrum to obtain an empirical distribution of C from which the goodness of fit can be estimated.
Fortunately, in most practical cases using the total mean and variance of C with a simple Gaussian approximation is accurate enough to assess the goodness of fit. We illustrate this with two practical examples in Sect. 3.
We implemented the above approach in the SPEX package^{1} (Kaastra et al. 1996). To help the user to see if a Cvalue corresponds to an acceptable fit, SPEX gives, after spectral fitting, the expected value of C and its rms spread, based on the bestfit model. Both quantities are simply determined by adding the expected contributions and their variances over all bins.
3. Simple approximations for the expected value and variance of C
We obtained simple approximations to the infinite series involved in (4) and (6), with relative errors better than 2.2 × 10^{4} for C_{e} and better than 1.6 × 10^{4} for C_{v}, as follows: With the help of the above equations, the goodness of fit for the model can be easily assessed. Given the other properties of cstat, such as unbiased parameter estimates, in almost all circumstances cstat is the preferred statistic to be used and the use of χ^{2}statistics in Xray spectral analysis should be avoided.
4. Two practical examples
Fig. 2 Simplified spectra of the Perseus cluster (top) with Hitomi and of Capella (bottom) with RGS. 

Open with DEXTER 
We tested our method to assess the goodness of fit with two examples. In the first example, we simulated a spectrum that has approximately the shape of the Perseus cluster as measured with the Hitomi SXS instrument (Hitomi Collaboration et al. 2016). For this demonstration purpose we simplified the model by adopting an isothermal spectrum in collisional ionisation equilibrium with a temperature of 4 keV, protosolar abundances, and an emission measure matching the flux of Perseus as measured by Hitomi. In our spectral fits, we used the temperature and emission measure of the source as free parameters. The spectrum has 5804 bins and is shown in Fig. 2, without noise, but with instrumental features included.
We scaled this spectrum in flux by factors of 10^{− k/ 2} with k ranging from 0–10, i.e. higher values of k corresponding to lower fluxes. For each flux value, we simulated 1000 spectra, performed a best fit, and determined C. We then produced a histogram of Cvalues and calculated the 90%, 95%, and 99% percentile points C_{90}, C_{95}, and C_{99} of this distribution. In an actual analysis of an observed spectrum, a model would be rejected at the 90% confidence level if C>C_{90}, and this is similar for the other confidence levels. We compare these percentile points, scaled with the expected mean C_{e} and variance C_{v}, as delivered by SPEX and described in Sect. 2, in Fig. 3.
Our second example is a simulation of the spectrum of Capella with the Reflection Grating Spectrometer (RGS) of XMMNewton. This spectrum is approximated by a single isothermal component with temperature 0.5 keV and appropriate flux for Capella. All further steps are the same as for the Perseus example. The main difference between both spectra is that while Perseus is dominated by continuum emission, owing to its high temperature, Capella is dominated by line emission, owing to its low temperature. The Capella spectrum has 1433 bins and is also shown in Fig. 2.
Fig. 3 Percentile points 90%, 95%, and 99% for the distribution of C for simulated Perseus (top) and Capella (bottom) spectra. The expected value C_{e} was subtracted and the difference is scaled with . The percentile points are shown as a function of the average number of counts in the spectra. The dotted lines show the expected percentile values if the distribution of C had been normal. These values are reached asymptotically for large numbers of counts in the spectra. 

Open with DEXTER 
From Fig. 3 we see that for more than about 30 counts in the spectrum the percentile points C_{90} and C_{95} (dots connected by solid lines) are close to the values calculated from the mean value and variance of C using a Gaussian approximation (dotted lines). This holds for both examples. For C_{95}, the Gaussian approximation even works for spectra with 10 counts.
But even going down to about 5 counts does not result in dramatic differences. In general one needs to be very cautious with spectra that have so few counts. For instance, 30 counts correspond to an average number of counts per bin of 0.005 and 0.02 only for the scaled Perseus and Capella spectra, respectively.
5. Conclusion
The Cstatistic can be used to estimate the goodness of fit of a model in the vast majority of all cases. When the spectrum has more than 10–30 counts, the distribution of C for the tested model is close enough to a Gaussian distribution for the highest confidence levels above 95%. This paper describes an algorithm that calculates the expected mean and variance of C that can be used to assess these confidence levels using a simple normal distribution.
Note added in proof. Equation (3), attributed here to Castor, was introduced also by Baker & Cousins (1984).
Acknowledgments
SRON is supported financially by NWO, The Netherlands Organization for Scientific Research. I thank the referee for useful comments.
References
 Arnaud, K. A. 1996, in Astronomical Data Analysis Software and Systems V, eds. G. H. Jacoby, & J. Barnes, ASP Conf. Ser., 101, 17 [Google Scholar]
 Baker, S., & Cousins, R. D. 1984, Nucl. Inst. Meth. Phys. Res., 221, 437 [NASA ADS] [CrossRef] [Google Scholar]
 Cash, W. 1979, ApJ, 228, 939 [NASA ADS] [CrossRef] [Google Scholar]
 Freeman, P., Doe, S., & Siemiginowska, A. 2001, in Astronomical Data Analysis, eds. J.L. Starck, & F. D. Murtagh, Proc. SPIE, 4477, 76 [Google Scholar]
 Gehrels, N. 1986, ApJ, 303, 336 [NASA ADS] [CrossRef] [Google Scholar]
 Hitomi Collaboration, Aharonian, F., Akamatsu, H., et al. 2016, Nature, 535, 117 [NASA ADS] [CrossRef] [Google Scholar]
 Kaastra, J. S., Mewe, R., & Nieuwenhuijzen, H. 1996, in UV and Xray Spectroscopy of Astrophysical and Laboratory Plasmas, eds. K. Yamashita, & T. Watanabe, 411 [Google Scholar]
 Mighell, K. J. 1999, ApJ, 518, 380 [NASA ADS] [CrossRef] [Google Scholar]
 Nousek, J. A., & Shue, D. R. 1989, ApJ, 342, 1207 [NASA ADS] [CrossRef] [Google Scholar]
 Wheaton, W. A., Dunklee, A. L., Jacobsen, A. S., et al. 1995, ApJ, 438, 322 [NASA ADS] [CrossRef] [Google Scholar]
All Figures
Fig. 1 Expected value of the contribution per bin to C, and its rms uncertainty as a function of the mean expected number of counts μ. 

Open with DEXTER  
In the text 
Fig. 2 Simplified spectra of the Perseus cluster (top) with Hitomi and of Capella (bottom) with RGS. 

Open with DEXTER  
In the text 
Fig. 3 Percentile points 90%, 95%, and 99% for the distribution of C for simulated Perseus (top) and Capella (bottom) spectra. The expected value C_{e} was subtracted and the difference is scaled with . The percentile points are shown as a function of the average number of counts in the spectra. The dotted lines show the expected percentile values if the distribution of C had been normal. These values are reached asymptotically for large numbers of counts in the spectra. 

Open with DEXTER  
In the text 