ERROR_PROPAGATION · IMPLICATIONS_FOR_CHISQ_CALCULATIONS,_AND_BACKWARD_COMPATIBILITY
errors -- description of PROS error computations for low count/bin data
This document describes how errors are computed and propagated in PROS tasks which deal with binned data. Particular attention is given to the case where the number of events per bin is small and Poisson rather than Gaussian statistics apply.
Many tasks in X-ray Astronomy require the analysis of x-ray events which have been binned in either time, energy, or space. The most obvious example is fitting a model spectrum to a background-subtracted pulse-height spectrum, but other examples, such as determining the net count rate and error in spatial region(s), or searching for periodicities by examining pulse profiles binned at different trial periods, are also common. In such tasks, an accurate estimate of the errors associated with the binned data is often crucial for the proper interpretation of the results.
The appropriate probability distribution for the binned data is, in general, the Poisson distribution. However, as the number of counts per bin, n, becomes large, the distribution approaches a Gaussian distribution with sigma = sqrt(n). Exactly what "becomes large" means depends, to some extent, on context, but typical values of n ~ 20 are often cited as sufficient for the use of the Gaussian limit. In fact, since n ~ 20 is often easily achieved in x-ray observations, the assumption of a Gaussian distribution with sigma = sqrt(n) is generally used for all n, with the warning that the user should adjust bin sizes to avoid n < ~ 20. This was the approach in PROS prior to version 2.3.
However, beginning with that version, we have modified our approach in recognition of the many cases when users are unable or simply unwilling to avoid small numbers of counts per bin by rebinning. Uncritical use of the sqrt(n) approximation in such cases may lead to problems. The case when n = 0, for example, is a particular problem, since the sqrt(n) approximation would imply 0 error, which is clearly incorrect. Moreover, since many procedures, such as CHISQ fitting, require the inverse of sigma**2, bins with n = 0 would have to be eliminated, in order for such procedures to run properly. More importantly, cumulative probabilities of the Poisson distribution between limits defined in terms of sigma = sqrt(n) differ significantly from the familiar values appropriate to the Gaussian distribution when n is small. For example, n + sqrt(n), which yields the 1 sigma (84% CL) single-sided upper limit in the Gaussian limit, significantly underestimates this quantity for n < ~ 20. Thus, proper quantitative interpretation of errors requires improved estimates for errors for small n.
In the remainder of this document, we discuss such an estimate, which approaches the Gaussian limit for large n, but which has the advantage of yielding a well-defined and consistent upper limit for all n. Since an important use of binned counts is in parameter estimation via CHISQ minimization, we end with a brief discussion of the impact of our approach on that technique.
The inadequacy of the sqrt(n) approximation for small n has long been known (Regener, 1951). A detailed treatment of the problem was given by Gehrels (1986), and we model our approach after his.
If n counts per bin are actually observed, then the upper limit U, to the true value per bin N (i.e. the parent mean) at the confidence level C is given by the solution to
n __ \ / U**i * exp(-U)/(i!) = 1 - C (1) -- i=0
for Poisson statistics. Note that U is well-defined when n=0. Values of U for various confidence levels appear in Table 1 of Gehrels (1986). The confidence level we're interested in is C = 0.8413, since this corresponds to the Gaussian 1 sigma limit. It's equally possible to deal with lower limits, but as it turns out, the error derived from the upper limit analysis in the Poisson case is the larger of the two, and hence the conservative choice. The error is then estimated to be "sigma" = U - n. In the Gaussian limit this just reduces to sqrt(n). Values of U, sqrt(n), and sigma are listed below, using Table 1 of Gehrels, for various values of n. Also listed are values to an approximate solution to eq. 1 for C=0.8413 (cf. Gehrels eq. 7). The approximation is
sigma' ~ 1 + sqrt(n + 0.75) (2)
and, as one can see from the %diff = (sigma'-sigma)/sigma, sigma' differs from sigma by less than 2% overall.
Table 1 n U sqrt(n) sigma sigma' %diff 0 1.841 0 1.841 1.866 1.4 1 3.300 1 2.300 2.323 1.0 2 4.638 1.414 2.638 2.658 0.8 5 8.382 2.236 3.382 3.398 0.5 10 14.27 3.162 4.27 4.28 - 20 25.55 4.472 5.55 5.56 - 40 47.38 6.325 7.38 7.38 - 100 111.0 10 11.0 11.0 -
Because eq. 2 is so easy and accurate, we have adopted it whenever we need to compute the "1 sigma error" on n counts.
It should be noted, however, that in citing a single error, we are still implicitly assuming a symmetric distribution, and thus continue to offer only an approximate solution. Our intent is not to provide a rigorous treatment but to provide a single error estimate which reduces to the traditional sqrt(n) limit for large n while providing a more consistent estimate of a particular statistic for small n. Users who require more rigorous treatments are encouraged to consider other, more sophisticated analysis techniques which can deal with unbinned data and the exact (Poisson) probability distribution.
The sigma' approximation is used only to calculate errors on raw counts. Calculation of errors on derived quantities, such as net count rate, is still carried out using standard error propagation formulae. For example, if net count rate in a region is given by
NET = [ T - B * Area(T)/Area(B) ] / TIME,
SIGMA(NET) = SQRT[SIGMA'(T)**2 + SIGMA'(B)**2 * (Area(T)/Area(B))**2]/TIME.
A number of tasks contain a "no-background" option, i.e., an option to zero the background value when calculating net quantities. Since the error on 0 counts is no longer 0, it is assumed that if the user selects such an option, then the background error should also be zeroed.
A quick examination of Table 1 reveals that sigma' is greater than sqrt(n) for all n. Since CHISQ depends on 1/sigma'**2, one should expect smaller CHISQ values in general. The effect should reduce CHISQ to ~ 65% of the previous value for n ~ 20 on average. Analysis of the standard test cases in the XTIMING package which compute CHISQ confirm the expected result but also demonstrate that the actual CHISQ minima, and hence best-fit parameter values, are not significantly affected, for reasonable values of n.
To further investigate the effect of the sigma' approximation on CHISQ, best-fit parameters, and their errors, we have performed a simple numerical experiment. We have generated sets of 10,000 random Poisson deviates for parent means of 5, 10, 20, 40, and 100, respectively, and divided each set into 400 samples with 25 members each. For each sample, we have calculated the sample mean, according to
25 __ \ / x(i)/sigma(i)**2 -- i=1 <X> = ---------------------- (3) 25 __ \ / 1/sigma(i)**2 -- i=1
which yields the minimum CHISQ for a fit to a constant. We have also calculated the minimum CHISQ and the deviation, DX, from <X> which increases CHISQ by 1 (corresponding to a 68% confidence interval for one interesting parameter). We have repeated this process using 3 different techniques for calculating sigma(i), namely, sigma(i)=sqrt(parent mean), sigma(i)=sqrt(x(i)), and sigma(i)=sigma'(i). For each technique, we tabulate below <<X>>, the average of the sample means, <CHISQ>, the average minimum CHISQ, and the percentage P68 of samples for which <X>+/-DX includes the parent mean.
Table 2 Sigma ------------------------------------------------------------------- sqrt(Parent Mean) sqrt(x(i)) sigma' Parent -------------------- -------------------- -------------------- Mean <<X>> <CHISQ> P68 <<X>> <CHISQ> P68 <<X>> <CHISQ> P68 5 5.01 23.67 71 3.88 28.95 12 4.55 10.91 52 10 9.99 23.35 72 8.93 26.32 28 9.28 13.30 55 20 20.11 24.50 68 19.08 25.71 44 19.31 16.27 65 40 39.98 24.49 69 38.97 25.22 53 39.12 18.34 64 100 100.00 24.21 68 99.02 24.45 60 99.12 20.03 67
These results indicate that the sigma' approximation provides a reasonable estimate of the 68% confidence interval for n > ~ 20, and a much better estimate than the sqrt(n) approximation for all n. The relatively poor performance of the latter in this respect is due, in large part, to the significant bias evident in calculating <X>. This bias has long been known (see, e.g., Bevington, 1969, pp 248-250) and has been recently demonstrated in a context specific to X-ray Astronomy (Nousek and Shue, 1989). It is due to an error estimate, sqrt(x(i)), which is correlated with the data, x(i) (Wheaton, et al. 1993). The sigma' approximation also suffers from this bias, although to a lesser extent; this, combined with the somewhat larger errors inherent in the latter approximation, result in its improved performance. Users should note, however, that the sigma=sqrt(Parent Mean) technique provides unbiased estimates of <X> and excellent estimates of the confidence interval for all values of n. A generalization of this approach for multi-parameter fits is beyond the scope of this document. Users are referred to Wheaton et al. (1993).
Finally, it is recognized that many users may still wish to use the sqrt(n) approximation, either to compare with earlier analysis results, to assess the effect of the sigma' approximation, or simply because they feel more comfortable with the sqrt(n) approximation. A technique is thus provided for forcing the use of this approximation. An additional keyword, POISSERR, has been added to QPOE and image headers, and defaults to TRUE for the sigma' approximation. Changing the value of this keyword to FALSE using qphedit (for QPOE files) or hedit (for image files) will force the use of the sqrt(n) approximation in all later error calculations using those files.
The following procedure will change the QPOE header keyword POISSERR from T (true for POISSON errors) to F (false for GAUSSIAN errors). a) verify current value > imhead rh110267n00.qp long+ | match POISSERR POISSERR= T b) change T to F (Note the value is boolean, so DO NOT use quotes) > qphedit rh110267n00.qp POISSERR F rh110267n00.qp,POISSERR (T -> F): y rh110267n00.qp,POISSERR: T -> F continue ? (yes): rh110267n00.qp updated b) verify update > imhead rh110267n00.qp long+ | match POISSERR POISSERR= F i.e. the same procedure using 'hedit' is valid for images.
Bevington, P.R. 1969, "Data Reduction and Error Analysis for the Physical Sciences",(McGraw-Hill:New York).
Gehrels, N. 1986, ApJ, 303, 336.
Nousek, J.A. & Shue, D.R. 1989, ApJ, 342, 1207.
Regener, V. 1951, Phys. Rev., 84, L161.
Wheaton, W.A., Dunklee, A.L., Jocobson, A.S., Ling, J.C., Mahoney, W.A. & Radocinski, R.G. 1993, submitted to ApJ.