Fast New EM-type Algorithms with Applications in Astrophysics
Submitted to The Astrophysical Journal
David A. van Dyk
Department of Statistics, Harvard University, Cambridge, MA
02138, U.S.A.
The Richardson-Lucy algorithm has long been a popular
tool for restoring degraded images obtained, for example, with
CCD detectors. This algorithm is actually a special case of
the EM algorithm, a well known method in the statistics literature,
which is popular for numerically stable model fitting while
accounting for missing data or latent model structure. The two-fold
objective of this paper is first to introduce the EM algorithm and
its many recent extensions and improvements and second to
illustrate a number of new astrophysical applications of
EM-type algorithms for robust chi square
fitting, (penalized) maximum likelihood estimation, and
the computation of Bayesian posterior modes. The recent extensions to
the EM algorithm aim to extend its wide range of applications and to
reduce required computational time. We show that these extensions
can significantly improve the convergence rate of an EM
algorithm for likelihood-based
inference in a new class of spectral models which allow for
high-resolution low-count Poisson data. These models will soon be
available in the
Chandra Interactive Analysis of Observations (CIAO) software
and account for absorption of photons, a
continuum, several line profiles, stochastic instrument response,
and background contamination. In particular, our method for
handling background counts is integrated into a model for
the observed counts, thus
avoiding the need for ad hoc solutions such as subtracting off counts
and the embarrassing resulting difficulty of negative counts. This spectral
model is especially useful with low-count per bin images as are
expected from the new generation of high-resolution space observatories
such as the Chandra X-Ray Observatory, XMM, Constellation X, GLAST,
etc. This application of EM not only demonstrates the power of EM to
handle highly
structured hierarchical models, but also illustrates how EM-type
algorithms can be used to solve outstanding problems in astrophysics.
Return to David van Dyk's homepage.