Issue 
A&A
Volume 616, August 2018



Article Number  A95  
Number of page(s)  21  
Section  Celestial mechanics and astrometry  
DOI  https://doi.org/10.1051/00046361/201732537  
Published online  31 August 2018 
Optimality of the maximum likelihood estimator in astrometry
^{1}
Information and Decision Systems Group, Department of Electrical Engineering, Facultad de Ciencias Físicas y Matemáticas, Universidad de Chile, Beauchef 850, Santiago, Chile
email: sebastian.espinosa@ing.uchile.cl, josilva@ing.uchile.cl
^{2}
Departamento de Astronomía, Facultad de Ciencias Físicas y Matemáticas, Universidad de Chile, Casilla 36D, Santiago, Chile
email: rmendez@uchile.cl
^{3}
Department of Electrical Engineering, University of Southern California, 90007 CA, USA
Received:
22
December
2017
Accepted:
4
May
2018
Context. Astrometry relies on the precise measurement of the positions and motions of celestial objects. Driven by the everincreasing accuracy of astrometric measurements, it is important to critically assess the maximum precision that could be achieved with these observations.
Aims. The problem of astrometry is revisited from the perspective of analyzing the attainability of wellknown performance limits (the Cramér–Rao bound) for the estimation of the relative position of lightemitting (usually pointlike) sources on a chargecoupled device (CCD)like detector using commonly adopted estimators such as the weighted least squares and the maximum likelihood.
Methods. Novel technical results are presented to determine the performance of an estimator that corresponds to the solution of an optimization problem in the context of astrometry. Using these results we are able to place stringent bounds on the bias and the variance of the estimators in close form as a function of the data. We confirm these results through comparisons to numerical simulations under a broad range of realistic observing conditions.
Results. The maximum likelihood and the weighted least square estimators are analyzed. We confirm the suboptimality of the weighted least squares scheme from medium to high signaltonoise found in an earlier study for the (unweighted) least squares method. We find that the maximum likelihood estimator achieves optimal performance limits across a wide range of relevant observational conditions. Furthermore, from our results, we provide concrete insights for adopting an adaptive weighted least square estimator that can be regarded as a computationally efficient alternative to the optimal maximum likelihood solution.
Conclusions. We provide, for the first time, closeform analytical expressions that bound the bias and the variance of the weighted least square and maximum likelihood implicit estimators for astrometry using a Poissondriven detector. These expressions can be used to formally assess the precision attainable by these estimators in comparison with the minimum variance bound.
Key words: astrometry / methods: statistical / methods: analytical
© ESO 2018
1. Introduction
Astrometry, which deals with the accurate and precise measurement of the positions and motions of celestial objects, is the oldest branch of observational astronomy, dating back at least to Hipparchus of Nicaea in 190 BC. Since, from its very beginnings, this branch of astronomy has required measurements over time to fulfil its goals, it could be considered the precursor of the currently fashionable “timedomain astronomy”, preceding it by at least 20 centuries. In recent years, astrometry has experienced a “coming of age” motivated by the rapid increase in positional precision allowed by the use of alldigital techniques and space observatories^{1}.
A number of techniques have been proposed to estimate the location and flux of celestial sources as recorded on digital detectors, such as chargecoupled devices (CCD). In this context, estimators based on the use of a least squares (LS) error principle have been widely adopted (King 1983; Stetson 1987; Alard & Lupton 1998). The use of this type of decision rule has been traditionally justified through heuristic reasons. First, LS methods are conceptually straightforward to formulate based on the observation model of these problems. Second, they offer computationally efficient implementations and have shown reasonable performance (Lee & Van Altena 1983; Stone1989; Vakili & Hogg 2016). Finally, the LS approach was the classical method used when the observations were obtained with analog devices (Van Altena & Auer 1975; Auer & Van Altena 1978), which are well characterized by a Gaussian noise model for the observations. In the Gaussian case the LS is equivalent to the maximum likelihood (ML) solution (ChunLin 1993; Gray & Davisson 2004; Cover & Thomas 2012), and, consequently, the LS method was taken from the analog to the digital observational (Poisson noise model) setting.
Considering astrometry as an inference problem (of, usually, point sources), the astrometric community has been interested for a long time in understanding the fundamental performance limits (or information bounds) of this task (Lindegren 2010 and references therein). It is well understood by the community that the characterization of this precision limit offers the possibility of understanding the complexity of the task and how it depends on key attributes of the problem, like the quality of the observational site, the performance of the instrument (CCD), and the details of the experimental conditions (Mendez et al. 2013, 2014). On the other hand, it provides meaningful benchmarks to define the optimality of practical estimators in the process of comparing their performance with the bounds (Lobos et al. 2015).
Concerning the characterization and analysis of fundamental performance bounds, we can mention some works on the use of the parametric Cramér–Rao (CR) bound by Lindegren (1978), Jakobsen et al. (1992), Zaccheo et al. (1995), Adorf (1996), and Bastian (2004). The CR bound is a minimum variance bound (MVB) for the family of unbiased estimators (Rao, 1945; Cramér, 1946). In astrometry and joint photometry and astrometry, Mendez et al. (2013, 2014) have recently studied the structure of this bound, and have analyzed its dependency with respect to important observational parameters under realistic (groundbased) astronomical observing conditions. In this context, closedform expressions for the Cramér–Rao bound were derived in a number of important settings (high pixel resolution and low and high signaltonoise (S/N) regimes), and their trends were explored across different CCD pixel resolutions and the position of the object in the CCD array. As an interesting outcome of those studies, the analysis of the CR bound has allowed us to predict the optimal pixel resolution of the array, as well as providing a formal justification for some heuristic techniques commonly used to improve performance in astrometry, like dithering for undersampled images (Mendez et al. 2013, Sect. 3.3). Recently, an application of the CR bound to moving sources has been done by Bouquillon et al. (2017), indicating excellent agreement between our theoretical predictions, simulations, and actual groundbased observations of the Gaia satellite, in the context of the Gaia Ground Based Optical Tracking (GBOT) program (Altmann et al. 2014). The use of the CR bound on other applications is also of interest, for example, in assessing the performance of star trackers to guide satellites with demanding pointing constraints (Zhang et al. 2016), or to meaningfully compare positional differences from different catalogs (for example involving the Sloan Digital Sky Survey (SDSS) and Gaia, see Lemon et al. 2017). Finally, a formulation of the (nonparametric) Bayesian CR bound in astrometry, using the socalled “VanTrees inequality” (Van Trees, 2004), has been presented by our group in Echeverria et al. (2016). This approach is particularly well suited for objects at the edge of detectability, and where some prior information is available, and has been proposed for the analysis of Gaia data for faint sources, or for those with a poor observational history (Michalik et al. 2015; Michalik & Lindegren, 2016).
From the perspective of astrometric estimators, Lobos et al. (2015) have studied in detail the performance of the widely adopted LS estimator. In particular extending the result in So et al. (2013), Lobos and collaborators derived lower and upper bounds for the mean square error (MSE) of the LS estimator. Using these bounds, the optimality of the LS estimator was analyzed, demonstrating that for high S/N there is a considerable gap between the CR bound and the performance of the LS estimator (indicating a lack of optimality of this estimator). This work showed that for the very low S/N observational regime (weak astronomical sources), the LS estimator is near optimal, as its performance closely follows the CR bound. The limitations of the LS method in the medium to high S/N regime proved in that work opens up the question of studying alternative estimators that could achieve the CR bound on these regimes, which is the main focus of this paper, as outlined below.
2. Contribution and organization
In this work we study the ML estimator in astrometry, motivated by its wellknown optimality properties in a classical parametric estimation setting with independent and identically distributed measurements (i.i.d.; Kendall et al. 1999). We know that in the i.i.d. case this estimator is efficient with respect to the CR limit (Kay 1993), but it is important to emphasize that the observational setting of astrometry deviates from the classical i.i.d. case and, consequently, the analysis of its optimality is still an open problem. In particular, we face the technical challenge of evaluating its performance, a problem that, to the best of our knowledge, has not been addressed by the astrometric community. Concerning the independent but not identically distributed case, Bradley & Gart (1962) and Hoadley (1971) give conditions under which ML estimators are consistent^{2} and asymptotically normal^{3}. Those conditions, however, are technically difficult to proof in this context.
The main challenge here is the fact that, as in the case of the LS estimator (Lobos et al. 2015), the ML estimator is the solution of an optimization problem with a nonconvex cost function of the data. This implies that it is not possible to directly compute the performance of the method. To address this technical issue, we extend the approach proposed by Fessler (1996) to approximate the variance and the mean of an implicit estimator solution of a generic optimization problem of the data through the use of a Taylor approximation around the mean measurement (see Theorem ^{1} below). Our extension considers high order approximations of the function that allow us not only to estimate the performance of the ML estimator through an explicit nominal value, but also provides a confidence interval around it. With this result we revisit the more general weighted least square (WLS) and ML methods providing specific upper and lower bounds for both methods (see Theorems ^{2} and ^{3}). The main findings from our analysis of the bounds are twofold: first we show that the WLS exhibits a suboptimality similar to that of the LS method for medium to high S/N regimes discovered by Lobos et al. (2015) and, second, that the ML estimator achieves the CR limit for medium to high S/N and, consequently, it is optimal on those regimes. This last result is remarkable because, in conjunction with the result presented in Lobos et al. (2015), we are able to identify estimators that achieve the fundamental performance limits of astrometry in all the S/N regimes for the problem.
The paper is organized as follows: Sect. 3 introduces the background, preliminaries, and notation of the problem. Section 4 presents the main methodological contribution of this work. Section 5 presents the application to astrometry considering the WLS and the ML schemes. Finally, numerical analysis of these performance bounds are presented in Sect. 6 and the final remarks and conclusion are given in Sect. 7.
3. Preliminaries and background
We begin by introducing the problem of astrometry. For simplicity, we focus on the onedimensional (1D) scenario of a linear array detector, as it captures the key conceptual elements of the problem^{4}.
3.1. Relative astrometry as a parameter estimation problem
The main problem at hand is the inference of the relative position (in the array) of a point source. This source is modeled by two scalar quantities, the position of object x _{c} ∈ ℝ in the array^{5}, and its intensity (or brightness, or flux) that we denote by F̃ ∈ ℝ^{+} . These two parameters induce a probability distribution over an observation space that we denote by 𝕏. Formally, given a point source represented by the pair (x _{c}, F̃), it creates a nominal intensity profile in a photon integrating device (PID), typically a CCD, which can be expressed by(1)
where ϕ(x−x _{c},σ) denotes the onedimensional normalized point spread function (PSF) and where σ is a generic parameter that determines the width (or spread) of the light distribution on the detector (typically a function of wavelength and the quality of the observing site, see Sect. 6; Mendez et al. 2013, 2014).
The profile in Eq. (1) is not measured directly, but it is observed through three sources of perturbations. First, an additive background noise that captures the photon emissions of the open (diffuse) sky, and the noise of the instrument itself (the readout noise and darkcurrent Janesick 2001, 2007; Howell 2006; McLean 2008), modeled by B̃ _{i} in Eq. (2). Second, an intrinsic uncertainty between the aggregated intensity (the nominal object brightness plus the background) and actual measurements, which is modeled by independent random variables that follow a Poisson probability law. Finally, we need to account for the spatial quantization process associated with the pixel resolution of the PID as specified in Eqs. (2) and (3). Modeling these effects, we have a countable collection of independent and nonidentically distributed random variables (observations or counts) {I _{i} : i ∈ ℤ, where I _{i} ∼ Poisson(λ _{i}(x _{c}, F̃)), driven by the expected intensity at each pixel element i, given by(2)
where 𝔼{} is the expectation value of the argument and {x _{i} : i ∈ ℤ} denotes the standard uniform quantization of the real line array with resolution ∆x > 0, that is, x _{ i+1} − x _{i} = Δx for all i ∈ ℤ. In practice, the PID has a finite collection of measurement elements (or pixels) I _{1},…,I_{n}, then a basic assumption here is that we have a good coverage of the object of interest, in the sense that for a given position x _{c} (4)
At the end, the likelihood (probability) of the joint observations I^{n} = (I_{1},…,I_{n}) (with values in ℕ^{ n }) given the source parameters (x _{c}, F̃) is given by(5)
where denotes the probability mass function (PMF) of the Poisson law (Gray & Davisson, 2004).
Finally, if F̃ is assumed to be known^{6}, the astrometric estimation is the task of defining a decision rule τ_{n}() : ℕ^{n} → Θ with Θ = ℝ being the parameter space, where given an observation I ^{n} the estimated position is given by x̂_{c}(I^{n} ) = τ_{n} (I^{n} ).
3.2. The Cramér–Rao bound
In astrometry the Cramér–Rao bound has been used to bound the variance (estimation error) of any unbiased estimator (Mendez et al. 2013, 2014). In general, let I^{n} be a collection of independent observations that follow a parametric PMF defined on ℕ. The parameters to be estimated from I^{n} will be denoted in general by the vector . Let τ_{n}(I^{n} ) : ℕ^{ n } → Θ be an unbiased estimator^{7} of , and be the likelihood of the observation given . Then, the Cramér–Rao bound (Rao 1945; Cramér 1946) establishes that if(6)
then, the variance (denoted by Var), satisfies that(7)
where is the Fisher information matrix given by, ∈i, j ∈ {1,…,m},(8)
In particular, for the scalar case (m = 1), we have that for all θ ∈ Θ, (9)
where is the collection of all unbiased estimators and . For astrometry, Mendez et al. (2013, 2014) have characterized and analyzed the Cramér–Rao bound, leading to
3.3. Achievability and performance of the LS estimator
Concerning the achievability of the CR bound with a practical estimator, Lobos et al. (2015, Proposition ^{2}) have demonstrated that this bound cannot be attained, meaning that for any unbiased estimator τ_{n} (·) we have that(11)
where I^{n} follows the Poisson PMF in Eq. (5).
This finding should be interpreted with caution, considering its pure theoretical meaning. This is because Eq. (11) does not exclude the possibility that the CR bound could be approximated arbitrarily closely by a practical estimation scheme. Motivated by this refined conjecture, Lobos et al. (2015) proposed to study the performance of the widely adopted LS estimator^{8} with the goal of deriving operational upper and lower performance bounds of its performance that could be used to determine how far this scheme could depart from the CR limit. Then, from this result, it was possible to evaluate the goodness of the LS estimator for concrete observational regimes. To bound the performance of the LS estimator, the challenge was that τ _{LS}(I^{n} ) is an implicit function of the data (where no closeform expression is available) and, consequently, Lobos et al. (2015) derived a result to bound the estimation error and the variance of τ _{LS}(I^{n} ). We can briefly summarize the main result presented in Lobos et al. (2015), Theorem ^{1}) by saying that under certain mild sufficient conditions (that were shown to be realistic for astrometry), there is a constant δ > 0 (that depends on the observational regime, in particular the S/N) and a nominal variance , which is determined in closedform in the result, from which it is possible to bound (τ _{LS}(I^{n} ) by the simple expression(12)
We note that when δ is small, tightly determines the performance of the LS estimator, and its comparison with can be used to evaluate the goodness of the LS estimator for astrometry. Based on a careful comparison, it was shown in Lobos et al. (2015, Sect. 4) that, in general, is close to for the small S/N regime of the problem. However, for moderate to hig h S/N regimes, the gap between and becomes quite significant^{9}.
These unfavorable findings for the LS method have motivated us to study alternative schemes that could potentially approach better the Cramér–Rao bound for the rich observational context of medium to high S/N regimes. This will be the focus of the following sections, where in particular we explore the performance of the ML and WLS estimators, thus extending and generalizing the analysis done for the LS estimator by our group presented in Lobos et al. (2015).
4. Bounding the performance of an implicit estimator
Before we go to the case of the WLS and the ML estimators, we present a general result that bounds the performance of any estimator that is the solution of a generic optimization problem. Let us consider a vector of observations I^{n} = (I_{1},…I_{n}) ∈ ℝ^{n} and a general socalled cost function J(α, I^{n} ). Then the estimation of x_{c} from the data is the solution of the optimization problem(14)
where α represents the position of the object in the context of astrometry. As in our previous work (Lobos et al. 2015), the challenge here is that this estimator is implicit because no closedform expression of the data which solves Eq. (14) is assumed. In particular, this implies that both the variance and the estimation error of τ_{J}(I^{n}) cannot be determined directly. To address this technical issue, we extend the approach proposed by Fessler (1996) to approximate the variance and the mean of an implicit estimator solution of a problem described by Eq. (14) through the use of a Taylor approximation around the mean measurement, that is, .
More precisely, we assume that J(α,I^{n} ) has a unique global minimum at τ_{J} (I^{n} ), and that it has a regular behavior, so its partial derivatives are zero, that is,(15)
Then we can obtain τ_{J} (I^{n} ), by a first order Taylor expansion around the mean by(16) (17)
where t ∈ [0,1] is fixed but unknown^{10}. For simplicity, Eq. (17) can be written in matrix form as(18)
where denotes the row gradient operator and is the residual error of the Taylor expansion. From Eq. (18) we can readily obtain the expression for its variance(19) (20) (21)
In Eq. (21) we recognize two terms: that captures the linear behavior of τ_{J} (·) around and γ_{J} (n) which reflects the deviation from this linear trend. It should be noted that the above expression does not depend on τ(I^{n}) itself, but on its partial derivatives evaluated at the mean vector of observations. Then in the adoption of this approach to estimate Var{τ_{J} (I^{n} )}, a key task is to determine .
In an effort to analyze both the linear and nonlinear aspects of a general intrinsic estimator solution to Eq. (14), the following result offers sufficient conditions to determine in closedform, and to bound the magnitude of the residual term γ_{J} (n) in Eq. (21).
The proof of this result and the expression for in Eqs. (24) and (25) are presented in detail in Appendix A.
Revisiting the equality in Eq. (21), Theorem ^{1} provides general sufficient conditions to bound the residual term γ_{J} (n) and by doing that, a way of bounding the variance of τ_{J} (I) which is the solution of Eq. (14). In particular, it is worth noting that if the ratio , then Eq. (25) offers a tight bound for . In this last context, (called the nominal value of the result) provides a very good approximation for .
On the application of this result to the WLS and ML estimators, we will see that the main assumption in Eq. (22) is satisfied in both cases (see Eqs. (B.8) and (C.6)), and from that plays an important role to approximate the performance of ML and WLS in a wide range of observational regimes. In addition, the analysis of the bias in Eq. (24) shows that these estimators are unbiased for any practical purpose and, consequently, contrasting their performance (estimation error ) with the CR bound is a meaningful way to evaluate optimality.
5. Application to astrometry
In this section we apply Theorem ^{1} to bound the variances of the ML and WLS estimators in the context of astrometry. Following the model presented in Sect. 3.1, I^{n} = (I_{1},...,I_{n})^{T} denotes the measurements acquired by each pixel of the array and where each of them follows a Poisson distribution given by(26)
as expressed by Eqs. (2) and (3).
5.1. Bounding the variance of the WLS estimator
The WLS estimator, denoted by τ_{WLS}(I^{n}) in Eq. (28), is implicitly defined through a cost function given by(27)
where is a weight vector, and α is a general source position parameter. Specifically we have that(28)
Applying Theorem ^{1} we obtain the following result:
The proof of this result and the expressions for ϵ_{WLS}(n) and β_{WLS}(n) in Eqs. (29) and (30), respectively, are elaborated in Appendix B. This result offers concrete expressions to bound the bias as well as the variance of the WLS estimator. For the bias bound in Eq. (29), it will be shown that ϵ_{WLS}(n) is very small (orders of magnitude smaller than x_{c}) for all the observational regimes explored in this work and, consequently, the WLS can be considered an unbiased estimator in astrometry, as would be expected. Concerning the bounds for the variance in Eq. (30), we will show that for high and moderate S/N regimes the ratio , and consequently in this context is a precise estimator of . For the very small S/N the results offers an admissible interval around the nominal value . Therefore in any context shows to be a meaningful approximation for the performance of the WLS.
5.2. Bounding the variance of the ML estimator
The ML estimator, denoted by τ_{WLS}(I^{n}) in Eq. (34), is implicitly defined through a cost function(33)
where α is a general source position parameter. Specifically, given an observation I^{n} we have that(34)
Applying Theorem ^{1} we obtain the following result:
The proof of this result and the expressions for ϵ_{ML}(n) and β_{ML}(n) in Eqs. (35) and (36), respectively, are elaborated in Appendix C.
6. Numerical analysis
In this section we evaluate numerically the performance bounds obtained in Sect. 5 for the WLS and ML estimators, and compare them with the astrometric CR bound in Proposition ^{1}. The idea is to consider some realistic astrometric conditions to evaluate the expressions developed in Theorems ^{2} and ^{3} and their dependency on important observational conditions and regimes. As we shall see, key variables in this analysis are the tradeoff between the intensity of the object and the noise represented by the S/N value, and the pixel resolution of the CCD.
6.1. Experimental setting
We adopt some realistic design and observing variables to model the problem (Mendez et al. 2013, 2014). For the PSF, analytical and semiempirical forms have been introduced (see for instance the groundbased model in King 1971 and the spacebased models by King 1983 or Bendinelli et al. 1987). In this work we will adopt a Gaussian PSF, that is, in Eq. (3), and where σ is the width of the PSF and is assumed to be known. This PSF has been found to be a good representation for typical astrometricquality groundbased data (Méndez et al. 2010). In terms of nomenclature, measured in arcsec, denotes the fullwidth at halfmaximum parameter, which is an overall indicator of the image quality at the observing site (Chromey, 2016).
The background profile, represented by {B̃_{i}, i = 1,…,n} in Eq. (2), is a function of several variables, like the wavelength of the observations, the moon phase (which contributes significantly to the diffuse sky background), the quality of the observing site, and the specifications of the instrument itself. We will consider a uniform background across pixels underneath the PSF, that is, B̃_{i} = B̃ for all i. To characterize the magnitude of B̃, it is important to first mention that the detector does not measure photon counts [e^{−}] directly, but a discrete variable in “Analog to Digital Units (ADUs)” of the instrument, which is a linear proportion of the photon counts (McLean 2008). This linear proportion is characterized by the gain of the instrument G (scaling value) in units of [e^{−}/ADU]. We can define F ≡ F̃/G and B ≡ B̃/G as the intensity of the object and noise, respectively, in the specific ADUs of the instrument. Then, the background (in ADUs) depends on the pixel size Δx arcsec as follows:(38)
where f_{s} is the (diffuse) sky background in ADU arcsec^{−1}, while D and RON^{2}, both measured in e^{−}, model the darkcurrent and readoutnoise of the detector on each pixel, respectively. The first component in Eq. (38) is attributed to the site, and its effect is proportional to the pixel size. On the other hand, the second component is attributed to errors of the PID (detector), and it is pixelsize independent. This distinction is central when analyzing the performance as a function of the pixel resolution of the array (see details in Mendez et al. 2013, Sect. 4). More important is the fact that in typical groundbased astronomical observation, long exposure times are considered, which implies that the background is dominated by diffuse light coming from the sky, and not from the detector Mendez et al. (2013), Sect. 4).
For the experimental conditions, we consider the scenario of a groundbased station located at a good site with clear atmospheric conditions and the specification of current sciencegrade CCDs, where f_{s} = 1502.5 ADU arcsec^{−1}, D = 0, RON = 5 e^{−} arcsec (equivalent to arcsec), and G = 2, e^{−} ADUe^{−1} (with these values we will have B = 313 ADU for Δx = 0.2 arcsec using Eq. (38)). In terms of scenarios of analysis, we explore different pixel resolutions for the CCD array Δx ∈ [0.1,0.65] measured in arcsec, and different signal strengths F̃ ∈ {1080; 3224; 20 004; 60 160}, measured in e^{−}, which corresponds to S/N ∈ ∼ {12,32,120,230}. Increasing F̃ implies increasing the S/N of the problem, which can be approximately measured by the ratio F̃/B̃. On a given detector plus telescope setting, these different S/N scenarios can be obtained by appropriately changing the exposure time (open shutter) that generates the image.
Fig. 1.
Relative performance of the bias (as measured by ) stipulated by Theorem ^{1} for the WLS estimator (left side, Eq. (29)) and the ML estimator (right side, Eq. (35)). Results are reported for different values of the source flux F̃ ∈ {1080, 3224, 20 004, 60, 160}, all in e^{−} (top to bottom symbols respectively), as a function of the detector pixel size. The 0% level corresponds to having achieved no bias. 
6.2. Bias analysis
Considering the upper bound terms ϵ_{WLS}(n) and ϵ_{ML}(n) for the bias error obtained from Theorems ^{2} and ^{3} for the WLS and ML, respectively, Fig. 1 presents the relative bias error for different S/N regimes and pixel resolutions. In the case of the ML estimator, the bounds for relative bias error are very small in all the explored S/N regimes and pixel resolutions meaning that for any practical purposes this estimator is unbiased as expected from theory (Gray & Davisson, 2004). For the case of the WLS, we observe that from medium to high S/N the relative error bound obtained is very small, meanwhile at low S/N unbiasedness cannot be fully guaranteed from the bound in Eq. (29). In general, our results show that both WLS and ML are unbiased estimators for astrometry in a wide range of relevant observational regimes (in particular from medium to high S/N) and, consequently, it is meaningful to analyze the optimality of these estimators in comparison with the CR bound in those regimes.
In the following sections, we move to the analysis of the variance of the WLS and ML with particular focus on the medium to high S/N regimes across all pixel resolutions using the performance bounds derived in Eqs. (30) and (36), respectively.
6.3. Performance analysis of the WLS estimator
In this section, we evaluate numerically the expression derived in Theorem ^{2} to bound the variance of the WLS estimator in Eq. (30). For that we characterize the admissible regime predicted for the variance of the WLS estimator, that is, the interval(39)
for S/N ∈ ∼ {12,32,120,230} and Δx ∈ [0.01,0.65] arcsec. In these bounds, we recognize its central value (or nominal value) in Eq. (31) and the length of the interval 2β_{WLS}(n) that is determined in closed form for its numerical evaluation in Eqs. (A.2) and (B.11). We note that 2β_{WLS}(n) can be considered an indicator of the precision of our result to approximate the variance of the WLS in astrometry.
Fig. 2.
Range of the square root of the variance performance (in miliarcsecond = mas) for the WLS method in astrometry using uniform weights (equivalent to the LS method) predicted by Theorems ^{1} and ^{2}, Eq. (30). Results are reported for different representative values of F̃ and across different pixel sizes (topleft to bottomright): F̃ ∈ {1080; 3224; 20 004; 60 160} e^{−}. 
6.3.1. Revisiting the uniform weight case
To begin the analysis, we consider the setting of uniform weights across pixels, that is, the case of the LS estimator and, without loss of generality, we locate the object in the center of the field of view^{11}, which can be considered a reasonable scenario to represent the complexity of the astrometry task. At this point, it is important to remind the reader that from the analysis of in Remark ^{2}, the nominal value is equivalent to the CR bound when the w_{i} are selected as a function of the true position in Eq. (32). In view of this observation, the selection of nonuniform weights can be interpreted as biasing the estimation to a particular area of the angular space, which goes in contradiction to the essence of the inference problem that estimates the position with no prior information, and only relies on the measured counts. From this interpretation, revisiting the LS estimator is an important first step in the analysis of the WLS framework.
On the specifics, the boundaries of the interval in Eq. (39) and its nominal values are illustrated in Fig. 2 for the different observational regimes. In addition, Fig. 2 shows the CR bound as a reference to evaluate the optimality of the LS scheme across settings. We observe that for the low S/N ∼ 12 regime, the nominal values precisely match the CR bound, however, our result is not conclusive as the interval around the nominal performance is significantly large. This is the regime where our result is not conclusive regarding the performance of the LS estimator. In the regime S/N ∈ (30,50) (top right panel on Fig. 2), we notice an important reduction in the range of admissible performance, and our result becomes more informative and meaningful. In this context, the nominal values are very close to the CR bound, and we could assert that the LS estimator offers sufficiently good performance in the sense that it is very competitive with the MVB. Importantly, when we move to the regime of relatively high S/N and very high S/N (from 100 to 300), our results are very accurate in predicting the performance of the LS method, and we find that the gap between the CR and the nominal value is very significant (the deviation from the MVB is 16% and 30% for S/N 120 and 230 at Δx = 0.2 arcsec, respectively). This last result confirms one of the main findings presented in Lobos et al. (2015), who showed that for medium to high S/N the LS estimator is suboptimal with respect to the MVB. In Fig. 2 we also show the square root of the empirical variance () with respect to the empiricallydetermined mean position x̂_{c} (using the WLS estimator), all as derived from the simulations, showing good consistency with our predictions.
Fig. 3.
Worstcase discrepancies in Eq. (40) for the WLS estimator using the weights set indexed by the positions . Results are reported for two S/N scenarios, namely F̃ = 20 004 e^{−} (left) and F̃ = 60 160 e^{−} (right) (right), and across different pixel sizes. 
6.3.2. Nonuniform weight case
The suboptimality of the WLS scheme from moderate to high S/N seen in Fig. 2 can be extended for any arbitrary selection of a fixed set of weights, as would be expected. Given that the space of weights’ selection is literally unlimited, we use the insight obtained from Remark ^{2} that states that a selection of weights can be interpreted as an specific prior on the position of the object where the optimum, but unfortunately unknown, selection (achieving the CR bound) of weights in Eq. (32) is an explicit function of the unknown position of the object. Then, we consider a finite set of positions {x_{c,1},…,x_{c,M}} that uniformly partition the field of view, and their corresponding weight sets {w_{i}(x_{c,1}) : 1 = 1,…,n}, {w_{i}(x_{c,1}) : 1 = 1,…,n}, {w_{i}(x_{c,M}) : 1 = 1,…,n} using Eq. (32) to cover a representative collections of weights for the problem of astrometry.
Then, for a specific selection of weights in our admissible collection (attributed to a prior belief regarding the position of the object in the field of view), we evaluate the worse discrepancy between (which is the most favorable expression for the variance predicted from Theorem ^{2}), and the CR bound in Proposition ^{1}, across a collection of presumed object positions in the range of positions
where denotes the center of the array (which, as indicated at the beginning of Sect. 6.3.1, is equal to the true object position x_{c}) and is the dispersion parameter of the PSF. The idea of using this worstcase difference is justified from the fact that in this parameter estimation problem we do not know the position of the object and, consequently, the optimality of any WLS estimator should be evaluated in the worsecase situation, as the scheme should be able to estimate the position of the object in any scenario (position). More precisely, for a given weight selection {w_{i}, i=1,…,n}, we use the worstcase discrepancy(40)
For this analysis, we note that both and are functions of the position x_{c}, Δx, and S/N.
Figure 3 illustrates the worstcase discrepancy in Eq. (40) for the medium and high S/N regimes where Theorem ^{2} provides an accurate and meaningful prediction of the performance of the WLS method, that is, for S/N ∈ {120,230}, and across Δx ∈ [0.05,0.7] arcsec. The discrepancy is quite significant, on the order of 37% and 60% in the range for Δx ∈ [0.1,0.3] arcsec for S/N 120 and 230, respectively.
To refine the worstcase analysis presented in Fig. 3, and to evaluate in more detail the sensitivity of the discrepancy indicator given by Eq. (40), we evaluate the discrepancy of WLS using the weights associated to (the center position of the array) with respect to the CR bound associated to the positions , to study how the discrepancy (measuring the nonoptimality of the method) evolves when the adopted position moves far from the prior imposed by the WLS in the center of the array. Figure 4 illustrates this behavior, where it is possible to see that the discrepancy is very sensitive and proportional to the misassumption of the object position, where the worstcase discrepancy in the maximum achievable location precision is on the order of 40% for pixel sizes in the range [0.1, 0.6] arcsec for S / N ~ 120, and about 60% for pixel sizes in the range [0.1, 0.6] arcsec for S / N ~ 230. These worstcase scenarios happens in both cases when the object is located the farthest from the prior assumption, as would be expected.
The main conclusion derived form this CR bound analysis is that, independent of the weight selection adopted, as long as the weights are fixed, the WLS estimator is not able to achieve the CR bound in all observational regimes. More precisely, the discrepancy (measuring the nonoptimality) in the less favorable case of an hypothetical and feasible position of the object is very significant, in the range of 40%–60% for the important regime of high and very high S/N.
Fig. 4.
Performance discrepancies (measuring the nonoptimality) of the WLS estimator using the center position as a prior for the weight selection with respect to the CR bound obtained for the true object positions . Results are reported for two S/N scenarios, namely F̃ = 20 004 e^{−} (left) and F̃ = 60 160 e^{−} (right), and across different pixel sizes. 
6.4. Performance analysis of the ML estimator
In this section we perform the same analysis done for the WLS in Sect. 6.3, but using the result in Theorem ^{3}. In particular, we consider the admissible regime for the variance of the ML estimator given by
in Eq. (36), where the nominal value in this case, in Eq. (37), is precisely the CR bound, while the length of the interval 2β_{ML}(n), is given by Eqs. (A.2) and (C.9) in closed form.
Considering an object located in the center of the array, that is, x_{c} = x_{0} the performance curves are presented in Fig. 5 for S/N ∈ {12, 32, 120, 230} and Δx ∈ [0.1, 0.65] arcsec. First, we note that there is a significant difference in the predictions of our methodology for the ML estimator in comparison with what we predict in the WLS case. In fact, the results of our approach are very precise for the determination of the variance of the ML estimator in all the regimes, from small to high S/N, which is remarkable. More important it is the fact that, from these findings, we observe that the performance deviation from the MVB in the worst case (small S/N) is very small (see Table 1, first row), while for any practical purposes the variance of the ML estimator achieves the CR limit for all the other regimes, from medium to high S/N, which is a numerical corroboration of the optimality of the ML estimator in astrometry, as predicted theoretically by Theorems ^{1} and ^{3}. In Fig. 5 we also show the square root of the empirical variance (Var(τM̂L)) with respect to the empiricallydetermined mean position x̂_{c} (using the ML estimator), all as derived from the simulations, showing good consistency with our predictions.
Complementing this analysis, we conducted the same comparison considering different positions for the object within the array obtaining the same trends and conclusions. To summarize these results, Fig. 6 shows the value , which is an indicator of the quality of the estimator (the smaller the better) for different scenarios of the position of the object . In particular, for all the evaluated positions, the relative discrepancy is bounded (relative to the CR bound) by values in the range of 0.025% and 0.012% for pixel resolution in the range Δx ∈ [0.1, 0.2] arcsec for S/N = 120 and S/N = 230, respectively. Finally, Table 1 presents the relative discrepancy for all the range of S/N values considered in this study for the case Δx = 0.2 arcsec.
We conclude from this analysis that the ML estimator in nearly optimal for the medium, high, and very high S/N regimes across pixel resolutions, achieving the MVB for astrometry. This is an interesting result, since it lends further support to the adoption of these types of estimators for very demanding astrometric applications, as has been done in the case of Gaia (Lindegren, 2008). We note that Vakili & Hogg (2016) reach the same conclusion regarding the optimality of the ML method in comparison with the MVB, through simulations of twodimensional (2D) CCD images using a broad set of Moffat PSF stellar profiles. While their results are purely empirical, it is interesting that they test the ML using a different PSF from ours, and in a 2D scenario, and yet they reach the same conclusions. More recently, Gai et al. (2017) have tested (also empirically) the ML method in a 1D scenario (similar to ours), but in the context of a Gaialike PSF. They find that the ML is unbiased (in agreement with our results, see Fig. 1) and by comparing two implementations of the ML they conclude that they predict selfconsistent and reliable results over a broad range of flux, background, and instrument response variations. It would still be interesting to compare the performance of those implementations against the CR MVB in order to further test our theoretical predictions.
Fig. 5.
Range of the square root of the variance performance (in miliarcsecond = mas) for the ML method in astrometry as predicted by Theorems ^{1} and ^{3}, Eq. (36). Results are reported for different representative values of F̃ and across different pixel sizes (topleft to bottomright): F̃ ∈ {1080, 3224, 20 004, 60 160} e^{‒}. 
Performance quality of the ML estimator relative to the Cramér–Rao bound expressed in terms of the indicator from the result in Theorem ^{3}.
6.5. Comments on an adaptive WLS estimator for astrometry
In Sect. 5.1 we presented results that offer a nominal prediction for the variance of the WLS method through Eq. (30), which turns out to be very accurate in the regime from medium to high S/N as shown in Sect. 6.3. Interestingly, Remark ^{2} tells us that these nominal values precisely match the CR limit for an optimal selection of weights given in Eq. (32), namely w_{i} ∼ 1/λ_{i}(x_{c}) for all i = 1, . . . , n (compare Eqs. (31) and (37)). As this selection is unfeasible, because it requires the knowledge of x_{c} (see the expression in Eq. (2)), we can approximate this value by a noisy version of it, considering the fact that the expected value of the observations I_{i} that we measure at pixel i is λ_{i}(x_{c}) using our model in Eq. (4). Therefore I_{i} can be interpreted as a noisy version of λ_{i}(x_{c}) and(41)
can be seen as a noisy version of the ideal weight . Adopting this datadriven weighting approach, we would have an adaptive WLS method as the weights are not fixed but instead they are a function of the data {I_{1} : i = 1, . . . , n}. This selection of weights can be interpreted as an empirical version of the optimal weights that achieves the CR bound. Then, the problem reduces to solve(42)
Figure 7 presents the performance of this scheme for the same regimes we have been exploring in this work, supporting the conjecture that this selection of weights resembles the optimal weight selection and in fact achieves MSE performances that are surprisingly close to the CR bound in all the observational regimes.
Our results in this subsection show that some commonly adopted weighting schemes, using the analogous idea to Eq. (41), are a very good datadriven choice for methods that employ a WLS scheme, such is the case, for example, of the wellknow PSF stellar fitting program (including astrometry) Dominion Astrophysical Observatory Photometry (DAOPHOT), described in Stetson (1987, Eq. (10)).
Fig. 6.
Performance optimality of the ML estimator (computed as ) for different positions of the target object in the array, as a function of pixel resolution. The left panel shows the case Fį = 20 004 e^{‒}, the right panel the case Fį = 60 160 e^{‒}. 
Fig. 7.
Performance comparison between the of the adaptive WLS estimator and σ_{CR}(n), both in mas. Results are reported for different F̂ and across different pixel sizes: (topleft to bottomright) F̃ ∈ {1080, 3224, 60 160}e^{‒}. 
6.6. Undersampled case
So far we have analyzed the performance of the ML and WLS estimators in regimes where the PSF (assumed to have a FWHM = 1 arcsec) is well sampled by the detector. In this section we analyze the undersampled case, where we have only a few pixels that capture most the flux of the source. This scenario would be typical of nearly diffractionlimited images taken, for example, from space, or with adaptive optics (AO) system from the ground.
Table 2 shows the performance of the bounds presented in Theorems ^{2} and ^{3} in the undersampled regime, where the PSF is concentrated in only a few pixels (order of five to nine pixels), and where the flux is relatively large. For the spacebased case, the background is reduced almost to the pure instrument noise (left numbers in the table), while in the groundbased case the background would be larger, but we would still have a very large Strehl ratio. From the table it can be seen that the performance is nearly optimal with the ML estimator, even in the severely undersampled regime (small FWHM values). On the other hand, the WLS deteriorates its performance very quickly as the image size decreases, reaching up to nearly 45% when the FWHM is on the order of the detector pixel size.
We emphasize that the above results are computed in the ideal situation of a uniform pixel response function (and a perfect flatfield correction), whereas in reality a nonuniform pixel response has a critical impact on the Cramér–Rao bound as demonstrated by Adorf (1996), Fig. 1). In addition, in the case of staremode spacebased images acquired under very small background, charge transfer inefficiency effects due to timevarying charge traps in the detector would also blur (and systematically shift) the images, deteriorating the ultimate astrometric precision above the Cramér–Rao limit (for details see, e.g., Bristow et al. 2006, especially their Fig. 10).
7. Summary, conclusions, and outlook
In this paper we study the performance of the WLS and ML estimators for relative astrometry on digital detectors subject to Poisson noise, in comparison with the best possible attainable precision given by the CR bound. Our study includes analytical results and numerical simulations under realistic observational conditions to help us to corroborate our theoretical findings.
By generalizing the proposal by Fessler (1996) we are able to obtain, for the first time, closeform expressions for the variance and the mean of implicit estimators (as is the particular case of the WLS and ML schemes), which can be computed directly from the data (see Theorem ^{1}, in particular Eqs. (24) and (25), and Appendix A). When applying this result to astrometry with digital detectors, we are able to bound both the bias and the variance of the relative position of a celestial source on a CCD array as a function of all the relevant parameters of the problem (see Eqs. (29) and (30) or Eqs. (35) and (36) for the WLS and ML estimators, respectively). We verified that the bias of the WLS and ML methods are negligible in all the observational regimes explored in this paper (see Fig. 1).
A careful analysis of our predictions confirms earlier results by (Lobos et al. 2015; for the LS method) in that the WLS method is, in general, suboptimal (in comparison with the MVB given by the CR result), especially at high and very high S/N (see the two bottom panels on Fig. 2). However, a judicious datadriven selection of weights (called “adaptive” WLS method by us, Sect. 6.5), improves the performance of the WLS substantially (see Fig. 7). This is an interesting result, given the widespread use and simple numerical implementation of the WLS method.
The ML method is found to have both a smaller bias than the WLS method (compare left and right panels of Fig. 1), although the bias on both methods is already quite small, and a tight correspondence to the MVB throughout the entire range of S/N regimes explored in this paper (Fig. 5). Therefore, the ML estimator for astrometry is consistently optimal, and should be the estimator of choice for highprecision applications. This paper, along with Lobos et al. (2015), completes an indepth study of the performance of commonly used estimators in astrometry using PIDs, and sets the stage for the development of codes that could efficiently implement astrometric ML estimators on 2D detectors, incorporating also the simultaneous measurement of fluxes, as explored by Gai et al. (2017).
Performance quality of the ML and WLS estimators relative to the nominal bound expressed in terms of the indicator and , respectively, in the undersampled regime (spacebased and groundbased are on the left and right, respectively).
See, e.g., Reffert (2009), Fig. 1a in Høg (2017) for an overview spanning more than 2000 years of astrometry, and Benedict et al. (2016) for a summary of the contributions from the Hubble Space Telescope (HST) (fine guide sensors), and, of course, the exquisite prospects from Prusti et al. (2016), with applications ranging from fundamental astrophysics (Van Altena, 2013; Cacciari et al. 2016), to cosmology (Lattanzi, 2012).
The analysis can be extended to the twodimensional case as presented in Mendez et al. (2013).
The joint estimation of photometry and astrometry is the task of estimating both from the observations (see Mendez et al. 2014).
This is the solution of , with λ_{i}(α) = F̃g_{i}α + B̃_{i}, (α) being a generic variable representing the astrometric position, g_{i}(·) is given by Eq. (3), and where arg min represents the argument that minimizes the expression. More details are presented in Lobos et al. (2015).
In particular, for the very high S/N regime and assuming , Lobos et al. (2015), Proposition ^{3}) show that this gap reaches the condition
It is important to remember that the CR bound is a function of the value of the parameter to be estimated, in this case the position x_{c} (see the Fisher information in Eq. (10)).
The derivation of this identity is presented in Appendix B.2
The derivation of this result is presented in Appendix C.2
Acknowledgments
SE acknowledges support from CONICYTPCHA/MagisterNacional/2016 – 22160840. The authors want to thank the support of the Advanced Center for Electrical and Electronic Engineering, AC3E, Basal Project FB0008, PIA ACT1405, and the Chilean Centro de Excelencia en Astrofisica y Tecnologias Afines (CATA) BASAL PFB/06 from CONICYT. The authors also acknowledge support by CONICYT/FONDECYT grants # 1151213, 1170044, and 1170854. RAM also acknowledges support from the Project IC120009 Millennium Institute of Astrophysics (MAS) of the Iniciativa Cientifica Milenio del Ministerio de Economia, Fomento y Turismo de Chile. Finally, we thank the anonymous referee for his or her careful reading of our manuscript, and for pointing out the exploration of the relevant spacebased undersampled regime that led to the discussion in Sect. 6.6.
References
 Adorf, H.M. 1996, in Astronomical Data Analysis Software and Systems V, 101, 13 [NASA ADS] [Google Scholar]
 Alard, C., & Lupton, R. H. 1998, ApJ, 503, 325 [NASA ADS] [CrossRef] [Google Scholar]
 Altmann, M., Bouquillon, S., & Taris, F. 2014, Proc. SPIE, 9149, 91490P [CrossRef] [Google Scholar]
 Auer, L., & Van Altena, W. 1978, AJ, 83, 531 [NASA ADS] [CrossRef] [Google Scholar]
 Bastian, U. 2004, GAIA Technical Note, 2004 BASNOCODE [Google Scholar]
 Bendinelli, O., Parmeggiani, G., Piccioni, A., & Zavatti, F. 1987, AJ, 94, 1095 [NASA ADS] [CrossRef] [Google Scholar]
 Benedict, G. F., McArthur, B. E., Nelan, E. P., & Harrison, T. E. 2016, PASP, 129, 012001 [Google Scholar]
 Bouquillon, S., Mendez, R., Altmann, M., et al. 2017, A&A, 606, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bradley, R. A., & Gart, J. J. 1962, Biometrika, 49, 205 [CrossRef] [Google Scholar]
 Bristow, P., Kerber, F., & Rosa, M. 2006, in The 2005 HST Calibration Workshop: Hubble After the Transition to TwoGyro Mode, 299 [Google Scholar]
 Cacciari, C., Pancino, E., & Bellazzini, M. 2016, Astron. Nachr., 337, 899 [NASA ADS] [CrossRef] [Google Scholar]
 Chromey, F. R. 2016, To Measure the Sky: An Introduction to Observational Astronomy (Cambridge: Cambridge University Press) [Google Scholar]
 ChunLin, L. 1993, IAU Symp., 156, 113 [NASA ADS] [Google Scholar]
 Cover, T. M., & Thomas, J. A. 2012, Elements of Information Theory (New York: Wiley) [Google Scholar]
 Cramér, H. 1946, Scand. Actuar. J., 1946, 85 [CrossRef] [Google Scholar]
 Echeverria, A., Silva, J. F., Mendez, R. A., & Orchard, M. 2016, A&A, 594, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fessler, J. A. 1996, IEEE Trans. Image Process., 5, 493 [NASA ADS] [CrossRef] [Google Scholar]
 Gai, M., Busonero, D., & Cancelliere, R. 2017, PASP, 129, 054502 [NASA ADS] [CrossRef] [Google Scholar]
 Gray, R. M., & Davisson, L. D. 2004, An Introduction to Statistical Signal Processing (Cambridge: Cambridge University Press) [CrossRef] [Google Scholar]
 Hoadley, B. 1971, Ann. Math. Stat., 1977 [CrossRef] [Google Scholar]
 Høg, E. 2017, ArXiv eprints [arXiv:1707.01020] [Google Scholar]
 Howell, S. B. 2006, Handbook of CCD Astronomy (Cambridge: Cambridge University Press), 5 [Google Scholar]
 Jakobsen, P., Greenfield, P., & Jedrzejewski, R. 1992, A&A, 253, 329 [NASA ADS] [Google Scholar]
 Janesick, J. R. 2001, Scientific ChargeCoupled Devices (Bellingham: SPIE Press), 83 [Google Scholar]
 Janesick, J. R. 2007, Photon Transfer (San Jose: SPIE Press) [CrossRef] [Google Scholar]
 Kay, S. M. 1993, Fundamentals of Statistical Signal Processing. Vol 1, Estimation Theory (Englewood Cliffs: PrenticeHall) [Google Scholar]
 Kendall, M., Stuart, A., Ord, J., & Arnold, S. 1999, Kendall’s Advanced Theory of Statistics. Vol. 2A (London: Hodder Arnold Publication) [Google Scholar]
 King, I. R. 1971, PASP, 83, 199 [NASA ADS] [CrossRef] [Google Scholar]
 King, I. R. 1983, PASP, 95, 163 [NASA ADS] [CrossRef] [Google Scholar]
 Lattanzi, M. 2012, Mem. Soc. Astron. It. 83, 1033 [NASA ADS] [Google Scholar]
 Lee, J.F., & Van Altena, W. 1983, AJ, 88, 1683 [NASA ADS] [CrossRef] [Google Scholar]
 Lemon, C. A., Auger, M. W., McMahon, R. G., & Koposov, S. E. 2017, MNRAS, 472, 5023 [NASA ADS] [CrossRef] [Google Scholar]
 Lindegren, L. 1978, IAU Colloq., 48, 197 [NASA ADS] [Google Scholar]
 Lindegren, L. 2008, Gaia DPAC Public Document GAIAC3TNLULL078 [Google Scholar]
 Lindegren, L. 2010, ISSI Scientific Reports Series, 9, 279 [NASA ADS] [Google Scholar]
 Lobos, R. A., Silva, J. F., Mendez, R. A., & Orchard, M. 2015, PASP, 127, 1166 [NASA ADS] [CrossRef] [Google Scholar]
 McLean, I. S. 2008, Electronic Imaging in Astronomy: Detectors and Instrumentation (New York: Springer Science & Business Media) [Google Scholar]
 Méndez, R. A., Costa, E., Pedreros, M. H., et al. 2010, PASP, 122, 853 [NASA ADS] [CrossRef] [Google Scholar]
 Mendez, R. A., Silva, J. F., & Lobos, R. 2013, PASP, 125, 580 [NASA ADS] [CrossRef] [Google Scholar]
 Mendez, R. A., Silva, J. F., Orostica, R., & Lobos, R. 2014, PASP, 126, 798 [NASA ADS] [Google Scholar]
 Michalik, D., & Lindegren, L. 2016, A&A, 586, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Michalik, D., Lindegren, L., Hobbs, D., & Butkevich, A. G. 2015, A&A, 583, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Prusti, T., De Bruijne, J., Brown, A. G., et al. 2016, A&A, 595, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rao, C. R. 1945, Bull. Calcutta Math. Soc., 37, 81 [Google Scholar]
 Reffert, S. 2009, New Astron. Rev., 53, 329 [NASA ADS] [CrossRef] [Google Scholar]
 So, H. C., Chan, Y. T., Ho, K., & Chen, Y. 2013, IEEE Signal Process. Mag., 30, 162 [NASA ADS] [CrossRef] [Google Scholar]
 Stetson, P. B. 1987, PASP, 99, 191 [NASA ADS] [CrossRef] [Google Scholar]
 Stone, R. C. 1989, AJ, 97, 1227 [NASA ADS] [CrossRef] [Google Scholar]
 Vakili, M., & Hogg, D. W. 2016 ArXiv eprints [arXiv:1610.05873] [Google Scholar]
 Van Altena, W. F. 2013, Astrometry for Astrophysics: Methods, Models, and Applications (Cambridge: Cambridge University Press) [Google Scholar]
 Van Altena, W. F., & Auer, L. 1975, in Image Processing Techniques in Astronomy (Berlin: Springer), 411 [NASA ADS] [CrossRef] [Google Scholar]
 Van Trees, H. L. 2004, Detection, Estimation, and Modulation Theory, Part I: Detection, Estimation, and Linear Modulation Theory (New York: Wiley) [Google Scholar]
 Zaccheo, T., Gonsalves, R., Ebstein, S., & Nisenson, P. 1995, ApJ, 439, L43 [NASA ADS] [CrossRef] [Google Scholar]
 Zhang, J., Hao, Y. C., Wang, L., & Long, Y. 2016, PASP, 128, 035003 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A
Proof of Theorem 1
We begin presenting the expressions for to complete the statement of the result.(A.1) (A.2)
Proof of Theorem ^{1}: using the chain rule in the cost function J(α, I^{n}) and taking the partial derivative of both sides in Eq. (15), we have that(A.6)
Thus, we have n equations with one unknown, and it holds for any I^{n}. Defining the operators(A.7)
of dimensions 1 × 1 and 1 × n, respectively, we can express Eq. (A.6) in matrix form as(A.8)
Assuming that the matrix ∇^{20} J(τ(I^{n}), I^{n}) is non singular, we can calculate ∇_{τ}(I_{n}) from Eq. (A.8),(A.9)
Finally, using Eq. (A.9), evaluating at Ī^{n}, and then replacing in Eq. (21), we have that(A.10)
Moving into the residual term γ_{J}(n) in Eq. (21) captured by β_{J}(n), we must consider the variance of the error function Var{e(I̅^{n}, I^{n} ‒ I̅^{n})} and the covariance . For the first, we have that(A.11)
On the other hand, for the covariance, using the main assumption in Eq. (22), it is clear that(A.12)
Finally, replacing Eqs. (A.11) and (A.13) in the definition of γ_{J}(n), we have that(A.14)
For the bias expression of the result in Eq. (24), using the hypothesis in Eq. (23), we can take expectation at both sides of Eq. (18) to obtain that(A.15)
Appendix B
Proof of Theorem 2
Proof: the proof and, in particular, the derivation of , β_{WLS}(n) and ϵ_{WLS}(n) simply reduces to a straightforward application of Theorem ^{1}. For that we need to first validate the assumptions of Theorem ^{1}. If we begin with Eq. (A.9),(B.1)
and then we calculate the gradient terms on the right hand side (RHS) of Eq. (B.1) for our WLS context, it follows that(B.2) (B.3) (B.4)
Following Eq. (B.1), we need to evaluate Eqs. (B.3) and (B.4) at α = τ_{WLS}(I̅^{n}). For that, we have the following^{12}:(B.5)
Then we will use the following result.
Proposition 2Under the assumption of a Gaussian PSF, τ_{WLS}(I̅^{n}) = x_{c}.
We note that this proposition is the second assumption used in Theorem ^{1}. Using this proposition, we obtain that (B.6) (B.7)
Finally, applying Eqs. (B.6) and (B.7) in Eq. (A.9), we have that(B.8)
which offers the decomposition needed for the application of Theorem ^{1} (Eq. (22)). For the value of in Eq. (A.5), since the observations are independent and follow a Poisson distribution, we have that(B.9)
Then if we replace Eqs. (B.6), (B.7), and (B.9) in Eq. (A.5), we have that(B.10)
On the other hand, the expression for β_{WLS}(n) and ϵ_{WLS}(n) can be determined from the evaluation of Eqs. (A.2) and (A.1), respectively. Looking at them, the problem reduces to determine the key term . For that, if we use Fessler (1996), Eq. (17) we can obtain the following identity^{13},(B.11)
which concludes the result.
B.1. Proof of Proposition 2
Proof: using the function , we need to show that the minimum is reached only at α = x_{c}. From this, we have that h(α) ≥ 0 and it achieves its minimum at x_{c}. To prove uniqueness, let us assume that there is another position at which h is zero. Then(B.12)
The last identity is not possible, because if we use a Gaussian PSF there is at least one i ∈ {1, . . . , n} such that .
Proof of Eq. (B.11)
Proof: if we recall Fessler (1996), Eq. (17) and consider J_{WLS}(α, I^{n}) as the cost function, we have that(B.13)
where from the definition of J_{ML} (α, I^{n}) we have that(B.14) (B.15) (B.16) (B.17)
Concerning is just the ith component of the gradient in Eq. (B.1), then we use Eqs. (B.3) and (B.4)(B.18)
Finally, replacing Eqs. (B.14)–(B.18) in Eq. (B.13), and evaluating in α = τ_{WLS}(I̅^{n} + t(I^{n} ‒ I̅^{n})), we obtain the desired result(B.19)
Appendix C
Proof of Theorem 3
Proof: again the proof and the derivation of , β_{ML}(n) and ϵ_{ML}(n) reduce to apply Theorem ^{1}. First, we need to validate the assumption of Theorem ^{1}. Beginning with the equality in Eq. (B.1), it follows that(C.1) (C.2)
To evaluate these two expression at α = τ_{ML}(I̅^{n}) as required in Eq. (B.1), we use that^{14} (C.3)
Then we will use the following result.
Proposition 3Under the assumption of a Gaussian PSF, α = τ_{ML}(I̅_{n}) = x_{c}.
We notice again that this proposition is the second assumption used in Theorem ^{1}. From this proposition, it follows that(C.4) (C.5)
Finally, we apply Eqs. (C.4) and (C.5) in Eq. (A.9) to obtain that(C.6)
which shows that the sufficient condition in Eq. (22) of Theorem ^{1} is satisfied. To compute the value in Eq. (A.5), we have that(C.7)
since the observations are independent and follow a Poisson distribution. Then, replacing Eqs. (C.4), (C.5), and (C.7) in Eq. (A.5), we have that(C.8)
which resolves the identity in Eq. (37). Finally β_{ML}(n) and ϵ_{ML}(n) come from evaluating Eqs. (A.2) and (A.1) in this ML context. For that we only need to determine . Using Fessler (1996), Eq. (17), we can obtain the identity^{15} (C.9)
which concludes the result.
Proof of Proposition 3
Proof: let us consider the function given by(C.10)
where applying the first order condition y = λ_{i}, ∀i ∈ {1, . . . , n}. Returning to our problem in Eq. (C.3) where λ_{i} = I̅_{i}(x_{c}) and y_{i} = λ_{i}(α), it is clear, considering the Gaussian profile in PSF, that(C.12)
which concludes the result.
Proof of Eq. (C.9)
Proof: if we recall Fessler (1996, Eq. (17)) and consider J_{ML}(α, I^{n}) as the cost function, we have that(C.13)
where from the definition of J_{ML}(α, I^{n}) we have that(C.14) (C.15) (C.16) (C.17)
Concerning , this expression is the ith component of the gradient in Eq. (B.1), then we use Eqs. (C.1) and (C.2)(C.18)
Finally, replacing Eqs. (C.14)–(C.18) in Eq. (C.13), and evaluating in α = τ_{ML}(I̅^{n} + t(I^{n} ‒ I̅^{n})), we obtain the desired result(C.19)
All Tables
Performance quality of the ML estimator relative to the Cramér–Rao bound expressed in terms of the indicator from the result in Theorem ^{3}.
Performance quality of the ML and WLS estimators relative to the nominal bound expressed in terms of the indicator and , respectively, in the undersampled regime (spacebased and groundbased are on the left and right, respectively).
All Figures
Fig. 1.
Relative performance of the bias (as measured by ) stipulated by Theorem ^{1} for the WLS estimator (left side, Eq. (29)) and the ML estimator (right side, Eq. (35)). Results are reported for different values of the source flux F̃ ∈ {1080, 3224, 20 004, 60, 160}, all in e^{−} (top to bottom symbols respectively), as a function of the detector pixel size. The 0% level corresponds to having achieved no bias. 

In the text 
Fig. 2.
Range of the square root of the variance performance (in miliarcsecond = mas) for the WLS method in astrometry using uniform weights (equivalent to the LS method) predicted by Theorems ^{1} and ^{2}, Eq. (30). Results are reported for different representative values of F̃ and across different pixel sizes (topleft to bottomright): F̃ ∈ {1080; 3224; 20 004; 60 160} e^{−}. 

In the text 
Fig. 3.
Worstcase discrepancies in Eq. (40) for the WLS estimator using the weights set indexed by the positions . Results are reported for two S/N scenarios, namely F̃ = 20 004 e^{−} (left) and F̃ = 60 160 e^{−} (right) (right), and across different pixel sizes. 

In the text 
Fig. 4.
Performance discrepancies (measuring the nonoptimality) of the WLS estimator using the center position as a prior for the weight selection with respect to the CR bound obtained for the true object positions . Results are reported for two S/N scenarios, namely F̃ = 20 004 e^{−} (left) and F̃ = 60 160 e^{−} (right), and across different pixel sizes. 

In the text 
Fig. 5.
Range of the square root of the variance performance (in miliarcsecond = mas) for the ML method in astrometry as predicted by Theorems ^{1} and ^{3}, Eq. (36). Results are reported for different representative values of F̃ and across different pixel sizes (topleft to bottomright): F̃ ∈ {1080, 3224, 20 004, 60 160} e^{‒}. 

In the text 
Fig. 6.
Performance optimality of the ML estimator (computed as ) for different positions of the target object in the array, as a function of pixel resolution. The left panel shows the case Fį = 20 004 e^{‒}, the right panel the case Fį = 60 160 e^{‒}. 

In the text 
Fig. 7.
Performance comparison between the of the adaptive WLS estimator and σ_{CR}(n), both in mas. Results are reported for different F̂ and across different pixel sizes: (topleft to bottomright) F̃ ∈ {1080, 3224, 60 160}e^{‒}. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.