MULTIPARAMETER LINEAR LEAST-SQUARES FITTING TO POISSON DATA ONE COUNT AT A TIME

A standard problem in gamma-ray astronomy data analysis is the decomposition of a set of observed counts, described by Poisson statistics, according to a given multicomponent linear model, with underlying physical count rates or fluxes which are to be estimated from the data. Despite its conceptual simplicity, the linear least-squares (LLSQ) method for solving this problem has generally been limited to situations in which the number n(i) of counts in each bin i is not too small, conventionally more than 5-30. It seems to be widely believed that the failure of the LLSQ method for small counts is due to the failure of the Poisson distribution to be even approximately normal for small numbers. The cause is more accurately the strong anticorrelation between the data and the weights w(i) in the weighted LLSQ method when square-rootn(i) instead of square-rootnBAR(i) is used to approximate the uncertainties, sigma(i), in the data, where nBAR(i) = E[n(i)], the expected value of n(i). We show in an appendix that, avoiding this approximation, the correct equations for the Poisson LLSQ (PLLSQ) problem are actually identical to those for the maximum likelihood estimate using the exact Poisson distribution. Since weighted linear least-squares involves a kind of weighted averaging, LLSQ estimators generally produce biased results when the data and their weights are correlated. We describe a class of weighted PLLSQ estimators which are linear functions of the observed counts. Such PLLSQ estimators are unbiased independent of nBAR(i), even when the average number of counts in an entire fit is much less than one. Their variance is a minimum when the weights are calculated from the true variances of the data, but in general these are not accurately known. Fortunately, the variance of the estimate is a very weak function of the weights near the optimum value, so for the PLLSQ problem it is easy in practice to find weights that are virtually ideal yet still completely unbiased. PLLSQ estimators which are linear in the data also allow fitting multiple data sets by the calculation of only a scalar product, without the need to repeat the accumulation and solution of the LLSQ equations. Owing also to the linearity of the estimates in the data, each count contributes to the answers independently of every other count, so that the results for small bins are independent of the particular choice of binning. This property makes possible PLLSQ methods which avoid binning the data altogether. Some alternatives to the approximation of the uncertainties in the data by the square root of the observed counts are discussed. We apply the method to solve a problem in high-resolution gamma-ray spectroscopy for the JPL High-Resolution Gamma-Ray Spectrometer flown on HEAO 3. Systematic error in subtracting the strong, highly variable background encountered in the low-energy gamma-ray region can be significantly reduced by closely pairing source and background data in short segments. Significant results can be built up by weighted averaging of the net fluxes obtained from the subtraction of many individual source/background pairs. Extension of the approach to complex situations, with multiple cosmic sources and realistic background parameterizations, requires a means of efficiently fitting to data from single scans in the narrow (almost-equal-to 1.2 keV, for HEAO 3) energy channels of a Ge spectrometer, where the expected number of counts obtained per scan may be very low. Such an analysis system is discussed and compared to the method previously used.