Photometric Biases in Modern Surveys

TL;DR

This paper demonstrates that maximum-likelihood photometric estimators systematically overestimate flux, especially for resolved sources and in multi-band fitting, with implications for survey data analysis and bias correction methods.

Contribution

It reveals the extent and behavior of ML flux biases in modern surveys and provides correction prescriptions for these biases.

Findings

ML estimators overestimate flux, especially for resolved sources.

Bias increases with lower signal-to-noise ratio and more model parameters.

Simulations and real survey data confirm the presence of these biases.

Abstract

Many surveys use maximum-likelihood (ML) methods to fit models when extracting photometry from images. We show these ML estimators systematically overestimate the flux as a function of the signal-to-noise ratio and the number of model parameters involved in the fit. This bias is substantially worse for resolved sources: while a 1% bias is expected for a 10$\sigma$ point source, a 10$\sigma$ resolved galaxy with a simplified Gaussian profile suffers a 2.5% bias. This bias also behaves differently depending how multiple bands are used in the fit: simultaneously fitting all bands leads the flux bias to become roughly evenly distributed between them, while fixing the position in "non-detection" bands (i.e. forced photometry) gives flux estimates in those bands that are biased low, compounding a bias in derived colors. We show that these effects are present in idealized simulations, outputs from the Hyper Suprime-Cam fake object pipeline (SynPipe), and observations from Sloan Digital Sky Survey Stripe 82. Prescriptions to correct for the ML bias in flux, and its uncertainty, are provided.

Equations

$$\hat{\mathbf{b}}\sim \mathcal{N}({\mathbf{b}}^{\ast },\mathbf{C})$$

$${\hat{b}}_i\sim \mathcal{N}(\mu =b^{\ast },{\sigma }^2={\sigma }^2)$$

$$\mathbf{f}(x,y)=f{\mathbf{p}}_{x,y}$$

$$\hat{\mathbf{f}}(x^{\ast },y^{\ast },{\mathbf{b}}^{\ast })={\mathbf{f}}^{\ast }(x^{\ast },y^{\ast })+\hat{\mathbf{b}}\sim \mathcal{N}(f^{\ast }{\mathbf{p}}_{x^{\ast },y^{\ast }}+{\mathbf{b}}^{\ast },\mathbf{C})$$

$$\ln \mathcal{L}(x,y,f,b)=-\frac{nm}{2}\ln (2\pi {\sigma }^2)-\frac{1}{2{\sigma }^2}\sum _i({\hat{f}}_i-f{p}_i(x,y)-b{)}^2$$

$${\partial }_f\ln \mathcal{L}(x,y,f_{\text{ML}},b)=\frac{1}{{\sigma }^2}\sum _i({\hat{f}}_i-b)p_i(x,y)-f_{\text{ML}}{p}_i^2(x,y)=0$$

$$\boxed{f_{\text{ML}}(x,y,b)=\frac{{\sum }_i({\hat{f}}_i-b)p_i(x,y)}{{\sum }_i{p}_i^2(x,y)}}$$

$${\tilde{\sigma }}_f^2(x,y)\equiv -{\left({\partial }_f^2\ln \mathcal{L}\right)}^{-1}=\frac{{\sigma }^2}{{\sum }_i{p}_i^2(x,y)}\equiv \boxed{A_{\mathrm{psf}}(x,y)\times {\sigma }^2}$$

$$A_{\mathrm{psf}}^{\mathrm{G}}\to 4\pi s^2$$

$${\partial }_x\ln \mathcal{L}(x_{\text{ML}},y,f,b)=\frac{f}{{\sigma }^2}\sum _i\left(({\hat{f}}_i-b)-f{p}_i(x_{\text{ML}},y)\right){\partial }_x{p}_i(x_{\text{ML}},y)=0$$

$${\tilde{\sigma }}_x^2(x_{\text{ML}},y,f,b)=\frac{{\sigma }^2}{f^2}{\left(\sum _i{\left({\partial }_x{p}_i(x_{\text{ML}},y)\right)}^2-\frac{1}{f}\left(({\hat{f}}_i-b)-f{p}_i(x_{\text{ML}},y)\right){\partial }_x^2{p}_i(x_{\text{ML}},y)\right)}^{-1}$$

$${\tilde{\sigma }}_x^2(x_{\text{ML}},y,f)\approx \frac{1}{f^2}\frac{{\sigma }^2}{{\sum }_i{\left({\partial }_x{p}_i(x_{\text{ML}},y)\right)}^2}\equiv \boxed{S_{\mathrm{psf}}(x_{\text{ML}},y)\times {\left(\frac{f^2}{{\sigma }^2}\right)}^{-1}}$$

$$S_{\mathrm{psf}}^{\mathrm{G}}\to 8\pi s^4=2s^2{A}_{\mathrm{psf}}$$

$${\sigma }_x^2\to \frac{2s^2{A}_{\mathrm{psf}}^{\mathrm{G}}{\sigma }^2}{f^2}=2s^2{\left(\frac{f^2}{{\tilde{\sigma }}_f^2}\right)}^{-1}$$

$${\partial }_b\ln \mathcal{L}(x,y,f,b_{\text{ML}})=\frac{1}{{\sigma }^2}\sum _i{\hat{f}}_i-f{p}_i(x,y)-b_{\text{ML}}=0$$

$$b_{\text{ML}}(x,y,f)=\frac{1}{nm}\sum _i{\hat{f}}_i-f{p}_i(x,y)$$

$$\boxed{{\tilde{\sigma }}_b^2=\frac{{\sigma }^2}{nm}\equiv \frac{{\sigma }^2}{A}}$$

$${\hat{f}}_i=f^{\ast }{p}_i(x^{\ast },y^{\ast })+b^{\ast }+\sigma Z_i\sim \mathcal{N}(f^{\ast }{p}_i(x^{\ast },y^{\ast })+b^{\ast },{\sigma }^2)$$

$$\ln \mathcal{L}(x^{\ast },y^{\ast },f,b^{\ast })=-\frac{A}{2}\ln (2\pi {\sigma }^2)-\frac{1}{2{\sigma }^2}\sum _i((f^{\ast }-f)p_i(x^{\ast },y^{\ast })+\sigma Z_i{)}^2$$

$$f_{\text{ML}}(x^{\ast },y^{\ast },b^{\ast })=f^{\ast }+\frac{\sigma }{{\sum }_i{p}_i^2(x^{\ast },y^{\ast })}\sum _i{Z}_i{p}_i(x^{\ast },y^{\ast })$$

$$\frac{\sigma }{{\sum }_i{p}_i^2(x^{\ast },y^{\ast })}\sum _i{Z}_i{p}_i(x^{\ast },y^{\ast })\sim \mathcal{N}\left(0,\frac{{\sigma }^2}{{\sum }_i{p}_i^2(x^{\ast },y^{\ast })}={\tilde{\sigma }}_f^2(x^{\ast },y^{\ast })\right)$$

$$\boxed{f_{\text{ML}}(x^{\ast },y^{\ast },b^{\ast })\sim \mathcal{N}(f^{\ast },{\sigma }_f^2(x^{\ast },y^{\ast }))}$$

$$\ln \mathcal{L}(x^{\ast },y^{\ast },f^{\ast },b^{\ast })=-\frac{A}{2}\ln (2\pi {\sigma }^2)-\frac{1}{2}\sum _i{Z}_i^2$$

$$\sum _{i=1}^{nm}{Z}_i^2\sim {\chi }_A^2$$

$$(\hat{\mathbf{f}}-{\mathbf{f}}_{{\mathbf{\theta }}_{\text{ML}}}{)}^{\text{T}}{\mathbf{C}}^{-1}(\hat{\mathbf{f}}-{\mathbf{f}}_{{\mathbf{\theta }}_{\text{ML}}})\sim {\chi }_{A-p}^2$$

$$({\mathbf{\theta }}^{\ast }-{\mathbf{\theta }}_{\text{ML}}{)}^{\text{T}}{\mathbf{C}}_{\mathbf{\theta }}^{-1}({\mathbf{\theta }}^{\ast }-{\mathbf{\theta }}_{\text{ML}})\sim {\chi }_p^2$$

$$(\hat{\mathbf{f}}-{\mathbf{f}}_{{\mathbf{\theta }}_{\text{ML}}}{)}^{\text{T}}{\mathbf{C}}^{-1}(\hat{\mathbf{f}}-{\mathbf{f}}_{{\mathbf{\theta }}_{\text{ML}}})+({\mathbf{\theta }}^{\ast }-{\mathbf{\theta }}_{\text{ML}}{)}^{\text{T}}{\mathbf{C}}_{\mathbf{\theta }}^{-1}({\mathbf{\theta }}^{\ast }-{\mathbf{\theta }}_{\text{ML}})\sim {\chi }_A^2$$

$$f_{\text{ML}}^{\ast }\approx f_{\text{ML}}\left(1-\frac{X_2^2}{2}\frac{{\tilde{\sigma }}_{f_{\text{ML}}}^2}{f_{\text{ML}}^2}\right)$$

$$\boxed{\frac{{\delta }_{f_{\text{ML}}}}{f_{\text{ML}}}\equiv 1-\frac{\mathbb{E}[f_{\text{ML}}^{\ast }]}{f_{\text{ML}}}\approx \frac{{\tilde{\sigma }}_{f_{\text{ML}}}^2}{f_{\text{ML}}^2}}$$

$$\frac{\mathbb{V}[f_{\text{ML}}^{\ast }]}{{\tilde{\sigma }}_{f_{\text{ML}}}^2}\approx \frac{{\sigma }_{f_{\text{ML}}}^2}{f_{\text{ML}}^2}$$

Full text

methods: statistical — methods: data analysis — catalogs — techniques: image processing

\turnoffedit

1 Introduction

For example, the weak-lensing studies, supernova cosmology, classifying potentially hazardous asteroids, and to separate out main sequence and giant stars to map the galaxy .

the gold standard for precision photometry was photoelectric measurement with a photomultiplier tube (e.g. Landolt, 1973, 1983). The photons entering an aperture of, e.g., an eight arcsecond radius were detected and counted. In this way, hundreds of stars could be measured with great precision, and these “standard stars” formed the basis of the commonly used photometric systems. measurements of these standards were impressively consistent, but even the “Landolt Faint Standards” (Landolt, 1992) are relatively bright compared to the saturation limit of modern wide-area surveys.

Modern surveys use arrays of pixelated sensors such as charge-coupled devices (CCDs) in the optical/UV bands and mercury-cadmium-telluride (HgCdTe) devices in the near infrared. Surveys then perform “aperture photometry” by adding up all (background-subtracted) counts inside of a circle in an image. This process works well for isolated stars well above the background noise, and forms the basis of, e.g., the Sloan Digital Sky Survey photometric calibration (Padmanabhan et al., 2008).

Aperture photometry, however, has several significant issues that hamper its usage.

Shrinking the size of the aperture can this problem but then results in some amount of light from most sources being excluded, necessitating an “aperture correction” that needs to be calibrated . Aperture photometry can also become quite sensitive to issues relating to background estimation: because all counts in the aperture are being added together, these systematic offsets contribute an increasing portion of the counts with increasing aperture size .

point sources the appropriate 2-D “filter” is the point-spread function. Put another way, PSF photometry involves a parametric generative model of the noiseless data as a function of some parameters $\boldsymbol{\theta}$. Combining this with an appropriate noise model then yields a likelihood function. Maximizing the likelihood then provides an estimate of the flux with the highest signal-to-noise ratio.

Even for point sources, this process is not always straightforward. While aperture photometry simply requires an appropriately-sized aperture, PSF photometry relies on having an accurate model of the underlying source. This requires knowing the PSF precisely across the image to avoid applying a mismatched filter. While for space-based telescopes the PSF can be quite stable, for ground-based surveys it varies as a function of time and position in the focal plane to such a degree that it can cause flux errors of $1-2\%$, dominating the error budget of bright stars (e.g. Padmanabhan et al., 2008).

Furthermore, while PSF photometry can perform better than aperture photometry for fainter objects by avoiding “overweighting” the background, it still struggles in crowded regions where the measured parameters for objects become covariant with those of close neighbors (or even bright neighbors that are not so close). In the extremely crowded limit, the challenge of estimating the PSF and the background level may become severe. These questions are explored elsewhere (Brewer et al., 2013; Portillo et al., 2017) and are not the focus of this paper.

This maximum likelihood approach has been implemented in many software packages that are widely used in astronomy, including DAOPHOT/ALLFRAME (Stetson, 1987, 1994), DoPHOT (Schechter et al., 1993), HSTphot/DOLPHOT (Dolphin, 2000), Photo (Lupton et al., 2001), psphot (Magnier et al., 2016), The Tractor (Lang et al., 2016) in optimization mode, the Hyper Suprime-Cam (HSC) pipeline (Bosch et al., 2018), and the LSST Stack (Jurić et al., 2017). While these software packages differ on implementation details (e.g. estimating the sky level, measuring the PSF, and deblending sources), they are similar in how they measure the positions and fluxes of sources. (PSF/model fitting photometry for point/extended sources), possibly alongside aperture photometry at the same positions.

photometry (in theory) maximizes the SNR (i.e. is most precise), that does not mean that it is unbiased (i.e. is most accurate) or gives proper errors (i.e. has appropriate coverage). Indeed, that maximum-likelihood estimators are generically biased to order ${\rm SNR}^{-2}$. While there are many biases/systematic errors that can appear in photometry (cf. Nyland et al. 2017), in this paper we focus on the biases that arise in photometry due to it being a . We show that this leads the flux density estimates to generally be biased high, and that it broadly arises from PSF photometry “over-fitting” the data in a way that breaks symmetries. While this bias is generally small (a $\sim 1\%$ effect at $10\sigma$ for a point source) and likely not of concern for an individual source, its impact may be magnified when using a large population of low signal-to-noise sources.

In this work, we present corrections for the bias and consider its implications for . We first introduce a simplified version of maximum-likelihood photometry with a single source and Gaussian noise in Section 2. We derive the bias in the maximum-likelihood flux for a point source and the corresponding in Section 3. In Section 4, we compute the bias for multi-band photometry in both the forced photometry and simultaneous fit case. In Section 5, we show that the bias is present in mock data pipelines and SDSS data.

All the data and code used to create the plots presented in the paper are available online at https://github.com/joshspeagle/phot_bias. We invite readers to (re-)create their own plots and investigate the nature of this bias for themselves.

2 Maximum-Likelihood Photometry

Consider an $n\times m$ footprint of a pixelated image containing only one point source at some true position $(x^{*},y^{*})$ with some true flux density $f^{*}$. Assume the pixel-convolved point spread function (PSF) is constant throughout the footprint and is exactly known. Let the value of the pixel-convolved PSF in pixel $i$ for a point source located at position $(x,y)$ be $p_{i}(x,y)$. For notational convenience, the values across the whole image will be modeled as the $nm\times 1$ column vector $\mathbf{p}_{x,y}\equiv\mathbf{p}(x,y)=\{\dots,p_{i}(x,y),\dots\}$.

Assume that an estimate of the sky background111We define the sky background to be any flux not coming from the sources that are being cataloged, eg. scattered light from the atmosphere and telescope, dark current, and . has been subtracted from the footprint and that (residual) sky background in the footprint is Normally distributed with $nm\times 1$ mean column vector $\mathbf{b}^{*}$ and $nm\times nm$ covariance matrix $\mathbf{C}$, where the noise $\mathbf{C}$ is known but the mean bias $\mathbf{b}^{*}$ is not. The observed background noise in our pixels $\hat{\mathbf{b}}\equiv\{\hat{b}_{1},\dots,\hat{b}_{nm}\}$ is then distributed as

[TABLE]

where $\mathcal{N}(\boldsymbol{\mu},\mathbf{C})$ is the multivariate Normal distribution with mean vector $\boldsymbol{\mu}$ and covariance matrix $\mathbf{C}$.

In the case where the background noise is independently and identically distributed (iid) with mean $b^{*}$ and variance $\sigma^{2}$ the noise in each pixel follows

[TABLE]

and the mean and covariance become $\mathbf{b}^{*}=\{b^{*},\dots,b^{*}\}$ and $\mathbf{C}=\textrm{diag}(\{\sigma^{2}\})$. Although we will often use vector/matrix notation for compactness, we will derive results assuming the iid case throughout the main text. Some corresponding results for the general case are included in Appendix A.

Excluding noise from the sky background, the value of the model image for a point source at location $(x,y)$ and with flux density $f$ is

[TABLE]

The $nm\times 1$ observed flux densities $\hat{\mathbf{f}}$ within the footprint for our object with true flux density $f^{*}$ at position $(x^{*},y^{*})$ are then distributed as

[TABLE]

The log-likelihood for a model consisting of a single point source at location $(x,y)$ with flux $f$ and background $\mathbf{b}$ is then

[TABLE]

When extracting photometry, most often a maximum-likelihood (ML) approach is used. While ML estimators have been widely studied in the statistics literature, we derive some basic results for their application to photometry in this section for completeness. These results are already known in the literature and appear in, for example, King (1983).

2.1 Flux Density

Denote the maximum-likelihood (i.e. best-fit) flux for a given position $(x,y)$ and background $\mathbf{b}$ as $f_{\textrm{ML}}(x,y,\mathbf{b})$. By definition, the partial derivative of the log-likelihood at $f_{\textrm{ML}}(x,y,\mathbf{b})$ with respect to $f$ is zero such that

[TABLE]

where $\partial_{f}\equiv\partial/\partial f$. This yields

[TABLE]

which is equivalent to using the PSF centered at $(x,y)$ as a matched filter against the background-subtracted image.

We can construct a naive estimate of the error/uncertainty in the maximum-likelihood flux density estimate $\tilde{\sigma}_{f}(x,y)$ by calculating the (negative inverse of) the second derivative with respect to flux density

[TABLE]

where $A_{\rm psf}(x,y)\equiv 1/\sum_{i}p_{i}^{2}(x,y)$ is the “effective area” of the PSF. This uncertainty is proportional to the background noise $\sigma$ and effective PSF area $\sqrt{A_{\rm psf}(x,y)}$ (i.e. is larger when the PSF is broader). The estimate is naive because it ignores the possible covariances between flux density and other parameters.

As an example, a circular Gaussian PSF with a standard deviation of $s$ pixels has an effective area of

[TABLE]

in the oversampled limit where $(n,m)\gg s\gg 1$ (i.e. the footprint is large compared to the size of the PSF which is also large compared to the size of a single pixel).

2.2 Position

We can define the maximum-likelihood positions $(x_{\textrm{ML}},y_{\textrm{ML}})$ the same way by setting the partial derivative of the log-likelihood with respect to $(x,y)$ to 0. Using $x_{\textrm{ML}}$ as a representative case we get

[TABLE]

Similarly, the naive error/uncertainty in $(x_{\textrm{ML}},y_{\textrm{ML}})$ can be found by calculating the second derivative with respect to position:

[TABLE]

This expression has two components. The first involves the square of the first derivative of the pixel-convolved PSF $(\partial_{x}p_{i})^{2}$, while the second involves the second derivative of the pixel-convolved PSF $\partial^{2}_{x}p_{i}$ weighted by (fractional) model residuals. Assuming that the residuals are sufficiently small relative to $f$ and the PSF varies sufficiently slowly across the footprint, we can ignore this term and approximate the error as

[TABLE]

where $S_{\rm psf}(x,y)\equiv 1/\sum_{i}\left(\partial_{x}p_{i}(x_{\textrm{ML}},y)\right)^{2}$ is the “effective smoothness” of the PSF, analogous to the effective area $A_{\rm psf}(x,y)$. Note that the position error is directly proportional the effective PSF smoothness $\sqrt{S_{\rm psf}(x,y)}$ and inversely proportional to the signal-to-noise ratio (SNR) $f/\sigma$.

As above, for the case with a circular Gaussian PSF with a standard deviation of $s$ pixels, the effective smoothness is

[TABLE]

in the oversampled limit for $(n,m)\gg s\gg 1$. This gives a corresponding position error of

[TABLE]

where $\tilde{\sigma}^{2}_{f}$ is the naive flux density error estimate from §2.1.

2.3 Background

For a given flux density $f$ and position $(x,y)$, we can solve for the maximum-likelihood background $b_{\textrm{ML}}$ by again taking the first derivative

[TABLE]

and setting it to [math]. The ML solution is

[TABLE]

The associated naive errors then are

[TABLE]

where $A=nm$ is the area of the footprint.

This result shouldn’t be surprising, since it implies that the maximum-likelihood background is just the mean residual between the model $f\mathbf{p}_{x,y}$ and the data $\hat{\mathbf{f}}$ in our given $n\times m$ footprint. Since we have assumed a fixed value $b$ across the footprint, this is summed over all the pixels. Results for the general case where the background is actually a function of a $k$ nuisance parameters $\boldsymbol{\beta}$ across the footprint are outlined in Appendix C.

3 Bias in Maximum-Likelihood Estimates

While ML estimators are consistent estimators (i.e. $f_{\textrm{ML}}\rightarrow f^{*}$ as $f/\sigma_{f}\rightarrow\infty$) assuming no model mismatch, they are not guaranteed to be unbiased at any finite SNR. In this section, we derive estimates of the expected bias between the ML flux density estimate $f_{\textrm{ML}}$ and the true flux density $f^{*}$ along with its associated variance in the ideal case where $(x^{*},y^{*},b^{*})$ are known (§3.1) and the general case where it is not (§3.2). A schematic outline of our results is illustrated in Figure 1.

3.1 Ideal Case

It is helpful to rewrite the noisy observed flux densities $\hat{\mathbf{f}}$ in terms of random variables such that

[TABLE]

where each $Z_{i}\sim\mathcal{N}(0,1)$ is an iid random variable drawn from the standard Normal distribution with mean $\mu=0$ and variance $\sigma^{2}=1$. This represents a re-framing of the underlying data generating process: start with our true underlying model $f^{*}\mathbf{p}(x^{*},y^{*})+\mathbf{b}^{*}$ and add on a particular value comprised of $Z_{i}$ drawn from $\mathcal{N}(0,1)$ scaled by the standard deviation $\sigma$.

Given equation (18), the likelihood at the true position $(x^{*},y^{*})$ and background $b^{*}$ then reduces to

[TABLE]

Setting the partial derivative $\partial_{f}\ln\mathcal{L}=0$ then allows us to write the ML flux density estimate as

[TABLE]

We can rewrite this by noticing that the latter term is actually normally distributed such that

[TABLE]

where the equality comes from equation (8). This implies that the ML flux density estimate is distributed as

[TABLE]

since the naive error estimate $\tilde{\sigma}_{f}^{2}(x^{*},y^{*})$ is equivalent to the true error $\sigma_{f}^{2}(x^{*},y^{*})$ in the single parameter case . In other words, given the true values of $(x^{*},y^{*},\mathbf{b}^{*})$, the ML flux density estimate $f_{\textrm{ML}}$ is an unbiased estimate of the true flux density $f^{*}$ with a variance of $\sigma_{f}^{2}(x^{*},y^{*})$ corresponding to the true variance.

3.2 General Case

Following our random variable notation above, at the true values $(x^{*},y^{*},f^{*},b^{*})$ of the position, flux density, and background, respectively, we can rewrite the likelihood of the noisy data as

[TABLE]

where the $Z_{i}\sim\mathcal{N}(0,1)$ are again iid normal random variables. The sum of their squares represents the sum of the error-normalized residuals. We can rewrite this by recognizing that

[TABLE]

which follows a chi-square distribution with $A=nm$ degrees of freedom.

In general, we expect our best-fit parameters to “absorb” some of the scatter present in the data since we allow them to vary when we are trying to maximize the likelihood. We can make this more rigorous using , which implies that the sum of error-normalized residuals for a fit with $p$ free parameters $\boldsymbol{\theta}$ will follow

[TABLE]

with the sum of error-normalized residuals for the ML solution distributed as

[TABLE]

such that their combined sum leaves us with

[TABLE]

since $\chi^{2}_{i}+\chi^{2}_{j}\sim\chi^{2}_{i+j}$.

3.2.1 Decoupled Background

There is an asymmetry between the parameters connected to modeling the object and those connected to modeling the background. When modeling an object, $(f,x,y)$ can all modify the model image and provide information on the scale of the PSF. Even as we increase the size of our footprint $(n\times m)\rightarrow(\infty,\infty)$, there is a minimum variance $\sigma_{f,\min}^{2}$ achievable set by $A_{\rm psf}^{\max}$ and $\sigma$.

In contrast, the background estimate can continually improve as the image becomes larger. This holds true for any finite-parameter background model (see Appendix C) as the size of the footprint $A\rightarrow\infty$ becomes infinitely large. This implies that there is a fundamental difference between “object-related” parameters and “background-related” parameters.

In the case where the area of the footprint is substantially larger than the effective area of the PSF (i.e. $A\gg A_{\rm psf}$), we can consider the object parameters effectively decoupled from background parameters. In Appendix D, we show that the freedom in the object position parameters leads to a bias222Note that this effect is completely independent of (and unrelated to) selection effects such as Malmquist bias that can also cause sources to naturally be biased high close to detection limits/cutoffs. It is also independent of Eddington bias, which causes the number of bright objects to be overestimated when faint objects are more common than bright ones. in the generalized ML flux density $f_{\textrm{ML}}\equiv f_{\textrm{ML}}(x_{\textrm{ML}},y_{\textrm{ML}})$ relative to the ideal (unbiased) ML flux density $f_{\textrm{ML}}^{*}\equiv f_{\textrm{ML}}(x^{*},y^{*})$ in §3.1 that to leading order goes as

[TABLE]

where $X^{2}_{2}\sim\chi^{2}_{2}$ is a random variable drawn from the chi-square distribution with 2 degrees of freedom (which is determined by the 2 parameters $(x,y)$ used to fit for the position). This leads to a fractional bias of

[TABLE]

where $\mathbb{E}[f_{\textrm{ML}}^{*}]$ is the expectation value of $f_{\textrm{ML}}^{*}$. This shows that $f_{\textrm{ML}}$ is biased high relative to the true underlying flux density $f^{*}$. See Appendix D for a more detailed derivation involving higher-order terms and E for an alternate derivation using bias tensors. This bias is the same as the “gradient bias” identified by Ivison et al. (2007) and the “noise bias” identified by Refregier et al. (2012). A related derivation was done by Guy et al. (2010), which neglected the correlation between noise and the maximum likelihood position and obtained a different result. The authors of that paper confirm that when they include the noise-position correlation, they agree with our result (private communication).

This result gives a straightforward procedure to approximately “de-bias” $f_{\textrm{ML}}$ using equation (28). Doing so, however, increases the variance in the measurement to first-order by

[TABLE]

since the exact bias for any individual measurement is not known. This is an example of the bias-variance trade-off and leads to an increase in the total effective error following

[TABLE]

The magnitude of this bias for PSF photometry is generally small but not totally negligible: a nominally $10\sigma$ source (i.e. $f_{\textrm{ML}}=10\tilde{\sigma}_{f_{\textrm{ML}}}$) incurs a 1% bias and 0.5% error underestimate. We will return to this in §6.4 when we examine the behavior of $f_{\textrm{ML}}$ when modeling extended objects.

We test our analytic predictions by creating a set of simulated point-source images and running maximum-likelihood photometry on them. Our simulated images have a point-source with a circular Gaussian PSF of $\sigma=2\;\rm{pixels}$ in the center of a $101\times 101$ pixel image with iid Gaussian noise in each pixel. We simulate sources of nine fluxes ranging from $4.0\sigma$ to $9.4\sigma$, evenly spaced in $1/SNR$. For each flux, we create 100,000 different simulated images. We solve for the ML parameters using the Trust Region Reflective algorithm (as implemented in scipy), initialized with the true parameters to ensure that the optimizer finds the global maximum. Figure 2 shows that the mean flux bias and flux errors from these simulated images agree well with our predictions.

It is crucial to note this bias does not by itself arise from the fact that ML estimators tend to “overfit” the data when they use more parameters than are justified by the data. If that were the case, we should not have found that the ML estimate at the true position $f_{\textrm{ML}}$ was actually unbiased because the fitted parameters are exactly those used in the data-generating process. Instead, this bias arises because of the way in which this overfitting occurs. At the true position, $f_{\textrm{ML}}$ is allowed to adjust to the noise, but it does so in a symmetrical way: the noise fluctuations in the image are symmetric, and so $f_{\textrm{ML}}$ is just as likely to fluctuate upwards relative to $f^{*}$ as it is to fluctuate downwards. Once the position $(x,y)$ is allowed to vary, however, the model can move the source to absorb nearby positive noise. This breaks the symmetry from earlier: the source will tend to stay in the correct position with an overestimated $f_{\textrm{ML}}$ when the noise fluctuates upwards around the true position $(x^{*},y^{*})$, but will try to move position and absorb nearby positive noise when the noise fluctuates downwards around $(x^{*},y^{*})$. This position-dependent behavior of $f_{\textrm{ML}}$ is broadly illustrated in Figure 1 and shown in more detail in Figures 3 and 4. Averaging over this behavior as a function of position then gives the results derived above. We note that for sources that can be confidently detected (ie. ${\rm SNR}>5$), the ML positions only move from the true positions by fractions of a pixel (see Equation (14)). Thus the bias arises because the noise modifies the global maximum in the likelihood function, not because it creates a new local peak.

3.2.2 Coupled Background

The result we derived above holds if the image is sufficiently large that background estimation is effectively decoupled from modeling the actual object. In the case where this is not true (i.e. $A\not\gg A_{\rm psf}$), we instead need to consider how background estimation is covariant with our model parameters. The mixed 2nd-order derivatives of the log-likelihood give

[TABLE]

assuming that $p_{i}(x,y)$ is oversampled and roughly symmetric in $x$. As shown in , the contribution from the mixed partial with respect to $f$ and $b$ is expected to contribute to a fractional underestimate of the variance proportional to the ratio $A/A_{\rm psf}$ of the area of the footprint versus the effective area of the PSF. This gives a modified variance of

[TABLE]

since by construction $A\geq A_{\rm psf}$. Note that we have dropped $\tilde{\sigma}_{f}$ since $\sigma_{f}$ is now the “true” uncertainty of our ML estimator that takes into account the relative coupling between the ML background estimate $b_{\textrm{ML}}$ and the parameters used to model our object $(f_{\textrm{ML}},x_{\textrm{ML}},y_{\textrm{ML}})$.

Substituting in $\sigma_{f}$ for $\tilde{\sigma}_{f}$ into our expressions from §3.2.1 then modifies the effective SNR to give

[TABLE]

When the effective size of the footprint is large compared to the PSF, then $A/(A-A_{\rm psf})\rightarrow 1$ and the background is effectively decoupled from the modeling of the object, leading to our results in §3.2.1. When the effective size becomes more comparable to the PSF, then $A/(A-A_{\rm psf})\rightarrow\infty$ and the covariance between the background and flux density dominate the error budget. This has the effect of increasing the expected bias by decreasing the effective SNR.

While a coupled background worsens the bias that arises when the source position is unknown, no bias arises if the true source position is known and the background is uniform. In this case, the model image is linear in all the free parameters ($f$ and $b$) and so no bias is incurred. The covariance between the background level and flux density of the source will still inflate the uncertainty in the flux density, as shown in equation (34).

We augment the set of simulation images used to generate Figure 2 by considering five different image sizes approximately evenly distributed in $1/A$: 11, 13, 15, 23, and 101 pixels. The effective SNR of a given flux will decrease as image size decreases. In Figure 5, we show that our analytic predictions in terms of the effective SNR hold in these simulated images.

We want to note that while this effect is conceptually useful going forward, in practice the impact is extremely small. For instance, in SDSS the background is determined in patches of $256\times 256$ pixels.333https://www.sdss.org/dr14/algorithms/sky/ With a median seeing of $1.32\arcsec$ in $r$ band444https://www.sdss.org/dr14/imaging/other_info/ and pixel size of $0.4\arcsec$ (Gunn et al., 1998), $A_{\rm psf}/A=4.4\times 10^{-4}$. Similarly, in LSST the background will be determined in patches of $256\times 256$ or $512\times 512$ pixels.555https://confluence.lsstcorp.org/display/LSWUG/Measurement+in+the+LSST+Stack With a median seeing of $0.7\arcsec$ in $r$ band and pixel size of $0.2\arcsec$ (LSST Science Collaboration et al., 2009), $A_{\rm psf}/A=4.3\times 10^{-4}$ if the background is determined using $256\times 256$ pixel patches. A significantly larger effect is present in unWISE (Schlafly et al., 2019), which estimates the sky background in much smaller ($20\times 20$) pixel regions: with a PSF FWHM of $6\arcsec$ and a pixel scale of $2.75\arcsec$ (Wright et al., 2010), $A_{\rm psf}/A=2.8\%$.

4 Extension to Multi-Band Fitting

We now examine the case where an object is modeled in multiple bands. We will consider two cases. The first (§4.1) is where the object is “detected” in a single band, after which the position is fixed across all the bands. This is analogous to surveys such as The second (§4.2) is where the object is modeled simultaneously across all bands. This is somewhat equivalent to

4.1 Single-band Detection

Let’s assume that our object is detected in band $D$, after which the ML position $(x_{\textrm{ML}},y_{\textrm{ML}})=(x_{D}^{\textrm{ML}},y_{D}^{\textrm{ML}})\equiv(x_{D},y_{D})$ is fixed. As shown in §3, we expect that the ML flux density estimate in this band to be overestimated by an amount based on the PSF-normalized SNR. Following our previous assumption that the likelihood is multivariate normal around the ML parameters $\boldsymbol{\theta}_{\textrm{ML}}$, we would expect the values for our ML flux density estimates $f_{j,\textrm{ML}}$ in other band $j\neq D$ to be

[TABLE]

where $Z_{i}$ is again an iid normally distributed random variable, $\sigma_{j}$ is the noise in band $j$, and we have again assumed that background $b_{j}$ is known and fixed to the true value $b_{j}^{*}$. This result is analogous to equation (20), except that we have assumed a mismatch in the position between our model at $(x_{D},y_{D})$ and the source at $(x^{*},y^{*})$.

The expectation value of $f_{j,\textrm{ML}}$ assuming $(x_{D},y_{D})$ is fixed is

[TABLE]

Since we expect the mismatched PSF term to be smaller than the matched PSF term $\sum_{i}p_{i}(x^{*},y^{*})p_{i}(x_{D},y_{D})\leq\sum_{i}p_{i}^{2}(x_{D},y_{D})$, this implies that $\mathbb{E}[f_{j,\textrm{ML}}\,|\,x_{D},y_{D}]\leq f_{j}^{*}$ so that our ML flux densities are underestimated. Note that this bias tends to zero as $(x_{\textrm{ML}},y_{\textrm{ML}})\rightarrow(x^{*},y^{*})$, again confirming that our ML estimator is consistent in the limit of infinite SNR.

Subsequently taking the expectation value over position then gives the general expression

[TABLE]

where $P(x_{D},y_{D}\,|\,\boldsymbol{\theta}^{*},\mathbf{C}_{\boldsymbol{\theta}})$ is a 2-D multivariate normal distribution for $(x,y)$ (see Appendix B). While this integral does not have an analytic solution, since we expect the bias to increase as $(x_{D},y_{D})$ becomes progressively more offset from $(x^{*},y^{*})$, performing an average over possible ML positions further away from the true position should not be able to change our overall bias from an underestimate to an overestimate. Therefore, we arrive at the general conclusion that

[TABLE]

In other words, while our ML flux density estimates tend to be overestimated in the detection band, they will tend to be underestimated in all other bands. The severity of this (reverse) bias depends on the exact properties of the PSF in each band relative to the detection band (which establishes the ML position and associated covariances).

In the particular case where we have a circular Gaussian PSF with a standard deviation of $s$ pixels, we can evaluate these biases explicitly. At a fixed offset $r_{D}^{2}=(x_{D}-x^{*})^{2}+(y_{D}-y^{*})^{2}$, the expected value of $f_{j,\textrm{ML}}$ is

[TABLE]

Assuming that the covariance between $x_{D}$ and $y_{D}$ is small such that we can approximate the errors as $\sigma_{D,x}^{2}=\sigma_{D,y}^{2}$, the error-normalized offset will be distributed as

[TABLE]

This implies that

[TABLE]

where ${\rm Rayleigh}(\sigma_{D,x})$ is the Rayleigh distribution with scale parameter $\sigma_{D,x}$. . Marginalizing over $r_{D}$ then yields

[TABLE]

Assuming a Gaussian PSF in the detection band with standard deviation $s_{D}$ and $\sigma_{x,D}^{2}=8\pi s_{D}^{4}\sigma_{D}^{2}/f_{D}^{2}\approx 2s_{D}^{2}\sigma_{f,D}^{2}/f_{D}^{2}$, we get

[TABLE]

This corresponds to a fractional bias of

[TABLE]

In practice, we find that to properly model the bias at lower SNR requires incorporating a slightly higher-order Taylor expansion of our results (i.e. going from 2nd to 4th-order). As above, this expansion in general is non-trivial. However, it can be evaluated explicitly for circular Gaussian PSFs, as shown in Appendix G. Including this additional term then gives

[TABLE]

We test our predictions for forced photometry by creating a set of simulated point-source images in two bands, running maximum-likelihood photometry on one band (the detection band) and forced photometry on the other (the forced band). Our simulated images have a point-source with a circular Gaussian PSF of $\sigma=2\;\rm{pixels}$ in the center of a $101\times 101$ pixel image in each band with iid Gaussian noise in each pixel. We simulate sources with nine fluxes ranging from $4.0\sigma$ to $9.4\sigma$ in the detection band, evenly spaced in $1/SNR$. For each detection band flux, we consider four different forced band fluxes: $1\times$, $2\times$, $4\times$, and $8\times$ fainter than the flux in the detection band. For each flux combination, we create 100,000 different simulated images. Figure 6 shows that the flux in the detection band is overestimated and the flux in the forced band is underestimated, both by a fraction depending on the SNR in the detection band, as well as the bias on the measured color $C=-2.5\log(f_{D}/f_{j,\textrm{ML}})$.

4.2 Unforced Photometry in All Bands

When the object is modeled simultaneously across bands (i.e. with a common $(x,y)$ position), using Cochran’s theorem reveals that the sum of fluxes across bands will be biased high but does not show how the flux in individual bands is biased. To calculate the bias in each band, we use the bias tensor formulation introduced by Cox & Snell (1968), which shows that the leading-order term in this bias for parameter $s$ is

[TABLE]

where

[TABLE]

is the bias tensor and $\mathbb{E}_{\mathbf{D}}[\cdot|\boldsymbol{\theta}_{\textrm{ML}}]$ is the expectation value with respect to the data $\mathbf{D}$ for $\boldsymbol{\theta}_{\textrm{ML}}$ fixed.

The derivation of the bias in $f_{\textrm{ML}}$ with respect to $f_{\textrm{ML}}^{*}$ using bias tensors in the single-band case is outlined in Appendix E. There, we show that

[TABLE]

where

[TABLE]

When using multiple bands, the bias for band $i$ has the same terms as equation (49) for the single band case, but $\sigma_{x}^{2}=(\mathbfcal{F}^{-1})_{xx}$ is smaller because all bands help constrain the position. However, in the bias tensor (equation 50), the $\sigma_{x}^{2}$ is the uncertainty in position that would have been obtained using only band $i$. If all bands have Gaussian PSFs with widths $s_{j}$, then the flux bias for band $i$ is

[TABLE]

where the sum over $j$ is taken over all bands used in the fit. If all bands have the same PSF size $s$, then all bands have the same fractional flux bias of $(\sum_{j}{\rm SNR}_{j}^{2})^{-1}$.

We test our analytic predictions for simultaneous fitting by creating a set of two-band simulated point-source images with the same parameters as those used to make Figure 6 for forced photometry. Figure 7 shows that the fractional bias in both bands depends on the combined SNR, as predicted.

5 Application

While the results discussed in the previous sections are present in our simulations, we now turn our attention to real-world datasets to demonstrate that these effects are likely present in most datasets currently used in the astronomy community.

5.1 HSC SynPipe

We first investigate whether this effect is present in more realistic mock catalogs processed by real pipelines. In particular, we use simulated data from Hyper Suprime-Came Subaru Strategic Program (HSC-SSP; Aihara et al. 2018) Synthetic Object Pipeline (SynPipe; Huang et al., 2018). In brief, SynPipe injects fake objects into real images which are then processed by the HSC-SSP Pipeline (Bosch et al., 2018) to test the accuracy/precision of various aspects of the pipeline. These objects are drawn from a realistic color and magnitude distribution based in part on data from COSMOS, and so these mock tests represent fairly realistic realizations of the data seen by the HSC pipeline. See Huang et al. (2018) for additional details.

We analyze the PSF magnitudes for artificial star tests from two tracts (8764 and 9699) with good/poor seeing, respectively, processed using the same SynPipe configuration presented in Huang et al. (2018). The corresponding magnitude offsets and the predicted analytic relations are shown in Figures 8 and 9 for good seeing and poor seeing data, respectively. The magnitude offsets show good agreement with our model predictions, but include an additional systematic that is linear with magnitude. Investigating the source of this additional systematic is beyond the scope of this work666We have contacted many of the main authors of the pipeline so they are now aware of this systematic..

One crucial aspect of these results is that the SynPipe tests nominally represent forced photometric extractions, with detection done in the $i$ band. However, our results are almost entirely consistent with unforced photometry, where each band is derived separately. After some investigation, we find that this effect can be accounted for within the forced photometry algorithm used by the HSC pipeline, which effectively allows for limited “re-centering” in each band to improve the fit. Because the allowed range of positions is much larger than the relative positional uncertainties suggested by, e.g., $\sigma_{x}$ (§2.2) in most cases, this process effectively undoes the forcing effect described in §4.1. This phenomenon – where “forced” photometry from a particular pipeline is not quite what its namesake suggests – highlights the importance of transparency when pipelines provide users with results for conducting detailed analysis.

5.2 Stripe 82

To show that this bias appears in real data, we also look at SDSS catalogs of Stripe 82 (Annis et al., 2014). Stars that are low signal-to-noise in individual “runs” should have magnitudes that are biased high relative to their true values. While we do not have access to those true values, we approximate them using measurements taken from the combined images constructed from all the runs, which give much higher SNR measurements (with negligible bias) relative to the individual runs. We expect stars to be brighter, on average, in the individual run catalogs than in the stacked image catalog. Furthermore, each run and band will have a different bias, due to differences in seeing and sky brightness.

Figure 10 shows the magnitude difference between the individual run catalogs and stacked image catalog for each band and over a range of seeing conditions. The faintest stars are biased brighter in the individual runs, in rough agreement with our predictions. It is interesting to note that the apparent systematic trends seen in these data mirror those in the HSC SynPipe tests, and that the photometry is also described as “forced” photometry from the SDSS pipeline. As the HSC pipeline is in part derived from the SDSS pipeline, these similarities bolster our suspicions that these results are most likely caused by the same algorithmic choices.

6 Discussion

6.1 Aperture Photometry

We have shown that ML methods exhibit a generic bias when estimating the flux density of any particular isolated object from a given footprint. This bias causes flux density to be overestimated, and arises because fitting for an unknown object position does not treat noise fluctuations symmetrically. When the noise at the true position fluctuates low, the fit will be drawn away from the true position by nearby high noise fluctuations; conversely, when the noise at the true position fluctuates high, the fit will tend to remain near the true position. As is expected for any ML estimation bias, it is most pronounced at low SNR, scaling as ${\rm SNR}^{-2}$. This bias becomes worse as the models become more complex (as is the case for extended sources such as galaxies), and also biases colors when fitting across multiple bands with forced photometry.

Given these apparent drawbacks, some astronomers might wonder whether a return to aperture-based methods might present a compelling alternative. We want to offer a few arguments for why ML photometry should still be preferred and offer advice where aperture photometry might be more appropriate.

First, ML photometry performs better in the ideal case, where it is unbiased and follows the true error distribution, i.e. $f_{\textrm{ML}}^{*}\sim\mathcal{N}(f^{*},\sigma_{f}^{2})$. By contrast, an aperture will “miss” part of an object’s flux, leading to an underestimate in all cases. This correction is not expected to be the same for all sources unless the aperture is adaptively adjusted to match (a few times) the size of the PSF-convolved object, which is rarely the case. The “aperture corrections” involved to capture the total flux subsequently almost always serves as a dominant systematic hindering precise analyses.

Aperture photometry might also not eliminate the “centering bias” described in this work. Since an aperture also requires a position $(x,y)$ to be centered on, determining a central position for the aperture will likely be subject to the same types of biases as the ML case (§3.2).777Indeed, apertures are often centered using positions determined from either model-fitting approaches or various simplistic heuristics (e.g., “peak hunting”). These expected centering offsets will result in variable amounts of flux being excluded from the aperture, likely biasing aperture photometry to a similar extent as ML photometry. Unlike in the ML case where these biases can be studied using statistical methods, however, apertures by nature make such studies much more difficult.

The derived errors from aperture photometry are also generally larger than those from ML photometry. In the ML case, we showed in §2 that the error for a point source is $\sigma_{f}^{2}=A_{\rm psf}\sigma^{2}$, where $A_{\rm psf}$ is the effective area of the PSF. For a Gaussian PSF with a standard deviation of $s$ pixels, this gives $A^{G}_{\rm psf}=\pi(2s)^{2}$. Since aperture photometry just sums all pixels within a given aperture, the equivalent error for an aperture with radius of $r=2s$ pixels is just $\sigma_{f}^{2}=\pi r^{2}\sigma^{2}=\pi(2s)^{2}\sigma^{2}$. Any aperture larger than “2-sigma” then has errors that are strictly larger than those estimated from ML photometry, and even this $2s$ aperture excludes roughly 5% of the flux, requiring a significant aperture correction.

Aperture photometry is also inherently unstable as the aperture increases. While increasing the size of the aperture ensures a greater amount of the total flux is captured, it also increases the variance proportional to the size. While the SNR from ML photometry strictly improves as more data is added (see §3.2.1), the SNR from aperture photometry strictly degrades (ignoring aperture corrections).

Finally, aperture photometry is unable to integrate information across multiple bands/images. As discussed in §4.2, simultaneously fitting a single model across multiple images strictly improves the SNR and reduces the effective bias. Because apertures assume no model, they are unable to improve their SNR across multiple bands. While this comparison appears to be irrelevant in the examples shown in §5 which all exhibit tendencies equivalent to single-band unforced photometry, it will likely become more relevant in future survey pipelines.

We illustrate the magnitude of the above effects in Figure 11, which shows the magnitude difference between PSF and aperture (aper6) measurements from the same sources in the individual run catalogs as a function of the PSF magnitude measured from the stacked image catalog for each band and over a range of seeing conditions. There are substantial differences between the two measurements even at bright magnitudes that also show complex dependencies as a function of seeing. These most likely arise from the increasing sensitivity of the aperture to issues like background mis-estimation and incorrect aperture corrections. Once these effects are removed, we see behavior suggestive of a SNR*-2* bias. Compared to PSF photometry (see Figure 10), this bias becomes significant at much brighter magnitudes due to substantially larger aperture flux errors.

Ultimately, aperture photometry is appealing because it is so simple: it assumes no model and is straightforward to apply to almost any isolated object. While this leads to many of the drawbacks mentioned above and shown in Figure 11, it can also be desirable in cases where modeling complex sources can be difficult and/or the systematics involved limit the effective SNR of an object below that achievable with ML photometry. It thus serves a valuable purpose in cases where a model for the PSF and/or source cannot be cleanly determined and the source is relatively isolated. Our results suggest that it should only be used judiciously in most other cases.

6.2 Stacked Catalogs

One direct corollary of our results is that users must be extremely careful what exactly “stacking” means when constructing catalogs and estimating photometry from sources. This is important because stacking can occur at multiple levels, ranging from the images themselves to the catalogs produced from them. In the former case, where images are combined before they are processed through a pipeline, modeling the results is somewhat equivalent to the simultaneous fitting approach discussed in §4.2. Assuming that each image has roughly the same error $\sigma_{f}$, the effective error from $N$ images is expected to decrease to $\sigma_{f}/\sqrt{N}$ and the bias to decrease accordingly. Stacking on the image level thus reduces both the error and the bias.

In the case where stacking is done on the catalog level, however, each observation will be biased, with a mean of $f_{\textrm{ML}}=f^{*}(1+{\rm SNR}^{-2})$ and error of $\sigma_{f}$. Averaging the results over many identical catalogs then will reduce the error to $\sigma_{f}/\sqrt{N}$, but will not decrease the bias in any meaningful way. This implies that any measurement constructed from a stacked catalog may have systematic biases that far exceed the quoted statistical uncertainties.

Put another way, making catalogs from a series of images and then averaging the measurements across catalogs will not remove the flux density bias, because each catalog is individually biased. If inverse variance weighting is used and all images have the same PSF size, the fractional bias of the average will then be the reciprocal of the average ${\rm SNR}^{2}$. In contrast, the fractional bias from image stacking or simultaneous fitting is the reciprocal of the total ${\rm SNR}^{2}$, allowing multiple images of comparable SNR to drive down the bias. If catalogs are to be averaged, each catalog should first be individually debiased so that the average flux across catalogs is also unbiased. This procedure increases the variance in each catalog’s flux since the exact bias for each is not known (the bias-variance trade-off; see §3); however, this increase in variance is similar to that incurred by debiasing the flux measured from a stacked image or simultaneous fit.

For example, Budavári et al. (2017) propose detecting faint sources by stacking at the catalog level rather than the image level. This approach requires making catalogs with low SNR thresholds from each image, which is precisely the regime where the ML flux bias is most important. The catalogs made from each individual image should be debiased before being stacked.

We illustrate this effect in Figure 12 by constructing “stacked catalogs”888We “stack” our observations in flux density (linear) space to compute the weighted arithmetic mean. Note that stacking in magnitude (logarithmic) space (which computes the weighted geometric mean) introduces additional biases. using the same SDSS Stripe 82 data used to generate Figure 10. As expected, the bias remains unchanged regardless of the number of runs used to generate the stack even as the estimated errors (and thus the bias we would predict from the stacked catalogs) decreases substantially.

6.3 Bayesian Inference

The maximum likelihood bias we described can be ameliorated by using Bayesian inference. We show in Appendix F that marginalizing over the possible positions of the source cancels out the ML bias. Bayesian methods using Markov chain Monte Carlo sampling are starting to be implemented in astronomy. The widely-used Tractor code (Lang et al., 2016) returns ML measurements but also allows users to sample the posterior distribution of source parameters. Probabilistic cataloging (Brewer et al., 2013; Portillo et al., 2017) allows users to sample from the posterior in catalog space, yields an ensemble of catalogs that may vary in the number of sources identified. This catalog ensemble is useful in crowded fields to capture uncertainties in source identification and deblending. The computational cost required to sample from a posterior rather than reporting ML parameters can be daunting. As such, it is likely the case that ML methods will continue to be applied in both the near and intermediate terms in modern astronomical image processing pipelines.

6.4 Galaxies

7 Conclusion

In this paper, we study a photometric bias that arises from the maximum-likelihood (ML) estimator in model fitting photometry. It arises (in part) because the ML estimate can “soak up” a small amount of noise such that the fitted position is drawn away from the true position. This leads to an overestimate of the flux density, with a bias scaling with the inverse signal-to-noise ratio (${\rm SNR}^{-2}$) and the number of free parameters in the model. For example, it is 1% for a $10\sigma$ point-source and 2.5% for a $10\sigma$ 2-D Gaussian galaxy.

While this leads to an overestimate in the detection band, because the derived position is offset from the true position the flux density will be underestimated by the same ${\rm SNR}^{-2}$ in any other bands where the position is forced to the same value. This can double the effective bias in derived colors. By contrast, when all bands are modeled simultaneously, all bands are biased high, but less so than if they had been fit individually. In the case where all of the PSFs are the same size, this bias goes as $(\sum{\rm SNR}^{2})^{-1}$. If an object’s position is already known to great precision (for example, from a deeper or higher-resolution dataset), then forced photometry using this fixed position also does not suffer this bias. Methods that consider the distribution of possible positions, like Bayesian inference, do not exhibit this bias (see Appendix F).

We provide formulae that can be used in a variety of situations to calculate this bias and thus correct for it:

• •

Single-band point source photometry (Equation 29)

• •

Multi-band forced point source photometry (Equation 45)

• •

Multi-band joint point source photometry (Equation 51)

• •

(Equation LABEL:eqn:bias_par_gal)

We then show that this bias is likely common in many astronomical datasets using both mock HSC-like data and real SDSS data. The results further illustrate the importance of pipelines being transparent about the exact algorithmic implementation, since both tests are consistent with unforced photometry even though the data have been extracted using “forced” methods.

Although maximum-likelihood estimators may be biased, we still strongly encourage using them over simpler aperture-based methods in most cases. While apertures are appealing because of their simplicity, they require difficult-to-model aperture corrections to account for missing flux and likely exhibit similar biases due to offsets in aperture centers relative to the true positions of objects. Apertures also cannot effectively incorporate information across multiple bands, which can substantially reduce any relevant biases by improving the effective signal-to-noise ratio.

Though we have shown derivations in a simplified case, with iid Gaussian noise, Gaussian PSFs, and Gaussian galaxy profiles, this bias is generic to maximum-likelihood photometry estimation and would still arise if these assumptions were relaxed. This analysis could be repeated under different assumptions, like Poissonian noise, real PSFs, and Sersic galaxy profiles. While calculating the relevant corrections will likely be more involved in these more realistic cases, they are likely still tractable through the bias tensor formalism (See Appendix E) or through numerical simulations.

The authors would also like to thank Jim Bosch, Charlie Conroy, Daniel Eisenstein, Song Huang, Željko Ivezić, Ben Johnson, Douglas Scott, Eddie Schlafly, Sandro Tacchella, and Catherine Zucker for discussion and feedback that helped to improve the quality of this work. JSS is eternally grateful to Rebecca Bleich for her patience and support. JSS is supported by the National Science Foundation Graduate Research Fellowship Program. SKNP acknowledges support from the DIRAC Institute in the Department of Astronomy at the University of Washington. The DIRAC Institute is supported through generous gifts from the Charles and Lisa Simonyi Fund for Arts and Sciences, and the Washington Research Foundation. The Hyper Suprime-Cam (HSC) collaboration includes the astronomical communities of Japan and Taiwan, and Princeton University. The HSC instrumentation and software were developed by National Astronomical Observatory of Japan (NAOJ), Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU), University of Tokyo, High Energy Accelerator Research Organization (KEK) in Japan, Academia Sinica Institute for Astronomy and Astrophysics (ASIAA) in Taiwan, and Princeton University in the United States. Funding was contributed by the FIRST program from Japanese Cabinet Office; Ministry of Education, Culture, Sports, Science and Technology (MEXT); Japan Society for the Promotion of Science (JSPS); Japan Science and Technology Agency (JST); Toray Science Foundation; NAOJ; Kavli IPMU; KEK; ASIAA; and Princeton University. This paper makes use of data products derived from software developed for the Large Synoptic Survey Telescope (LSST). We thank the LSST Project for making their code available as free software at http://dm.lsstcorp.org. Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS web site is www.sdss.org. SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, the Chilean Participation Group, the French Participation Group, Harvard-Smithsonian Center for Astrophysics, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatório Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

Appendix A Results with Correlated Noise

Assume the noise in our footprint to be Normally distributed such that

[TABLE]

with mean vector $\mathbf{b}^{*}$ and covariance matrix $\mathbf{C}$. The corresponding log-likelihood for a point source within our footprint is

[TABLE]

where $\det$ is the determinant (i.e. the dimension-independent analog of area/volume).

A.1 Flux Density

At the ML flux $f_{\textrm{ML}}(x,y,\mathbf{b})$ the derivative with respect to $f$ is zero such that

[TABLE]

which yields

[TABLE]

The naive estimate (see §2.1) of the error/uncertainty $\tilde{\sigma}_{f}(x,y)$ is then

[TABLE]

A.2 Position

The maximum-likelihood positions $(x_{\textrm{ML}},y_{\textrm{ML}})$ can likewise be defined via

[TABLE]

with a naive error/uncertainty of

[TABLE]

A.3 Background

The maximum-likelihood background $\mathbf{b}_{\textrm{ML}}$ can likewise be defined using

[TABLE]

which gives

[TABLE]

The associated naive errors then are

[TABLE]

This result shouldn’t be entirely surprising. In §2.3, we noted that the maximum-likelihood background is just the mean residual between the model $f\mathbf{p}_{x,y}$ and the data $\hat{\mathbf{f}}$ in our given $n\times m$ footprint. In the iid case where we have assumed a fixed value $b$ across the footprint, we therefore take the average over all the pixels. In the case where every pixel has a separate possible value $b_{\textrm{ML},i}$ for the background, however, this leads to averaging done on a per-pixel basis for $\mathbf{b}_{\textrm{ML}}$. Since fitting the background separately in every pixel always over-fits the data by construction, we also derive results for the more realistic case where $\mathbf{b}_{\boldsymbol{\beta}}\equiv\mathbf{b}(\boldsymbol{\beta})$ is actually a function of a $k$ nuisance parameters $\boldsymbol{\beta}$ across the footprint in Appendix C.

A.4 Bias

As in §3, we can rewrite our likelihood in terms of random variable notation such that

[TABLE]

where each $Z_{i}\sim\mathcal{N}(0,1)$ is an iid random variable drawn from the standard Normal distribution and $\mathbf{C}^{1/2}$ is the (symmetric) square root of the covariance matrix. The likelihood at the true position $(x^{*},y^{*})$ and background $b^{*}$ then is

[TABLE]

This gives a ML flux density estimate of

[TABLE]

Since

[TABLE]

we then recover

[TABLE]

following §3.1.

Appendix B Results with Correlated Parameter Uncertainties

B.1 Error Underestimation in Two-Parameter Models

B.2 Background Covariance

B.3 Position Covariance

Appendix C Errors with General Background Models

In §A.3, we showed that the maximum-likelihood (ML) solution for a background model $\mathbf{b}$ across all $nm$ pixels is

[TABLE]

which had an error estimate of

[TABLE]

This result is singularly uninformative, because it implies that the “best” background model is exactly equal to the model residuals across the entire image.

In most cases, we often seek to parameterize the background as a smooth function $\mathbf{b}_{\boldsymbol{\beta}}\equiv\mathbf{b}(\boldsymbol{\beta})$ of $k$ nuisance parameters $\boldsymbol{\beta}$ across the footprint. This gives us

[TABLE]

which we can use to solve for $\boldsymbol{\beta}_{\textrm{ML}}$. Following Appendix B, the error estimates can be derived by inverting the Fisher Information Matrix (FIM) whose elements are

[TABLE]

As before, if we assume that our overall residuals $\hat{\mathbf{f}}-f\mathbf{p}_{x,y}-\mathbf{b}_{\boldsymbol{\beta}_{\textrm{ML}}}$ are small and that our background model varies sufficiently slowly with respect to $\boldsymbol{\beta}$, we can approximate this as

[TABLE]

where $\partial_{\boldsymbol{\beta}}\mathbf{b}_{\boldsymbol{\beta}_{\textrm{ML}}}$ is the $nm\times k$ Jacobian. This implies that we can use the Jacobian to linearly project from the $nm$-dimensional “data space” into the $k$-dimensional parameter space.

Appendix D Maximum-Likelihood Biases using Cochran’s Theorem

As discussed in §3.2, given the true values $(x^{*},y^{*},f^{*},\mathbf{b}^{*})$ of the position, flux density, and background, respectively, we can rewrite the likelihood of the noisy data as

[TABLE]

where $Z_{1},\dots,Z_{i}\sim\mathcal{N}(0,1)$ are again iid normal random variables and

[TABLE]

follows a chi-square distribution with $A=nm$ degrees of freedom. Assuming that the data are normally distributed and our ML parameters are also approximately Normally distributed, we can apply Cochran’s theorem to get

[TABLE]

Assuming the background is known (see §3.2.2), we note that we can relate the distribution of the sum of error normalized residuals around $f_{\textrm{ML}}\equiv f_{\textrm{ML}}(x_{\textrm{ML}},y_{\textrm{ML}})$ in the decoupled-background case to those around $f_{\textrm{ML}}^{*}\equiv f_{\textrm{ML}}(x^{*},y^{*})$ for a constant background model $b$ via

[TABLE]

where $X^{2}_{2}\sim\chi^{2}_{2}$ incorporates the noise “absorbed” by $(x_{\textrm{ML}},y_{\textrm{ML}})$ and we have exploited the fact that $X^{2}_{2}+X^{2}_{nm-2}\sim\chi^{2}_{nm}$. Exploiting the fact that

[TABLE]

then allows us to rewrite the above result as

[TABLE]

Although the true position is not known, in the interest of deriving the impact on the flux we can take the approximation that $(x^{*},y^{*})\approx(x_{\textrm{ML}},y_{\textrm{ML}})$, etc. for all terms that don’t explicitly involve the flux density $f$. This leaves us with an equation of the form

[TABLE]

where $c_{1}=\tilde{\sigma}_{f_{\textrm{ML}}}^{-2}$ and $c_{2}=\hat{\mathbf{f}}^{\textrm{T}}\mathbf{C}^{-1}\mathbf{p}_{x_{\textrm{ML}},y_{\textrm{ML}}}$ are roughly constant. This has a positive solution at

[TABLE]

Since $c_{1}$ and $c_{2}$ are known for a given $(x_{\textrm{ML}},y_{\textrm{ML}})$ and the distribution of $X^{2}_{2}\sim\chi^{2}_{2}$ is known exactly, this gives an expression for the distribution of the unbiased ML estimator $f_{\textrm{ML}}^{*}$. In general, assuming that the residuals are sufficiently small such that $c_{2}\approx 0$, this reduces to

[TABLE]

This can be immediately generalized to a model with $p$ model parameters $\boldsymbol{\theta}$ (excluding the background $b$) to get

[TABLE]

We can write this in a slightly more intuitive form by Taylor expanding around small $\tilde{\sigma}_{f_{\textrm{ML}}}/f_{\textrm{ML}}$ to get

[TABLE]

This gives a fractional bias of

[TABLE]

At lower SNR it is necessary to include the second-order term $\frac{X_{p-1}^{4}}{8}\frac{\sigma^{4}_{f_{\textrm{ML}}}}{f_{\textrm{ML}}^{4}}$ from the Taylor expansion to properly model behavior. Including this term gives a modified fractional bias of

[TABLE]

For a typical point source model with $p=3$ parameters $(x,y,f)$, this becomes

[TABLE]

Appendix E Maximum-Likelihood Biases using Bias Tensors

ML estimators have a bias $\delta$ which tends to zero as the signal-to-noise ratio (SNR) increases. Cox & Snell (1968) found that the leading-order bias term for any parameter $s$ can be found with

[TABLE]

where

[TABLE]

is the bias tensor and $\mathbb{E}_{\mathbf{D}}[\cdot|\boldsymbol{\theta}_{\textrm{ML}}]$ is the expectation value with respect to the data $\mathbf{D}$ for $\boldsymbol{\theta}_{\textrm{ML}}$ fixed.

With the background $b=b^{*}$ fixed, the non-zero terms in the flux density bias $\delta_{f}$ are

[TABLE]

since the off-diagonal elements of the FIM with respect to position $(x,y)$ are zero (§B.3) and we have substituted in for $(\mathbfcal{F}^{-1}(\boldsymbol{\theta}_{\textrm{ML}}))_{ff}=\tilde{\sigma}_{f_{\textrm{ML}}}^{2}$ and $(\mathbfcal{F}^{-1}(\boldsymbol{\theta}_{\textrm{ML}}))_{xx}=(\mathbfcal{F}^{-1}(\boldsymbol{\theta}_{\textrm{ML}}))_{yy}=\sigma_{x}^{2}$. Under similar assumptions as §3.2, is it straightforward to show that

[TABLE]

and thus:

[TABLE]

Substituting in for the definition of $\delta_{f}$ allows us to rewrite this as

[TABLE]

which reproduces equation (28).

With the background $b$ free, the flux bias has the same terms as Equation E3 except that $(\mathbfcal{F}^{-1}(\boldsymbol{\theta}_{\textrm{ML}}))_{ff}=\sigma_{f_{\textrm{ML}}}^{2}\geq\tilde{\sigma}_{f_{\textrm{ML}}}^{2}$ due to the covariance between the flux and background. Solving and rearranging as above then gives

[TABLE]

which reproduces equation (29).

Appendix F Unbiasedness of Marginalized Mean Flux

While the ML solution is the mode of the likelihood distribution, the mean flux density $f_{\rm mean}$, marginalizing over position $(x,y)$, need not be the ML flux density $f_{\textrm{ML}}$. Here, we show that the posterior mean flux is unbiased to order $SNR^{-2}$ in a specific case (flat priors and Gaussian PSF) for illustrative purposes.

As discussed in §2, the maximum likelihood flux $f_{\textrm{ML}}(x,y)$ varies with position and is maximized at $(x_{\textrm{ML}},y_{\textrm{ML}})$ assuming the background $b$ is known. While first derivative of $f_{\textrm{ML}}\equiv f_{\textrm{ML}}(x_{\textrm{ML}},y_{\textrm{ML}})$ with position is zero (by definition), the second derivative is

[TABLE]

under the same oversampled/smoothness assumptions discussed in §2.2. For a circular Gaussian PSF with standard deviation of $s$ pixels, the second derivative of the PSF is

[TABLE]

where $\Delta x_{i}$ is the difference in the x-coordinate of the center of pixel $i$ and the source. Taking

[TABLE]

and approximating the counts around the source as $\hat{f}_{i}\approx f_{\textrm{ML}}p_{i}$ then gives

[TABLE]

using the definition of $f_{\textrm{ML}}$ from equation (7).

Near the ML position, the ML flux density is approximately

[TABLE]

The likelihood can then be approximated (up to a constant factor $Z$) as

[TABLE]

where we have taken $\tilde{\sigma}_{x}^{2}=\tilde{\sigma}_{y}^{2}$ since our PSF is circular. The mean flux density $f_{\rm mean}$ after marginalizing over position is defined as

[TABLE]

Integrating over the flux density $f$ yields

[TABLE]

and then using equation (F5) for $f_{\textrm{ML}}(x,y)$ gives

[TABLE]

Substituting in $\sigma_{x}^{2}(f)=2\sigma_{f}^{2}s^{2}/f^{2}$ for our circular Gaussian PSF then gives

[TABLE]

where $\sigma_{f_{\textrm{ML}}}\equiv\sigma_{f}(f_{\textrm{ML}})$ and $f^{*}_{\textrm{ML}}=f_{\textrm{ML}}(x^{*},y^{*})$ are defined as in §3.

Since we showed §3 that the ML estimator at the true position $f_{\textrm{ML}}^{*}$ is unbiased, this explicitly demonstrates that the mean flux density, marginalized over position, is also unbiased.

Appendix G Second-Order Expansion in Flux Density

As discussed in Section §4.1, fixing the position of an object to the best-fit from a detection band $D$ leads to an underestimate in the flux density in all other bands. If the detected position is distributed as a 2-D Gaussian of width $\sigma_{x}$ about the true position, the average best-fit flux is given by Equation 43:

[TABLE]

At low signal to noise, the position error obtained by inverting the second derivative with respect to position (Equation 14) needs to be corrected by higher-order terms. The maximum-likelihood is found by setting the partial derivatives of the log-likelihood to zero. We Taylor expand the partial derivatives to leading order about the true parameters $(x^{*},y^{*},f^{*})$:

[TABLE]

with $\alpha,\beta\in\{x,y,f\}$. Evaluating the partial derivatives at the true parameters yields

[TABLE]

[TABLE]

with $Z_{i}=(\hat{f}_{i}-f^{*}p_{i}(x^{*},y^{*}))$ being independent, normally distributed variables with mean zero and unit variance. The sum over pixels of the second term of Equation G4 can be approximated in the oversampled limit with integrals over the PSF and its derivatives. Using Equations G3 and G4, Equation G2 can be written as a matrix equation:

[TABLE]

with $\zeta=\sum_{i}Z_{i}p_{i}(x^{*},y^{*})$, $\zeta_{\alpha}=\sum_{i}Z_{i}\partial_{\alpha}p_{i}(x^{*},y^{*})$, and $\zeta_{\alpha\beta}=\sum_{i}Z_{i}\partial_{\alpha\beta}p_{i}(x^{*},y^{*})$. Solving this matrix equation for $x_{\textrm{ML}}-x^{*}$ and expanding to second order in $\sigma$ yields:

[TABLE]

The sums $\zeta,\zeta_{\alpha},\zeta_{\alpha\beta}$ are random variables with mean zero just as $Z_{i}$ are. It can be shown that the expectation value of $x_{\textrm{ML}}-x^{*}$ is zero (ie. the position is unbiased) by evaluating expectation values like

[TABLE]

with the second equality following from the independence of the $Z_{i}$ and the third equality holding in the oversampled limit. Similarly,

[TABLE]

Since $x_{\textrm{ML}}-x^{*}$ has an expectation value of zero, the variance in position can be found by squaring Equation G6 and taking the expectation value. To evaluate the variance in position, these expectation values are needed:

[TABLE]

which yield:

[TABLE]

Approximating the position errors as being distributed as a 2-D Gaussian with this variance and using Equation 43 gives:

[TABLE]

Appendix H Gaussian Galaxies