Abstract
We introduce a new galaxy image decomposition tool, galphat (GALaxy PHotometric ATtributes), which is a front-end application of the Bayesian Inference Engine (bie), a parallel Markov chain Monte Carlo package, to provide full posterior probability distributions and reliable confidence intervals for all model parameters. The bie relies on galphat to compute the likelihood function. galphat generates scale-free cumulative image tables for the desired model family with precise error control. Interpolation of this table yields accurate pixellated images with any centre, scale and inclination angle. galphat then rotates the image by position angle using a Fourier shift theorem, yielding high-speed, accurate likelihood computation.
We benchmark this approach using an ensemble of simulated Sérsic model galaxies over a wide range of observational conditions: the signal-to-noise ratio S/N, the ratio of galaxy size to the point spread function (PSF) and the image size, and errors in the assumed PSF; and a range of structural parameters: the half-light radius re and the Sérsic index n. We characterize the strength of parameter covariance in the Sérsic model, which increases with S/N and n, and the results strongly motivate the need for the full posterior probability distribution in galaxy morphology analyses and later inferences.
The test results for simulated galaxies successfully demonstrate that, with a careful choice of Markov chain Monte Carlo algorithms and fast model image generation, galphat is a powerful analysis tool for reliably inferring morphological parameters from a large ensemble of galaxies over a wide range of different observational conditions.
1 INTRODUCTION
The formation and evolution of galaxies is an outstanding problem in astronomy, and galaxy morphology remains a key observational attribute in the quest to increase our understanding of galaxy evolution. The increasing sensitivity and resolution of planned surveys will make possible tests of evolution theories from the epoch of formation to the present. However, selection effects and features peculiar to one’s choice of models will affect any interpretation. Therefore, to exploit the promise of survey data, we need to verify that our conclusions are reliable. The tools described in this paper are a step in this direction.
Early ΛCDM hierarchical galaxy formation theory and simulations placed galaxies in the Hubble sequence by following a combination of merger histories and gas accretion (White & Frenk 1991; Steinmetz & Navarro 2002). Later, ‘zoom-in’ resimulations of individual galaxies have produced more realistic galaxy morphologies (Abadi et al. 2003a,b; Sommer-Larsen, Götz & Portinari 2003; Governato et al. 2004; Robertson et al. 2004; Zavala, Okamoto & Frenk 2008). Although challenging, combinations of semi-analytic models and direct simulations (Croft et al. 2009; Parry, Eke & Frenk 2009; Benson & Devereux 2010; Scannapieco et al. 2010) have quantified the distribution of galaxy morphology with redshift.
These recent theoretical studies have been motivated by large-scale spectroscopic and image surveys. In the local Universe, millions of galaxies have been detected in the Sloan Digital Sky Survey (SDSS, York et al. 2000) and the Two Micron All Sky Survey (2MASS, Skrutskie et al. 1997; Skrutskie 2006). Recent analyses have used a range of models from single-component Sérsic profiles (Sérsic 1963) to more sophisticated bulge and disc–bar models to characterize the structural properties of local galaxy morphology (Blanton et al. 2003; Allen et al. 2006; Gadotti 2009). In the more distant Universe, COSMOS (Scoville et al. 2007) provides an ample collection of multiwavelength galaxy images and spectroscopy to study the evolution of galaxy morphological structure for z≲ 1 as a function of mass and environment (Capak et al. 2007; Cassata et al. 2007; Kovač et al. 2010). Gas accretion, disc instability, mergers and supernova and black hole feedback have been modelled to explain morphological evolution and it is possible to assess their relative importance by quantitatively comparing observed galaxy morphological structures to those predicted from theory. Furthermore, future large-scale multiband imaging surveys, such as the Large Synoptic Survey Telescope (LSST Science Collaboration 2009), combined with accurate photometric distances, will provide a uniform and consistent data set to study the evolution of the galaxy population.
Algorithmic approaches to measuring galaxy morphology are recent inventions and are usually based on mixture models of parametric surface brightness distributions (Byun & Freeman 1995; Simard 1998; Wadadekar et al. 1999; Peng et al. 2002; Simard et al. 2002; MacArthur, Courteau & Holtzman 2003; de Souza, Gadotti & dos Anjos 2004; Pignatelli, Fasano & Cassata 2006; Méndez-Abreu et al. 2008). However, systematic biases owing to an ignorance of uncertainties in the sky background and covariances between model parameters complicate their interpretation. To circumvent these difficulties, we advocate embedding the galaxy morphology analysis into the broader context of inference and hypothesis testing. In this paper, we present such a Bayesian approach using the Markov chain Monte Carlo (MCMC) technique, facilitated by embedding it within the Bayesian Inference Engine (bie; Weinberg & Moss, in preparation). To motivate this approach, we first illustrate the inherent difficulties in galaxy image decomposition in Sections 1.1 and 1.2.
1.1 Case studies
The following three examples explore the limitations of conventional model fitting and the improvements gained using Bayesian inference when inferring the photometric attributes of galaxies.
1.1.1 Posterior distributions versus best-fitting parameters
We pick two galaxies from our pool of Sérsic profile simulated galaxy images (see Section 1.1.2). One has a high signal-to-noise ratio, S/N = 100.01, and the other has a low signal-to-noise ratio, S/N = 10.46. The value of S/N is defined by the ratio of the galaxy signal to the noise within the half-light radius (see Section 4.1). Since we know the Sérsic model parameter values used to generate these galaxy images, we fix all the parameters to their true values, except for the Sérsic index n, and calculate the χ2 likelihood for different values of n.
Fig. 1(a) shows the likelihood as a function of n for each galaxy. The upper panel plots the high S/N galaxy and the lower panel the low S/N one. For the high S/N galaxy the likelihood has a very strong mode around n = 5 with a change by a factor of 4 in log, as n varies from 4.5 to 5.5. However, for the low S/N galaxy the likelihood is very broad, smoothly changing by 0.4 in log as n varies from 3 to 5, and the likelihood profile is not symmetric around the maximum. In addition, the precise location of the global maximum is not informative; a local analysis of the profile using standard inverse-Hessian analysis would not give an accurate estimate of the profile shape. Finally, in general, a typical multidimensional likelihood surface will have an even more complex landscape.
Panel (a): the likelihood as a function of Sérsic index n for two example galaxies with different S/N. Each likelihood value is normalized to the maximum. The distribution for the high S/N galaxy (upper) is sharply peaked and the distribution for the low S/N galaxy (lower) is broadly peaked. Panel (b): the posterior probability density of n for the two galaxies. The black dot with error bar is the best-fitting parameter from galfit (Peng et al. 2002). The shaded region is the 68.3 per cent confidence interval and the vertical dotted line is the true value of n. The conventional error estimate based on the second derivative of the likelihood is much too large for low S/N.
Fig. 1(b) shows the posterior probability of each galaxy’s n for the given data, marginalized over the other parameters as computed by galphat. The solid curve is the posterior probability of n and the shaded region corresponds to a 68.3 per cent confidence interval.1 The true value is indicated by the vertical dotted line, and the error bar shows the result using galfit (Peng et al. 2002). galfit is a widely used galaxy image decomposition program, based on a maximum likelihood (ML) approach, implemented as a χ2 minimization using the Levenberg–Marquardt algorithm (Press et al. 1998). The posterior mode of n is offset from the true value by 0.2 for the high S/N galaxy and by 1.0 for the low S/N galaxy. Such a bias always occurs owing to random photon counting errors. Although the bias is large for low S/N, the 68.3 per cent confidence interval of n encloses the true value. The best-fitting value of n from galfit for the high S/N galaxy is close to the posterior mode and the associated error bar, corresponding to a 68.3 per cent confidence interval, encloses the true value. However, for the low S/N galaxy the best-fitting parameter from galfit, using its simple minimization algorithm, is more subject to small-scale variations of the likelihood profile owing to sampling and the best-fitting parameter may also depend on the initialization of the downhill solver. Furthermore, the inverse-Hessian estimate of the variance in one dimension is simply the inverse of the second derivative of the likelihood; geometrically, the faster the slope varies with n, the smaller the variance. This makes the error reported by galfit unrealistically large because there is no significant variation in the tangent slope at the best-fitting value of n = 5.16. A reliable error estimate is crucial to quantifying trends in the derived properties of galaxies. The Bayesian MCMC approach adopted here samples the full posterior distribution and yields reliable error estimates of each model parameter given the data. Hence, this provides a solid statistical base for an analysis of galaxy morphology.
1.1.2 Bias and prior assumptions
We now explore the distribution of the best-fitting n for 325 galaxies with their Sérsic model index n sampled from a normal distribution2
. Fig. 2(a) plots the residual value of the n using the posterior median of n from galphat (blue diamonds) with a uniform prior for 
and the best-fitting value from galfit (red circles). For the galfit estimates, we used the true values as the centroids of the distribution for the initial guess and randomly perturbed the magnitude, galaxy radius re, axial ratio, position angle and sky background by ±0.5, ±20, ±10, ±15 and ±1 per cent, respectively, about these true values, assuming a uniform distribution within these ranges. In galphat we assumed uniform priors for these parameters within these same ranges. The initial galfit guess for the Sérsic index is n = 2.5 for all the galaxies. As an image’s S/N decreases, the variance in the residual value of n becomes larger and n becomes preferentially overestimated owing to the asymmetric shape of the likelihood profile (see Fig. 1). For low S/N, the galfit values are sensitive to the initial guess in contrast to galphat, which samples the parameter space using MCMC and hence is insensitive to the initial guess. The variance in the best-fitting n values from galfit is slightly larger than the variance in the galphat posterior medians owing to its insufficient accuracy in finding the correct global likelihood minimum for low S/N data.
The distribution of residual Sérsic index n for 325 galaxies with 
. Panel (a): galphat posterior median values (blue diamonds) and best-fitting galfit values (red circles). galphat uses uniform priors with the same parameter range used in galfit for a fair comparison. Panel (b): same as (a) but with a non-uniform prior for n in galphat. The prior distribution has a σ roughly 13 times larger than the input distribution. By introducing an informative prior, the parameter recovery is significantly improved with decreasing S/N.
However, if we allow the parameters to have an informative prior distribution, the galphat-derived posterior distribution improves dramatically (Fig. 2b). We use a Weibull distribution3 for the prior probability of n with λ = 7.0 and k = 1.5, whose deviation (σ = 4.3) from the mean is still 13 times larger than the deviation of the true distribution but embodies our astronomical experience from the literature. The variance in the posterior median of the residual of n for low S/N images decreases significantly. With an appropriately chosen prior distribution, the Bayesian approach can dramatically increase the quality of the inference over the entire input catalogue. Many people have aesthetic objections to the Bayesian approach, because they view the selection of a prior as being subjective and, therefore, arbitrary. The statistics and astrostatistics literature abound with philosophical discussions of this issue. It is our point of view that prior information will be used by a researcher inevitably, so why not let the Bayes theorem tell us how to use it quantitatively? Conversely, wholesale adoption of uniform priors, e.g. in the maximum likelihood method, is both arbitrary and self-defeating, since we know that most parameters have both physical constraints and previously measured distributions.
These issues are less significant for data with high S/N. However, astronomical surveys have widely ranging S/N values, and the total number of galaxies for flux-limited samples will always be dominated by images with low or moderate S/N. Therefore, improving our estimates of structural parameters in the low S/N regime best uses the available data and optimizes scientific return.
1.1.3 Scatter plots versus joint posterior distributions
Here we explore the joint distribution of galaxy half-light radius re and Sérsic index n for the same galaxy image sample from the last subsection. Fig. 3(a) is a conventional scatter plot showing the joint distribution of best-fitting values of re and n from galfit. Fig. 3(b) shows the joint posterior distribution of re and n from galphat. The contour levels are the 30, 50, 68.3 (green line), 95 and 99 per cent confidence values (white to black). Since the galaxy sample is generated without any correlations between re and n, we should not see any covariance between re and n. Comparing these two panels, the scatter plot from galfit shows a (spurious) systematic trend of n with re while the joint posterior distribution from galphat does not. Some Sérsic-model parameters, e.g. re and n, are reported to exhibit a true covariance (Trujillo, Graham & Caon 2001) and this parameter covariance is exacerbated by a degeneracy with the sky background, often subject to systematic bias owing to observational conditions, galaxy profile shape and image size, so such systematic trends in the distribution of best-fitting parameters as exhibited by galfit will obscure or contaminate any intrinsic covariance.
The joint distribution of galaxy half-light radius re and Sérsic index n for 325 galaxies. The simulated galaxies were generated without any correlation between re and n. Panel (a) shows the conventionally used scatter plot using the best-fitting model parameter from galfit and panel (b) shows the joint posterior of re and n from galphat with 30, 50, 68.3 (green line) 95 and 99 per cent confidence level (white to black).
The qualitative difference between the maximum likelihood (ML)-inferred scatter plot and the Bayesian-inferred joint posterior underlines our assertion that understanding parameter covariance, the use of thoughtful prior distributions and a thorough error analysis are essential to reliably testing hypotheses based on the morphologies of a large number of galaxies. Even so, the choice of a parametric family induces a correlation between parameters and, therefore, a single ‘best-fitting’ value does not adequately characterize the knowledge acquired from the data. Examples and conclusions such as these motivate our using the entire posterior distribution in parameter space for all of the galaxies that we wish to study.
1.2 Bayesian MCMC to the rescue
The main disadvantage of the Bayesian approach is its computational expenses. Over the last 15 yr, MCMC techniques have continued to improve and we believe that these techniques are now suitable as mainstream tools. Here, we introduce a computationally tractable Bayesian MCMC approach that overcomes the difficulties outlined in the previous sections. A significant improvement over previous attempts to understand galaxy evolution using morphology can be achieved using a Bayesian approach for the following three reasons.
The ML method assigns the best-fitting parameter value to the model that has the highest probability of generating the observed data. With sufficient data, this yields the correct result (assuming that some model in the family generates the data). However, Nature gives us one realization of a particular galaxy and this leads to competing models that can generate the same data through different processes. Rather than the ML question, we want to know the best model among the possible candidates, or what is the probability of the model for the observed data? Written in the language of conditional probability, this question gives us Bayes theorem and requires an estimate of the probability of the model before acquiring the data (the prior probability). The Bayesian formulation of the inference problem provides a natural foundation for astronomy, where every single event is unique and observers cannot test theories by changing the initial conditions of the Universe. If data have high S/N and strongly support a particular model, the inference should not be subject to the bias introduced by the prior distribution. Conversely, if data have low S/N, the inference may be influenced by the prior assumptions, as intuitively expected.
We have seen that inferences based on best-fitting analyses can be contaminated by intrinsic covariance in the chosen model family. In addition, the topology of the likelihood function in a high-dimensional parameter space with a large number of free parameters is very complex, in general. Therefore, one needs to both find the true global extremum and assess the significance of this extremum with respect to nearby and possibly unanticipated extrema in the probability space. Bayesian MCMC provides the full posterior and the confidence levels of inferences for each model parameter. Hence one can investigate correlations or perform hypothesis tests with quantifiable confidence.
The adoption of specific functional families, e.g. Sérsic profiles, may not provide an adequate explanation for the observed data or have sufficient power to classify differences between data sets. In the Bayesian paradigm, one may consistently assess the explanatory power of different models even if they are not nested. For example, one may rigorously pose the question ‘how strongly is the assumption of a single-component Sérsic model supported by the galaxy image data compared to a two-component bulge and disc model?’ This provides a natural way of probabilistically classifying galaxies.
To summarize, our Bayesian approach uses galaxy morphology as an intermediate step in an overall inference problem for theories of galaxy evolution. In contrast, by focusing on the best-fitting parameters as the data-reduction goal, and comparing the implied correlations to those predicted by theories of galaxy formation, one runs the risk of interpreting false correlations and one decreases the information content of one’s data by using only the best-fitting parameter as a summary value. Motivated by the promise of dramatic improvement, we have developed a novel image decomposition software package, galphat (GALaxy PHotometric ATtributes), based on a Bayesian Markov chain Monte Carlo approach. The general application to astronomical image analysis is not new: Bayesian MCMC has been used for X-ray surface brightness estimation of galaxy clusters (Andreon et al. 2008), for object detection (Hobson & McLachlan 2003; Savage & Oliver 2007; Carvalho, Rocha & Hobson 2009; Guglielmetti, Fischer & Dose 2009), for dynamical modelling of a galaxy (Puglielli, Widrow & Courteau 2010) and for gravitational weak lensing (Miller et al. 2007; Kitching et al. 2008) and strong lensing analyses (Vegetti & Koopmans 2009). galphat is designed for performing morphological analyses of galaxy images similar to several widely used galaxy image decomposition packages such as budda (de Souza et al. 2004), galfit (Peng et al. 2002), gasp2d (Méndez-Abreu et al. 2008) and gasphot (Pignatelli et al. 2006), but explicitly samples the Bayesian posterior distribution and, therefore, provides a comprehensive framework for the statistical inference of galaxy morphology by manipulating the sampled posterior. For example, marginals of the posterior distribution can be used to compute the confidence intervals on individual parameters that incorporate the full covariance. In addition, the galphat-produced posterior distribution can be used for a variety of inference problems such as model comparison, hypothesis testing and correlation analyses. We note that gim2d (Simard 1998; Simard et al. 2002) also uses a Metropolis algorithm to explore the posterior parameter space and could be used similarly for computing the parameter confidence intervals if the posterior distribution were unimodal.
One of our original motivations for developing galphat is the large-scale analysis of galaxy morphological structures in 2MASS (Skrutskie et al. 1997; Skrutskie 2006). The interpretation of 2MASS-galaxy properties has suffered from using unreliable best-fitting parameters obtained using conventional fitting algorithms. An ensemble of posterior distributions for a complete sample of local galaxies becomes a rich data base for hypothesis testing in two ways. First, we may characterize the distributions of galaxy morphological structures, e.g. the luminosity, the size and the shape, as a function of environment, with rigorous statistical confidence levels. Secondly and more generally, we may compare galaxy formation theories based on the morphological evidence.
This paper introduces galphat and emphasizes the methods, features and performance issues. We demonstrate the feasibility of a large-scale statistical inference based on galaxy morphology. A more detailed exploration of the influence of the prior distribution, explicit examples of model comparisons between single Sérsic and two-component bulge-and-disc models, and inferences using an ensemble of posterior distributions will be reported in a follow-up paper (Yoon, Weinberg & Katz, in preparation). The paper is organized as follows. In Section 2, we describe the basic formalism of Bayesian inference for galaxy image data and introduce the bie, which is used to sample the posterior distributions of model parameters. We describe the structure of galphat in Section 3 and our ensemble of simulated galaxy images for calibrating galphat in Section 4. Comprehensive test results are presented in Section 5. We summarize our findings and conclude in Section 6.
2 BAYESIAN MARKOV CHAIN MONTE CARLO
2.1 Theoretical background
, we recover the ML method.The goal is to characterize the posterior distribution by sampling P(θ|D). For simple cases with a few model parameters, one can analytically calculate it or explore the probability space by evaluating the posterior probability over a grid in parameter space. However, for our Sérsic model with eight parameters this approach is not feasible. Fortunately, owing to rapid improvements in computational methodology, MCMC is feasible in high-dimensional parameter spaces. MCMC generates states by a first-order Markovian process and the distribution of states asymptotically converges to the target distribution P(θ|D) after a large number of iterations.
. This quantity is a measure of how well the evidence supports the theory. In other words, the more probable the evidence given the theory, the more the evidence supports the theory. Of course, one needs to know what the theory predicts to know how well the evidence supports it and this is the job of the MCMC simulation. Now, let P(M) be the prior probability of the theory. Then, analogous to the Bayes theorem from equation (1), we may write the probability of the theory given the data as 
since P(M1) =P(M2). Therefore, the Bayes factor describes the increase of the odds in favour of one theory over another in the light of the data. The Bayesian model comparison does not depend on the particular parameters used by each model. Instead, it considers the probability of the model considering all possible parameter values.However, it should be clear from equation (2) that the Bayes factor depends on the choice of the prior distribution. If the likelihood dominates, the effect of the prior is negligible, but if the prior contributes to the posterior significantly as well as the likelihood, an incorrect prior leads to a biased inference (e.g. Kass 1993). The choice of prior must be considered carefully. The effects of prior choice on galphat will be addressed in detail in our follow-up paper (Yoon, Weinberg & Katz, in preparation).
2.2 The Bayesian Inference Engine (bie)
The ability to realize the promise of this approach depends on an accurate computation of the posterior distribution. Although the Markov chain will approach its steady-state distribution almost certainly, the number of steps required, the mixing time, is not known. In addition, the exploration of parameter space suffers from the curse of dimensionality.4 We need assurance that the Markov chain is in a steady state beyond the local region in parameter space. For example, consider a posterior distribution with discrete, separated modes; many chains will not be able to move between these modes, resulting in an infinite mixing time and an incomplete posterior distribution.
Various MCMC algorithms have been proposed to improve the convergence of MCMC, and each of these has its own advantages and disadvantages. Beginning in 2000, a multidisciplinary investigator team from the Departments of Astronomy and Computer Science at the University of Massachusetts designed and implemented the Bayesian Inference Engine, an MCMC parallel software platform for performing statistical inference over very large data sets. The bie uses a scalable multiprocessor software architecture designed to operate on modest cost, with generally available hardware. MCMC algorithms and Bayesian computation in general are ideally suited to multiprocessor computation. The bie uses standard MPI and POSIX threads and, therefore, will run in a broad spectrum of parallel or scalar environments and can be easily ported to high-performance hardware for production analysis. Fundamentally, the bie is a library, but the package provides a command-line interpreter (CLI) with access to nearly all of the import object classes. This CLI was originally intended for interactive or script-based prototyping with subsequent stand-alone hard-coding. However, most users simply use the CLI with scripts. The bie’s object-oriented design allows a researcher to apply a wide variety of MCMC algorithms to the same target application by changing several lines in a program or a script. The bie currently includes: the standard Metropolis–Hastings algorithm (Metropolis et al. 1953; Hastings 1970), the simulated tempering algorithm (Neal 1996), the parallel tempering algorithm (Geyer 1991), the parallel hierarchical sampler (Rigat 2008, PHS hereafter), the differential evolution algorithm (Braak 2006) and an independent multiple chain algorithm. For convergence testing, the bie implements the algorithm proposed by Giakoumatos et al. (1999) for single-chain simulations and the Gelman–Rubin (Gelman & Rubin 1992) convergence diagnostic for multiple-chain simulations.
The object-oriented design makes the bie extensible; new MCMC algorithms, new convergence algorithms and new models or likelihood functions can be implemented and added to the bie at any time. For example, a typical user will typically develop new models for a specific scientific problem. At a later time, the user may use a newly available MCMC algorithm without changing this model code. Since MCMC computations are computationally expensive, the bie provides a full serialization and persistence system. This system saves the entire state of all the objects in the simulation. For example, the bie automatically saves checkpoint images. Therefore, galphat can restart from the very last MCMC step to sample the posterior further for obtaining more MCMC samples when needed, significantly saving computational resources. Moreover, the results of previously performed simulations can be restored on the fly and compared or reused in new ways. See Weinberg & Moss (in preparation) for additional details.
3 galphat: ALGORITHMS AND FEATURES
galphat is implemented as a user-contributed likelihood function to the bie. It reads two-dimensional galaxy image data from an FITS file, generates a model image and then computes the likelihood. Because the posterior simulation requires a large number of likelihood evaluations, optimization of the model image generation is essential to making the analysis of an entire image catalogue feasible. In this section, we describe the implementation details and features of galphat.
3.1 Overview of our galphat implementation
The bie controls the MCMC algorithm, the convergence testing and logging the sampled posterior distribution. As needed, the bie requests a likelihood evaluation from galphat. The flow chart for galphat is as follows.
galphat reads the input FITS files [the galaxy image and the point spread function (PSF)] and the tabulated model images (see below) for later interpolation.
As a part of the MCMC simulation, the bie calls the likelihood function with a parameter vector. Using these parameters, galphat interpolates and scales the image table using the scale radius and the minor/major axial ratio, and generates a model image in the principal-axis coordinates.
galphat convolves the model image with input PSF image in Fourier space using the fftw package5 and adds the sky background.
Finally, the model image is rotated by the position angle, using a Fourier shift algorithm, in pixel coordinates.
galphat returns the likelihood evaluation to the bie.
Steps (ii)–(v) are repeated as necessary.
We will describe the details for each important step below.
3.2 Model generation
Panel (a): Sérsic surface-brightness profiles for n = 0.5, 1, 2, 4 and 8 (equation 4). The profiles are normalized to have equal flux density at r=re. Panel (b): the fraction of light within r for Sérsic profiles with n = 0.5, 1, 2, 4 and 8. For n = 8, a few per cent of the flux has r > 100re!
The fraction of light within r for different values of n is shown in Fig. 4(b). As summarized in Graham & Driver (2005), for an exponential disc (i.e. n = 1) profile, 99.1 and 99.8 per cent of the flux are within the inner 4re and 5re, respectively. For a de Vaucouleur profile (i.e. n = 4), 84.7 and 88.4 per cent of the flux are within the inner 4re and 5re, respectively. The sky background and index n become strongly covariant for small images and large-n profiles, and this biases estimates of n.
3.2.1 Image tables
. Similarly, the pixel scale can be non-isotropically scaled to obtain an arbitrary axial ratio q=b/a: 
. In the end, we must interpolate over our model grid. For Sérsic models indexed by ni, we use linear interpolation to obtain an approximation of C for ni≤n≤ni+ 1: 
3.2.2 Rotation of the model galaxy
x and 
y are shear operators in the x and y directions, respectively. Consider the function f(x, y) sheared in the x direction by a: f(x, y) →f(x+ay, y). Using the Fourier shift theorem, 
Since the galaxy model image is smooth and the flux values go almost to zero at the edges, aliasing should not cause any significant problems but we pad the images with zeros for added safety. In practice, we increase the image size by 
in each dimension to provide a sufficient margin for image trimming after rotation. We convolve the interpolated, unrotated image with the PSF before rotating the image in Fourier space. We also reduce the dynamic range of the surface brightness by a logarithmic mapping for large values of n. Then we rotate this PSF convolved image using the three-shear algorithm described above and apply the inverse logarithmic mapping if necessary. galphat uses the fftw package version 3.1.2 (Frigo & Johnson 2005) for all its FFTs.
Since convolving with the PSF before rotation smooths out high-frequency features in the profile, this image generation method can produce very accurate images without any significant aliasing introduced by the FFT rotation. Furthermore, one could subdivide the pixels of the image to perform the rotation computation and aggregate the pixels afterwards to increase the accuracy of the rotation. The error decreases exponentially with the number of subdivisions. For the test case described here, a subdivision by a factor of 2 decreases the error by a factor of 10. Of course, the computation time increases as the square of the subdivision factor.
As an example, we choose a 200 × 200 pixel n = 4 Sérsic model with re = 10 pixels. We generate Fig. 5(a) by direct numerical integration with a θ = 30° rotation and we generate Fig. 5(b) using the image generation method in galphat. The two images appear the same. Differences between the two methods only become obvious when one looks at a relative residual image. Figs 6(a)–(c) show a comparison of the same galaxy (i.e. n = 4 and re = 10) generated using the methods of galphat to the direct integrated image for three different axial ratios, q = 0.1, 0.4 and 0.7, respectively. We zoom in on the central region: [− 2re, 2re]. The left-hand columns in Figs 6(a)–(c) show the face-on views of the relative residual image and the right-hand columns are the views of the relative residual surface, with the magnitude of the residual plotted in the z direction. The maximum relative difference decreases with increasing q and remains much less than 1 per cent except for the extreme case of a concentrated galaxy with n = 4 and q = 0.1 (Fig. 6a), which is a very unrealistically small axial ratio for an observed galaxy with an n this large and still it only has a 5 per cent maximum error at the centre. The errors in the outer region are negligible and the total flux is still conserved to better than one part in 106 in all three cases. Galaxies with smaller Sérsic indices have smaller errors even when the axial ratio is smaller; such small axial ratios are more realistic for observed galaxies when the Sérsic index is small. The errors would be further reduced if we generated the images using pixel subsampling.
A comparison between the integrated and the interpolated and rotated image of a normalized Sérsic profile with n = 4. The axial ratio is 0.4 and the position angle is 30° counterclockwise from the x-axis. Panel (a) is an image produced by direct numerical integration of the profile and panel (b) is an image interpolated from the table and rotated using the Fourier shift theorem as in galphat.
![Differences between an image produced through direct integration and that produced using the method in galphat for a galaxy with n = 4, re = 10 pixels, and a position angle of 30°. The central region, [− 2re, 2re], is shown to highlight the differences. Panels (a), (b) and (c) are for q = 0.1, q = 0.4 and q = 0.7, respectively. In each panel, the left-hand column shows the face-on view of the relative residual image and the right-hand column shows a surface plot with the magnitude of the relative residual indicated on the z-axis.](https://oup.silverchair-cdn.com/oup/backfile/Content_public/Journal/mnras/414/2/10.1111/j.1365-2966.2011.18501.x/2/m_mnras0414-1625-f6.jpeg?Expires=1716070336&Signature=PXNdfExlwXaoob-FdbLw9lp1-d6bFCE3HHShYvHxts7b51KQOHzzGAa8DtRDkCIWIYcMAmZOX16sONyXi57CLgT2ObQEcZhMUIq7JDZiGhiyPYqqR0v50bG3D77aaj~swzMfX4hiE8V7a9wiPu0yAcJOZJd7Isp3pKjHO~6oR2eoPOS1EVI~XTDP8J0XtsCN4othu5XYWP5h77aAcR1wAJBx5x8lg0SZqlmIDRWaj-Pyiwhuk0uhwoYtliaON4i2qY4L3eULINtwsRyT6nZlDmZZnC4uCFMsIN8Onair2RSnYbwfwc22l1Lyr16X699P3VrZVWHySUOc8IUQfpHNBQ__&Key-Pair-Id=APKAIE5G5CRDK6RD3PGA)
Differences between an image produced through direct integration and that produced using the method in galphat for a galaxy with n = 4, re = 10 pixels, and a position angle of 30°. The central region, [− 2re, 2re], is shown to highlight the differences. Panels (a), (b) and (c) are for q = 0.1, q = 0.4 and q = 0.7, respectively. In each panel, the left-hand column shows the face-on view of the relative residual image and the right-hand column shows a surface plot with the magnitude of the relative residual indicated on the z-axis.
3.2.3 Computation time for model generation
The wall clock time for a posterior simulation depends on the model, the image and the MCMC algorithm. A typical run using the PHS algorithm requires 
evaluations (see Section 5.3). Here we provide a CPU time estimate for the generation of a single Sérsic model (Table 1) for an example image with a size of 190 × 190 pixels. The generation of a single Sérsic model requires 0.092 s (on a single processor). This is more than an order of magnitude faster than would be required using direct numerical integration over the pixels, and we still obtain a high numerical accuracy. Most of the CPU time is spent interpolating the image tables to obtain Ijk; the FFT rotation is a minor contribution. Since this image is larger than the mean galaxy size in our 2MASS target sample, we conclude that galphat is both sufficiently fast and sufficiently accurate to be suitable for a large-scale analysis over a large catalogue (see Section 5.3 for a discussion of the total MCMC run time).
galphat model-generation time
| Image/model | CPU | CPU time (total) | CPU time (interpolation) | CPU time (FFT rotation) |
| 190 by 190 | Quad-core AMD Opteron | 0.092 s | 0.079 s | 0.013 s |
| Sérsic | 2613 MHz |
| Image/model | CPU | CPU time (total) | CPU time (interpolation) | CPU time (FFT rotation) |
| 190 by 190 | Quad-core AMD Opteron | 0.092 s | 0.079 s | 0.013 s |
| Sérsic | 2613 MHz |
galphat model-generation time
| Image/model | CPU | CPU time (total) | CPU time (interpolation) | CPU time (FFT rotation) |
| 190 by 190 | Quad-core AMD Opteron | 0.092 s | 0.079 s | 0.013 s |
| Sérsic | 2613 MHz |
| Image/model | CPU | CPU time (total) | CPU time (interpolation) | CPU time (FFT rotation) |
| 190 by 190 | Quad-core AMD Opteron | 0.092 s | 0.079 s | 0.013 s |
| Sérsic | 2613 MHz |
3.3 The likelihood function and the prior distribution
. Assuming that each pixel is independent, our likelihood function L(D|θ) is 
As previously described, the prior distribution of model parameters affects the inference. A preliminary study with non-uniform prior distributions confirms that the posterior maximum and confidence values can be significantly changed for low S/N data that dominate the galaxy population in a flux-limited survey. However, for simplicity, we have used uniform prior distributions for all the parameters for the tests in this paper. Again, the use of an informative prior without careful consideration will almost certainly lead to a biased result. We present a detailed study of prior distributions for inference and model selection for galaxy profiles in our follow-up paper (Yoon, Weinberg & Katz, in preparation).
3.4 Sampling the posterior probability
It is difficult to ensure that the Markov chain correctly samples from the entire posterior distribution in a high-dimensional parameter space. To combat these difficulties, the bie uses a variety of algorithms, each with features effective for different problems that impede the efficiency of sampling and convergence. One should characterize the MCMC simulation on representative data using different algorithms before starting any production runs. We find that the PHS algorithm performs best for our problem; it more efficiently samples parameter space and more quickly reaches a steady state compared to the other MCMC algorithms. Here, we will briefly introduce the PHS algorithm (Rigat 2008).
3.5 Chain convergence
Monte Carlo simulations of the posterior distribution may suffer from two classes of difficulties: (1) the Markov chain is mixing poorly in a particular mode in the posterior distribution, and this leads to a large number of dependent states; and (2) the posterior distribution may have two or more discrete modes with similar posterior probability and the Markov chain cannot move between them. The first difficultly is easily diagnosed by observing very low or very high acceptance rates and is addressed by tuning the Metropolis–Hastings transition probability. The second difficulty is addressed by a variety of hybrid MCMC algorithms, implemented in the bie, designed to move between modes. The parallel chain algorithm, and tempering in general, decreases the contrast of the hills and valleys in the posterior distribution, which allows occasional large excursions between modes.
To test for convergence of the cold chain in the PHS algorithm, i.e. the chain with T0, we use a subsampled convergence diagnostic (Giakoumatos et al. 1999) for this single chain. The chain is cleaved and used to compare in-chain and interchain variances, similar to the Gelman–Rubin (Gelman & Rubin 1992) test. The ratio of these two variances should approach unity for a converged steady-state chain. We have tested galphat over a wide variety of synthetic image data typical of observed galaxies, using the empirically determined transition probability, and have confirmed that our MCMC algorithm samples the posterior with a reasonable acceptance rate of ≥25 per cent for good mixing and a swapping rate of ≥25 per cent for efficient mode exploration.
4 SIMULATED GALAXY IMAGES
To characterize the performance of galphat, we generated an ensemble of 3000 isolated, synthetic Sérsic galaxy images representative of survey observations. We vary both the S/N from 5 to 100 and re to probe the extremes of barely resolved galaxies and that of galaxies that extend beyond the image frame. The PSF is a Gaussian with a 2.96 FWHM in pixels, which we convolve with the model images. We use a Poisson noise model and a gain factor of 8.0 [e− ADU−1]. Both of these choices are motivated by 2MASS images. We describe the details of our choices below.
4.1 Varying S/N
For each S/N ∈{5, 10, 20, 50, 100}, we generate 100 Sérsic model galaxies with a randomly chosen combination of uniformly distributed re∈[6, 20], n∈[0.7, 7.0], axial ratio q∈[0.1, 1.0] and position angle PA ∈[0°, 90°]. We fix the sky background to 300 [ADU]. Once we choose re and q, we determine the galaxy’s total magnitude for the given S/N value and magnitude zero-point using equation (25).
4.2 Varying re
The inference of re is biased if the galaxy is small compared to the PSF. To test this, we generate 100 Sérsic model galaxies for all combinations of S/N ∈{5, 10, 20, 50, 100} and re∈{0.5, 1.0, 2.0, 4.0, 8.0} (in units of 
FWHM of the PSF) using the same distributions of n, q and PA as in Section 4.1. We compute the Poisson counting errors and the sky background as in Section 4.1.
5 galphat PERFORMANCE TESTING
Parametric surface brightness models invariably result in parameter covariance, and this covariance can be exacerbated by instrumental and selection errors. The Bayesian MCMC approach explicitly incorporates parameter covariance, noise sources and other selection effects including data S/N, PSF convolution and the sizes of the galaxy and the image, to yield a reliable inference, as illustrated in Sections 1.1 and 1.2.
Using the simulated galaxies from Section 4, we investigate the effect of observational attributes such as the S/N, the galaxy’s re compared to the PSF FWHM, the image size compared to galaxy’s re, and errors in the assumed PSF FWHM. Although we generated 40 000 converged MCMC states for each image, we do not use the full 40 000 MCMC states to construct the posterior since we want a more conservative estimate of the burn-in period than that diagnosed by the convergence test. As previously described, we use the PHS algorithm (Rigat 2008) for all of our tests and determine convergence using the subsampled convergence test (Giakoumatos et al. 1999). We tune the width of the Metropolis–Hastings transition distribution to yield, roughly, a 25 per cent acceptance rate and 25 per cent chain swapping rate for each chain. The prior distribution of the galaxy centroid is normal with a mean at the image centre and a σ = 3 pixels. The sky background prior is also normal with a mean of the input sky background and a σ = 0.5 ADU. The prior distributions for the other model parameters have a uniform probability within a finite range.
5.1 Examples of single fits
Before we present our results for an ensemble of galaxies, we present the results of galphat fits to four galaxies. From the simulated galaxy sample, we picked four example galaxies, two that represent elliptical galaxies (n≈ 4) and two that represent exponential disc (n≈ 1) galaxies. One galaxy of each type has a low S/N of 5 and one has a high S/N of 100.
In Figs 7–10 we show the marginalized posteriors of galphat model parameters: the total magnitude (MAG), the half-light radius (re), the Sérsic index (n), the axial ratio (q), the position angle (PA) and the sky background (sky). Figs 7 and 8 show the results at low S/N for an exponential disc-like galaxy and an elliptical-like galaxy, respectively, and Figs 9 and 10 show the results at high S/N for an exponential disc-like galaxy and an elliptical-like galaxy, respectively. In each figure, the full marginal distribution for each model parameter residual is shown on the diagonal with the vertical dotted line indicating the zero residual, and the joint marginal distributions of parameter pairs are shown on the off-diagonals with the seven colour contours corresponding to the 10, 30, 50, 68.3 (green solid line), 80, 95 and 99 per cent confidence levels. The locations of the zero residuals are indicated by the ×.
Posterior marginals for a galaxy with S/N = 5 and n = 1.57 (i.e. an exponential-disc galaxy). The full marginal distribution for each model parameter residual is shown on the diagonal with the vertical dotted line indicating the zero-point. Joint marginal pairs of parameter residuals are shown on the off-diagonals. The seven colour contours represent the 10, 30, 50, 68.3, 80, 95 and 99 per cent confidence levels and the green solid line is the 68.3 per cent confidence level. The locations of zero-points are indicated by × symbols. Posteriors are not unimodal and are spread over a large range in the parameter space. Although very weak covariances exist among magnitude, half-light radius (re), Sérsic index (n) and sky background, parameter uncertainties are largely dominated by Poisson random noise.
Posterior marginals for a galaxy with S/N = 5 and n = 4.25 (i.e. an elliptical galaxy). See the caption to Fig. 7. As for the exponential disc galaxy, parameter covariance is not significant and the posteriors spread out owing to noise.
Posterior marginals for a galaxy with S/N = 100 and n = 0.95 (i.e. an exponential-disc galaxy). See the caption to Fig. 7. Unlike the low S/N case, strong parameter covariances exist among magnitude, re, n and the sky background, and the parameter posteriors are confined in a narrow region close to the true value in parameter space.
Posterior marginals for a galaxy with S/N = 100 and n = 3.80 (i.e. an elliptical galaxy). See caption to Fig. 7. While parameter posteriors are compact and parameter uncertainty is dominated by parameter covariance, as for the high S/N exponential disc galaxy, there is a stronger parameter covariance among magnitude, re, n, and the sky background than for the exponential disc. This leads to inflated errors for individual parameters determined from the marginalized distribution.
When the S/N is small, for both an exponential disc galaxy (see Fig. 7) and an elliptical galaxy (see Fig. 8), the posteriors are not unimodal and are spread out over a large range in parameter space. Although there is a weak covariance between the magnitude, re, n, and the sky background, uncertainties in the parameter inferences are largely dominated by the Poisson random noise present in the image. However, the situation becomes very different at high S/N as shown in Fig. 9 for an exponential disc galaxy and in Fig. 10 for elliptical galaxy. Note that these figures have the same scale as Figs 7 and 8. Now the posterior forms a strong mode close to the true value and the morphology of the posterior is determined by the parameter covariance present in the Sérsic model, which strongly depends on the Sérsic index. There is a stronger parameter covariance among magnitude, re, n, and the sky background for the high S/N elliptical galaxy than for the high S/N exponential disc galaxy, which can lead to larger errors when marginalizing the posterior distribution.
The Bayesian-based galphat procedure explicitly incorporates the parameter covariance present in the Sérsic model. Furthermore, it enables us to utilize the entire posterior distribution for a galaxy population to reliably test hypotheses based on that population. This is practically feasible using galphat, as we will demonstrate in following sections.
5.2 Model covariance and bias
Using the posterior distribution from the converged Markov chain, we illustrate the inherent covariance between parameters by showing joint distributions of selected parameter combinations for the Sérsic model. To emulate a catalogue analysis, we generate an ensemble of galaxies with astronomically representative model parameters for each bin in observational conditions. The dependence on observational conditions is explored in the following several sections.
5.2.1 Effects of galaxy S/N
We characterize the posterior distributions of 500 images for Sérsic models (100 in each S/N bin with a range of structural parameters; see Section 4.1). Recall that these Sérsic models have eight free parameters: the centroid coordinates (x, y), the magnitude (MAG), the half-light radius re, the Sérsic index n, the axial ratio q, the position angle PA and the sky background (sky). We use the last 25 000 states to characterize the posterior. We constructed an ensemble posterior distribution by pooling the sampled distribution for all the images in each S/N bin.
We show all the marginalized distributions for errors in the magnitude, re, n, q, PA and the sky background for each S/N bin in Figs 11–15. Hereafter, we use the superscript i to denote the galphat inferred value. The parameters re and q are plotted as fractional changes scaled by their input values and the values of the other parameters are the differences from their input values. We quote the sky background error as the fractional per cent error, i.e. the fractional change scaled by the input sky background multiplied by 100. Each diagonal subplot is the full marginalized ensemble posterior error distribution of the corresponding parameter with a vertical dotted line indicating the location of the input value. Each off-diagonal subplot is the joint distribution of ensemble error posteriors for the corresponding parameter pair. The seven contour levels are the 10, 30, 50, 68.3, 80, 95 and 99 per cent confidence levels, and the green solid line marks the 68.3 per cent confidence level, corresponding to a ‘one-sigma’ normal confidence. The black crosses indicate the locations of the input values.
Posterior error distributions for the ensemble of galaxies with S/N = 5. The full marginal error distribution for each model parameter is shown on the diagonal. Joint marginal pairs of parameters are shown on the off-diagonals. galphat inferred parameters rie and qi are scaled by their input values and the other parameters are differences from their input values. The fractional sky background errors are percentages. The seven colour contours represent 10, 30, 50, 68.3, 80, 95 and 99 per cent confidence levels and the green solid line is the 68.3 per cent confidence level. The locations of the input values are indicated by vertical dotted lines for the diagonal and × symbols for the off-diagonals. The values of magnitude, re, n and sky background are strongly correlated. Although the constraints are tighter with increasing galaxy S/N, the strength of the parameter covariance increases with increasing galaxy S/N (see Figs 12–15).
Here, we use galaxies with 6 < re < 20 pixels, larger than the PSF that has an HWHM of 1.48 pixels to reduce the effects of resolution (the posterior distributions of small galaxies are described in Section 5.2.3). For this sample, q and PA are not covariates with the other parameters. However, the values of the magnitude, re, n, and the sky background are obviously covariate, and the covariance becomes stronger with increasing S/N. The origin of this covariance is straightforward to understand. For a given surface brightness profile, a magnitude underestimate (i.e. a luminosity overestimate) results in an overestimate of re to better match the observed brightness distribution. Since the Sérsic model parameters re and n have a positive correlation (e.g. Trujillo et al. 2001), n is also overestimated. Similarly, the sky background is underestimated to help compensate for the underestimated magnitude. This argument also holds exactly in the opposite direction for an overestimated magnitude. The shape of the joint posterior distributions in Figs 11–15 shows that the strength of the parameter covariance depends on galaxy S/N (and Sérsic index n as will show in Section 5.2.2).
In concert with our intuition, the confidence regions for q and PA shrink with increasing galaxy S/N. For example, the asymmetric heavy-tailed residual posterior distribution in ΔPA clearly seen in Fig. 11 becomes symmetric as S/N increases (Figs 12–15). This tail has its origin in the ambiguity of PA for q≈ 1.
The covariance of the magnitude, re, n and the sky background also changes with galaxy S/N. For low S/N (S/N = 5, Fig. 11), pairs of these parameters exhibit a clear covariance; the sky background is strongly covariant only with the magnitude. Also notice that the marginalized distributions of the errors in the magnitude and n, and of fractional errors in re, are not normal as is conventionally assumed and that the 68.3 per cent confidence region is not elliptical. This non-normal behaviour results from the lack of information at low S/N to constrain one or more of the covariate parameters.
As the S/N increases, the covariance of the magnitude, re, n and the sky background increases while the confidence regions decrease (see Figs 12–15) and the asymmetry of the marginalized posterior distributions vanishes. The posterior distribution is dominated by the likelihood function, which is sharply peaked and approaching its asymptotic form. In addition, the strength of the correlation depends on n owing to the strong correlation of n with the sky background. This can be seen in the joint posteriors of n and the sky background in Figs 13–15 where the confidence level contours appear to be a mixture of different covariance ellipses from groups of galaxies with different n, as will be shown in the next section.
5.2.2 Effects of Sérsic index n
To investigate trends with the Sérsic index n, we selected four model parameters to study in depth, magnitude, re, n and sky background, for S/N = 5, 20, 100. We divide the samples in each S/N bin into two groups: n > 2.0 and n≤ 2.0. We show the marginalized posterior error distributions for each group in Figs 16–18. The contours and curves for the samples with n≤ 2.0 and n > 2.0 are shown by the blue and red colours, respectively. The contour levels for each group corresponds to the 68.3, 95.4 and 99.7 per cent confidence levels. The black crosses are the locations of the input values and the grey contours and grey curves show the total sample as in Figs 11, 13 and 15.
Posterior error distribution for magnitude, re, n and sky background with S/N = 5 from Fig. 11 separated by n > 2.0 (red) and n≤ 2.0 (blue). Confidence levels are 68.3, 95.4 and 99.7 per cent and the input value is marked by a × and a vertical dotted line for joint and marginal distribution, respectively. Grey contours and curves represent the total sample. The parameter covariance is strong when n is large.
For S/N = 5 (Fig. 16), the marginalized error posterior for n≤ 2.0 has a sharp truncation on the left-hand side, owing to the prior distribution boundaries of 0.5 and 11.99 on n. Otherwise, the posterior error distributions of these two groups are similar. As expected for low S/N, the errors are dominated by random statistical errors and parameter covariance is not significant. However, for higher S/N (see Figs 17 and 18), the differences between low and high values of n are significant. When the S/N ≥ 10, the confidence regions for n > 2.0 galaxies (red) are larger than for those with n≤ 2.0 (blue). This is a consequence of Sérsic parameter covariance. The covariance is exacerbated by the degeneracy between the sky background and the extended profile for larger values of n. One can see the larger covariance for n > 2.0 in the joint posterior error distributions shown in Figs 17 and 18.
Moreover, covariance with the sky background significantly affects the inference of the magnitude and re. For example, when S/N = 100, a ±0.04 per cent variation in the sky background estimate may lead to an uncertainty in the magnitude of up to ±0.2 for elliptical galaxies, i.e. those with n > 2.0. Conversely, ignoring the sky background uncertainty can induce a significant bias in the other parameters that would be difficult to assess only using the best-fitting parameters from a simple χ2 minimization. Our approach incorporates the parameter covariance and random statistical uncertainty in any subsequent inference. Finally, we note that the parameter distributions of both large and small n galaxies are weakly biased at worst. This demonstrates that the galphat-inferred posterior maximum reliably recovers the true input parameters. We expect that careful attention to prior distributions could reduce the bias for lower values of S/N as well.
5.2.3 Effects of galaxyre
The intrinsic shapes of small galaxies are obscured by the PSF. Although any pixellization method induces some bias for any sized galaxy, if re is comparable to or smaller than the PSF width, the axial ratio q will be unrecoverable since the central pixels will contain most of the flux. Moreover, numerical techniques without explicit error control when computing I(x, y) (see equation 9) may fail to produce an accurate theoretical prediction for the model flux in the central pixels. galphat naturally addresses both problems. First, the Bayesian inference produces an estimate conditioned on the true value of the PSF and, therefore, will produce a posterior distribution consistent with all the possible models that yield the observed profile after convolution with the PSF. Secondly, the galphat model images are generated by interpolating a numerically integrated table that is accurate over the entire area of the galaxy with a numerical error tolerance set by the user. In particular, both the inner and outer profiles are well resolved because the tabulated grids are independent of scale. In other words, it is not possible for galphat to produce a poor estimate of the central pixel values for any value of re.
Even when accurately determined, the true value of q affects the other structural parameters because a small galaxy with a small q takes up fewer image pixels. Similar to Section 5.2.2, we investigate these effects in Figs 19–22 by dividing the galaxy sample into two groups: q > 0.3 (red) and q≤ 0.3 (blue), and by examining the posterior error distributions of magnitude, n, re, and q. We also select two bins in size: re = 0.5 times the HWHM of the PSF and re = 8.0 times the HWHM of the PSF, for two different values of S/N: 20 and 100.
Posterior error distributions for magnitude, scaled re, scaled q and residual n with S/N = 20 and re = 0.5 PSF HWHM, separating samples with q > 0.3 (red) and q≤ 0.3 (blue). Input values are shown by a ×. Contours are 68.3, 95.4 and 99.7 per cent confidence levels. The grey-scale contours and curves represent the total sample. When a galaxy is close to edge-one (i.e. q < 0.3) and smaller than the PSF, the marginalized posterior of q is almost flat and the posterior maximum of magnitude is negatively biased (i.e. the luminosity is overestimated).
Posterior error distributions for selected parameters with re = 0.5 PSF HWHM and S/N =100. See caption to Fig. 19.
Posterior residual distributions for selected parameters with re = 8.0 PSF HWHM and S/N =20. See caption to Fig. 19.
Posterior error distributions for selected parameters with re = 8.0 PSF HWHM and S/N =100. See caption to Fig. 19.
Figs 19 and 20 show the posterior error distributions for galaxies with re = 0.5 times the PSF HWHM for an S/N of 20 and 100, respectively. The three contours on the background grey contours, which indicate the total sample, denote the 68.3, 95.4 and 99.7 per cent confidence levels for galaxies with q > 0.3 (red) and q≤ 0.3 (blue). Since the PSF dilutes any evidence of the intrinsic shape inside the PSF FWHM, the posterior error distributions of re and q have significant probability density at larger values than the input values. Moreover, the magnitude is covariant with re and q. As a consequence, the posterior distribution of magnitude has significant probability at values smaller than the input value. Similarly, the bias in the Sérsic index n is exacerbated by its covariance with re.
In Fig. 19, the posterior error distribution for re has a tail to large re for both q > 0.3 and q≤ 0.3. As a result, the maximum of the posterior distribution in n is also similarly biased to large n for both groups. The bias in the two groups for q, however, is dramatically different: for q≤ 0.3 (blue) the distribution is approximately uniform over a large range, while for q > 0.3 (red) the distribution has a mode centred near 1. This is an artefact of the PSF convolution; the PSF naturally makes the galaxy appear rounder. Therefore, the bias in q for rounder galaxies is modest but the bias for flat, intrinsically edge-on galaxies is large.
The residual magnitude distribution for the entire sample (grey) is slightly negatively biased owing to the excess posterior probability of re at larger parameter values than the true value. The bias differences in q for the two groups leads to a bias difference in magnitude: the bias for q≤ 0.3 (blue) is larger than the bias for q > 0.3 (red). As described above, the posterior distributions for q are biased upwards as q decreases. On the other hand, because q is poorly constrained, the luminosity can be adjusted in a variety of ways to achieve the same surface brightness for a given re by making both q and the luminosity either large or small (see equation 8). This results in both a large spread in magnitude as well as a magnitude underestimate owing to a luminosity that is overestimated to compensate for the upward bias in the posterior distribution of q.
Fig. 20 shows the posterior error distributions for S/N = 100. As anticipated, the posterior distributions are better confined in parameter space and the biases of the posterior maxima are significantly reduced, e.g. the error posterior for n has a mode at 0. The posterior maxima for q > 0.3 galaxies (red) show no strong biases. The posterior error distribution for q when q≤ 0.3 (blue) reveals an extended tail and, in contrast to Fig. 19 for low S/N, it has a mode below 1.
As re increases beyond the PSF size, we expect these biases to decrease. Figs 21 and 22, which show the results when re = 8.0 HWHM of the PSF for an S/N of 20 and 100, respectively, confirms this expectation. The bias of the posterior maximum is modest and the differences between the q > 0.3 and q≤ 0.3 groups disappear. Also, for fixed re, the confidence regions decrease with increasing S/N, as expected.
5.2.4 Effects of image size
We have seen that Sérsic model parameters, such as the magnitude, n and re are covariant with the sky background. Therefore, an inaccurate sky background determination may bias the inferred galaxy structural parameters. To circumvent this, researchers have measured the sky background independently from their model fits. If the sky background measurement is not more precise than the model fit itself, this procedure has two disadvantages: (1) because of the covariance, the subsequent parameter posterior distribution will be biased; (2) characterization of the covariance will no longer be part of the posterior distribution.
The Bayesian MCMC approach implemented by galphat enables us to take all the parameter covariances into account in our galaxy modelling. In particular, we may characterize the influence of parameters such as the sky background. In this section, we explore the influence of the blank-sky fraction in the image. Assume that one must analyse a particular source on a large image. How much of the frame should one keep for one’s parameter inference? Retaining a large fraction of blank sky area is better to accurately determine the sky background. However, a larger fraction of blank sky implies a larger image size and more computation time. Of course, a truly blank sky is rare in Nature. None the less, an understanding of the trade-off between an accurate sky background determination and a fast model image generation is necessary to design an efficient analysis. We test the dependence on image size by generating 10 simulated galaxies with S/N = 50, re = 10 pixels, n = 4, q = 1 and a 500 [ADU] Poisson sky background and model them using different size image regions specified by 4, 8, 12, 16 and 20 times re. We add 30 000 converged MCMC states from each galaxy to obtain the posterior error distribution.
Fig. 23 shows the marginalized error posterior distributions for the sky background of the Sérsic model as a function of image size. The cross indicates the posterior maximum and the 68.3, 95.4 and 99.7 per cent confidence levels are indicated by the grey, red and blue boxes, respectively. The minimum and maximum data values are indicated by the error bars. Although the galphat inference of the sky background has a symmetric posterior error distribution around zero, regardless of the blank sky fraction in the image, the posterior maximum is below the true sky value as the blank sky fraction decreases. Since any galaxy surface brightness model, including a Sérsic profile, is degenerate with the sky background in the outer regions, a small image containing a relatively large galaxy may introduce a bias in the posterior distribution of the sky background and hence also affect the inference of the other parameters. In other words, owing to the Sérsic model surface brightness and sky background degeneracy, either increasing the galaxy luminosity and decreasing the sky background or decreasing the galaxy luminosity and increasing the sky background can match the background level of the image. However, this flexibility in a Sérsic profile to produce the sky background allows a galaxy to include part of the sky background flux and thus the sky background posterior tends to be biased downwards as we show in Fig. 24. This sky background bias, which stems from not enough information about the blank sky in the image, becomes more significant with increasing n. The galphat inference of the sky background for this particular case of n = 4 is not biased unless the image size is smaller than 12re and the bias in the sky background posterior maximum is very small, only 0.1 per cent, even when the image only extends to 2re.
Posterior error distributions of the sky background for images with different sizes. Ten galaxies with the same Sérsic model parameters but different random noise are modelled with different size image regions, from 4 to 20 times the input re. The galaxies have S/N = 50, n = 4, re = 10, q = 1 and a Poisson distributed sky background with 500 ADU. Within each bin, an ensemble parameter error posterior from the 10 galaxies is shown with the posterior maximum (cross symbol), 68.3 per cent (grey box), 95.4 per cent (red box) and 99.7 per cent (blue box) confidence levels. Minimum and maximum data values are indicated by the error bars. Although the sky background posterior is nearly symmetric regardless of the blank sky fraction in the image, the posterior maximum becomes slightly biased downwards with increasing confidence intervals as the blank sky fraction decreases.
Error posteriors of magnitude, re, n and sky background when the image region is 4re. The green line highlights the 68.3 per cent confidence level. Note that a 0.1 per cent bias in the sky background posterior maximum leads to biases in the other parameters, e.g. −0.3 in magnitude, and the inference for re is very weak.
However, the inferred distribution of the other parameters, i.e. the magnitude, re, and n, are remarkably sensitive to this tiny sky background bias, as shown in Fig. 24 where we plot the error posteriors of the magnitude, re, n and the sky background for a galaxy with n = 4 and an image size of 4re. In addition to the strong parameter covariances already illustrated in Sections 5.2.1 and 5.2.2, notice that the posterior maximum of the magnitude is biased low by 0.3 owing to the very small (0.1 per cent) bias in the sky background posterior maximum and that the posterior distribution of re is very broad. If the image size increases, these biases and weak inferences become less significant as shown in Fig. 25, which shows the case when the image size is 20re.
Error posteriors of magnitude, re, n and sky background when the image size is 20re. See caption to Fig. 24. Note that the parameter inference becomes more reliable with increasing image size.
Even this simple test using a small number of galaxies reconfirms the importance of an accurate sky background. It is tempting to subtract an independently measured sky background and then analyse small postage stamp images. However, the results of our simulations indicate that a variance in the sky background strongly biases other covariant parameters. Therefore, we advocate modelling the galaxy image including the sky background and characterizing the full posterior distribution of the model parameters. Alternatively, one could integrate the posterior distribution for a fixed sky estimate over the posterior distribution of independently determined background sky estimates to recover the full posterior distribution. This may be worth investigating in the future. Although galphat can handle the interactions of random uncertainties and parameter covariance over a wide range of image sizes, the bias is reduced with an image size of at least 10re for a galaxy with n = 4. This minimum image size will increase if the galaxy has n > 4.
The main goal for these simulations is to characterize galphat performance. In production analysis, one may use the survey catalogue or a prior estimate of the galaxy morphology distribution given the survey depth to estimate the relevant sizes for galaxies. In turn, this estimate allows one to choose an image size distribution a priori that appropriately limits the bias in the sky background. In addition, the consistency of one’s prior estimates could be checked a posteriori and repeated if necessary.
5.2.5 Effects of assumed PSF errors
We show the results in Fig. 26 where we plot the posterior maximum of relative error in re and n with the error bar indicating the 50 per cent confidence interval for the four different ΔFWHMs. If the PSF width were overestimated, we would expect that the inferred galaxy size would be smaller than its true size and that the inferred profile would be more concentrated, i.e. a larger n than the intrinsic value, and vice versa. Fig. 26 confirms these expectations. Also, as expected, the bias of the posterior maximum is less if the galaxy size becomes larger than the PSF FWHM. The bias of the re posterior maximum is larger when the PSF is overestimated than when it is underestimated as previously observed by MacArthur et al. (2003). Moreover, notice that there is a systematic offset of the ensemble posterior maximum of n even if the galaxy’s re is large, indicating that the bias of the n posterior maximum owing to errors in the assumed PSF FWHM also depends on n itself.
The posterior maximum errors of re and n for errors in the assumed PSF widths with ΔFWHM = 15, 5, −5 and −15 per cent as labelled. The error bars show the p = 0.5 confidence intervals. The biases in the posterior maximum become smaller with increasing galaxy size.
We characterize the n dependent bias by investigating the posterior error distribution of n for the galaxies with re = 8 HWHM of the true PSF. Fig. 27 plots the relative errors in the median values of the galphat posterior samples (ni) of n as 
for each galaxy in each of the PSF error samples: ΔFWHM =−15, −5, 0, +5 and 15 per cent. The median values of ni are biased and this bias becomes large as n increases. When the assumed PSF FWHM is overestimated, n is overestimated with increasing n since the larger PSF artificially extends the profile and conversely when the assumed PSF FWHM is underestimated, n is underestimated with increasing n since the smaller PSF artificially contracts the profile. The biases and uncertainties in ni for PSF FWHM overestimation are larger than for PSF FWHM underestimation. For example, the inferred value ni for an elliptical galaxy with n = 4 will be smaller than the true value by 20 per cent if the assumed PSF FWHM is smaller than the correct PSF by 15 per cent, but it will be larger by 30 per cent with considerable scatter if the assumed PSF FWHM is larger than the correct PSF by 15 per cent. If the correct PSF is used then the inference of n is unbiased but the scatter still increases with n.
Systematic bias in the inferred galaxy Sérsic index n owing to assumed PSF width errors based on the posterior error distributions for 100 galaxies with S/N = 100 and re = 8 HWHM of the PSF using four different assumed PSF FWHM errors with δHWHM =−15, −5, +5 and +15 per cent. Each diamond with error bar is the median and 50 per cent confidence range for the relative difference between the inferred value and the input value. The bottom-left panel shows results for the correct PSF.
5.3 Run time
The galphat run time depends on the MCMC algorithm, the complexity of model, the size of model image, the desired number of converged samples and the number of model components. Although it is difficult to characterize the galphat run time for each dependency, we find that the dependence on the number of converged MCMC states and on the size of the image, i.e. number of pixels, are approximately linear. However, the dependence on the MCMC algorithm (i.e. the temperature ladder and the number of chains, which depends on the parameter dimensionality) and the model complexity (i.e. the number of dimensions, the number of model components, and the parameter covariance in the chosen model family) are non-linear. Therefore, it requires some experimentation to find the best strategy for each application. By cross-checking with different MCMC algorithms, e.g. simulated tempering (Neal 1996) or differential evolution (Braak 2006) with tempering, for a subset of galaxy samples, we were able to tune the parameters for the PHS algorithm used in our study.
Table 2 shows an example of galphat run time for two galaxy images: one is typical of a 2MASS galaxy image and the other is a relatively large image, which we choose intentionally to demonstrate galphat’s feasibility in a marginal case as could occur for a 2MASS galaxy image, since we plan to analyse 2MASS galaxies in the future. Table 2 lists the total wall clock times of two example simulations for 40 000 converged states in 226 × 226 and 100 × 100 pixel images for the Sérsic model. The number of chains used for Sérsic modelling is eight. The bie assigns one processor per chain. The run time is the total time for obtaining the 40 000 converged samples using the PHS MCMC algorithm. The necessary number of converged samples depends on the application. The tests were performed using the University of Massachusetts Astronomy department HPC cluster.
galphat wall clock time.
| Image/model | Samples | CPU | Processors | Wall clock time | MCMC algorithm |
| 226 × 226 | 40 000 | Quad-core AMD Opteron | 8 | 1.5 h | Hierarchical tempering parallel chain |
| Sérsic | 2613 MHz | ||||
| 100 × 100 | 40 000 | Quad-core AMD Opteron | 8 | 0.4 h | Hierarchical tempering parallel chain |
| Sérsic | 2613 MHz |
| Image/model | Samples | CPU | Processors | Wall clock time | MCMC algorithm |
| 226 × 226 | 40 000 | Quad-core AMD Opteron | 8 | 1.5 h | Hierarchical tempering parallel chain |
| Sérsic | 2613 MHz | ||||
| 100 × 100 | 40 000 | Quad-core AMD Opteron | 8 | 0.4 h | Hierarchical tempering parallel chain |
| Sérsic | 2613 MHz |
galphat wall clock time.
| Image/model | Samples | CPU | Processors | Wall clock time | MCMC algorithm |
| 226 × 226 | 40 000 | Quad-core AMD Opteron | 8 | 1.5 h | Hierarchical tempering parallel chain |
| Sérsic | 2613 MHz | ||||
| 100 × 100 | 40 000 | Quad-core AMD Opteron | 8 | 0.4 h | Hierarchical tempering parallel chain |
| Sérsic | 2613 MHz |
| Image/model | Samples | CPU | Processors | Wall clock time | MCMC algorithm |
| 226 × 226 | 40 000 | Quad-core AMD Opteron | 8 | 1.5 h | Hierarchical tempering parallel chain |
| Sérsic | 2613 MHz | ||||
| 100 × 100 | 40 000 | Quad-core AMD Opteron | 8 | 0.4 h | Hierarchical tempering parallel chain |
| Sérsic | 2613 MHz |
For a Sérsic model, 40 000 converged samples for one galaxy with a typical image size of 100 × 100 can be obtained within 25 min of wall clock time using eight processors. This of course means it would take 3.2 h on a single-quad core processor or 0.5 d on a single core processor. This may seem computationally impractical but multicore processors and multiprocessor machines are becoming the norm and of course processor speeds and machine sizes are increasing all the time. For example, even now, each of our cluster nodes has 16 processors. Employing 32 nodes of our cluster for 1 d would yield the posterior distributions of about 3000 galaxies. This means that 10 000 2MASS galaxies can be analysed using a single Sérsic model in less than 4 d using 32 of our 16 core nodes, demonstrating the current feasibility of this approach.
6 SUMMARY
We introduce galphat, a Bayesian galaxy morphology analysis package, designed to efficiently and reliably generate the probability distribution of galaxy surface brightness model parameters from an image. We emphasize that the morphological analysis is a stepping stone in a larger chain of inference on theories of galaxy formation and evolution. Therefore, we believe that it is productive to consider the determination of galaxy morphology in the context of hypothesis testing and model comparison, and this demands the full posterior probability distributions for each galaxy in the ensemble.
In this section, we recap the history of morphological analyses, summarize the technical improvements offered by galphat, list major findings from our performance tests and briefly describe our future work.
6.1 Recap
Our approach offers a number of significant advantages for estimating surface brightness profile parameters. First, the topology of the likelihood function is almost certain to be multimodal in the high-dimensional space. Downhill optimization techniques demand precise prior information that is not generally available. In addition, one should have a way of assessing the significance of this global extremum with respect to nearby local and possibly unanticipated modes. Using the various tempering algorithms available in the bie, our tests have demonstrated that we can achieve a steady-state distribution and the simulated posterior will include any possible multiple modes supported by the prior distribution. Given the posterior distribution, we may then precisely estimate the confidence levels. Secondly, the model will often have correlated parameters. As illustrated in Section 1.1, any hypothesis testing that uses the ensemble of best-fitting parameters will be affected by these correlations. The full posterior distributions from galphat identify these correlations and incorporate them in subsequent inferences.
We can use posterior simulations over ensembles of images to test the significance of cluster and field environments on galaxies as evidenced in their photometric parameters. A more elaborate model might include angular harmonics of the light distribution (e.g. galactic bars and spiral arms); we could use Bayes factors to assess the support in the data for various harmonics. We could do this for an entire set of galaxy images in multiple bands simultaneously. This is much more natural and likely to be more powerful than the standard practice of cataloguing parameters and attempting to look for correlations in scatter plots.
Furthermore, the adoption of specific functional families that resemble galaxy profiles, e.g. Sérsic, may not provide the best discrimination in attempting to interpret the correlation between derived parameters and hypotheses. Our extensive study of Sérsic profiles undertaken for this paper convinces us that the strong interparameter covariances weaken the inference. This suggests that families of orthogonal functions conditioned to match the outer galaxy profile might be a more productive choice (e.g. Kelly & McKay 2004; Spergel 2010). A complete set of orthogonal functions might be straightforwardly transformed to match a fiducial profile, opening up a wide range of possible applications for characterizing galaxy properties. With carefully chosen prior distributions, we can use Bayes factors to test the significance of multiple component models and models with different functional forms. This analysis can straightforwardly answer questions such as: is the standard model (Sérsic and exponential disc profiles) preferred to all competing models (models with cores, dark-matter profiles) and vice versa? Is there a particular alternative model that is supported more significantly by the data? We will explore some of these possibilities in our follow-up paper (Yoon, Weinberg & Katz, in preparation).
6.2 galphat features
galphat is built on bie (Weinberg & Moss, in preparation), an object-oriented optimized parallel platform for Bayesian computation. The bie implements a variety of algorithms and tools that can be chosen at run time to match the problem at hand. For example, in the tests presented here, we have found that the PHS algorithm (Rigat 2008) provides the fastest convergence while efficiently exploring the parameter space. Other algorithms may excel for different model families or for those with additional components. In addition, the bie enables collaboration: new algorithms developed for the bie become available for the community.
These methods require more frequent likelihood evaluations than competing approaches, so optimization is essential. galphat incorporates the following key features: (1) pre-computed images using a scale-free relocatable interpolation scheme with strong error control; subsequent image generation accurately represents the surface brightness integrated over every pixel and is very fast; and (2) a Fourier shift-theorem-based rotation algorithm that is both accurate and much faster than the often-used interpolation methods. Although we explore Sérsic models in this paper, galphat can be applied to a wide variety of parametric families, limited ultimately by the physical memory available to store the lookup tables. The likelihood calculation time for a Sérsic model is less than 0.1 s for an image size of 190 × 190 pixels.
6.3 galphat performance on Sérsic profiles: parameter covariance and bias
A summary of our major findings are as follows:
Galaxy S/N. Computation of the posterior distribution enables assessment of parameter covariance that includes the full error model. For the Sérsic model, the magnitude, re, n and the sky background are strongly correlated, and the strength of the correlation increases with S/N. The covariance of the sky background with the other model parameters, i.e. the magnitude, re, and n, shows that a reliable inference requires an accurate characterization of the background noise. We fully expect that analogous trends will obtain for most parametric models.
Sérsic indexn. The parameter correlation increases with increasing Sérsic index n. As a consequence, the confidence intervals for magnitude, re, and n for galaxies with n > 2.0 are roughly three times larger than those for n≤ 2.0. However, the marginalized error posteriors of the sky background for these two groups do not have significantly different widths. Again, this underscores the need for the entire posterior distribution for subsequent inferences based on morphological quantities.
Galaxyre. If re is smaller than the PSF HWHM, most parameters are poorly constrained and the posterior maximum is biased. For example, the posterior distribution of re has a significant probability at larger re than the true value for a non-informative prior. Similarly, the inference of axial ratio q will be positively biased for intrinsically small q. Both biases lead to an overestimate of the total flux. As a galaxy’s re becomes larger compared to the PSF HWHM, this bias decreases. Also, the bias decreases with increasing S/N even if re is smaller than the PSF HWHM.
Image size. We find that the inferred value for the sky background is nearly symmetric about the true value even for an n = 4 Sérsic profile with an image size of 8 re. The background bias disappears and the uncertainty drops dramatically for image sizes larger than approximately 10 re. The covariance of sky background with galaxy surface brightness parameters motivates including the sky background in the overall model to ensure that the correlations are represented in the posterior probability distribution.
PSF variation. Errors in the assumed PSF width introduce a bias in the posterior distributions of galaxy size re and Sérsic index n. If the PSF FWHM is overestimated, the galaxy’s re is underestimated and its n is overestimated, and vice versa. The bias in the inferred value of n increases with n. The bias increase is larger if the PSF width is overestimated rather than underestimated.
Run time.galphat’s run time depends on the Monte Carlo algorithm, the desired size of the converged sample, the image size and the model complexity (i.e. the number of model parameters and the parameter covariance) and the computational complexity. Based on extensive experiments, we found that the PHS MCMC algorithm efficiently sampled the posterior for our Sérsic model tests. The typical wall clock times for generating 40 000 converged posterior MCMC samples of galaxies with image sizes of 226 × 226 and 100 × 100 pixels, using a Sérsic model are about 1.5 and 0.4 h, respectively, using eight 2-GHz AMD quad-core Opteron processors.
Although the optimized algorithms used in galphat significantly improve the likelihood computation time, it is still much slower than other conventional galaxy image decomposition algorithms. However, the existence of posterior probability distributions for an ensemble of galaxies enables reliable inferences for models of galaxy formation and evolution, and we feel that this more than compensates for the increased computational overhead. Moreover, the overhead will continue to decrease with the increasing availability and performance of HPC-class facilities.
6.4 Future work
We are investigating a two-component bulge-and-disc model composed of a Sérsic bulge with varying n and an exponential disc component with n = 1 using mutual independent prior distributions. This naive model exhibits strong correlations between multiple dimensions in parameter space, and these occasionally lead to poorly mixing Markov chains. From this preliminary study, we are convinced that a thoughtful prior distribution that removes astronomically unnatural regions of parameter space is required for production work. Yoon, Weinberg & Katz (in preparation) investigate galphat two-component bulge-and-disc model fits and introduce an informative prior based on typical 2MASS survey data. Using both bulge-and-disc two-component and single Sérsic models, we study the sensitivity of prior-distribution choice on the inference of galaxy parameters; demonstrate the practicality of model selection, the preference of a two-component bulge-and-disc model versus a single Sérsic model, using Bayes factor analyses; and illustrate the application of posterior distributions from an ensemble of galaxy images to large-scale inference problems.
In addition, we are currently using galphat to study the morphology of a Ks-band magnitude-limited sample of 2000 2MASS galaxies in the SDSS footprint, and will accurately characterize the luminosity functions of each component, i.e. the bulge and the disc, separately, and investigate any intrinsic correlations between the model parameters and any correlations with external galaxy properties such as the star formation rate and the galaxy environment in a forthcoming paper. We hope that galphat will become a well-used tool to aid in our understanding of galaxy properties in the near future and plan to release galphat to the community soon, when the companion methods paper (Yoon, Weinberg & Katz, in preparation) is published.
Here the confidence interval is estimated from the cumulative distribution function F(θ) of the marginalized parameter posterior probability density for the parameter θ. The 100(1 −α) per cent confidence interval [θ1, θ2] has 
and 
, where 0 < α < 1.
In all that follows we will denote the normal or Gaussian distribution with mean μ and variance σ2 as 
and the uniform distribution between a and b as 
. Other distributions will be introduced as needed.
The Weibull distribution is
The curse of dimensionality is the exponential growth of hypervolume as a function of dimensionality (Bellman 1961).
We thank our referee, Luc Simard, for encouraging and helpful comments. This work was supported in part by the NSF IIS Program through award 0611948 and by the NASA AISR Program through award NNG06GF25G. We thank Craig West for his support of the UMASS Astronomy HPC Linux cluster. We also thank Daniel McIntosh and Yicheng Guo for providing sample images and doing useful discussions in the early stages of this work.
REFERENCES
