Gaussian profile estimation in one dimension
Nathan Hagen*, Matthew Kupinski, and Eustace L. Dereniak
College of Optical Sciences, University of Arizona, Tucson, Arizona 85721, USA*Corresponding author: nhagen@optics.arizona.edu
Received 23 February 2007; revised 10 May 2007; accepted 11 May 2007;
posted 21 May 2007 (Doc. ID 80302); published 23 July 2007

We present several new results on the classic problem of estimating Gaussian profile parameters from a
set of noisy data, showing that an exact solution of the maximum likelihood equations exists for additive
Gaussian-distributed noise. Using the exact solution makes it possible to obtain analytic formulas for thevariances of the estimated parameters. Finally, we show that the classic formulation of the problem is
actually biased, but that the bias can be eliminated by a straightforward algorithm. © 2007 OpticalSociety of America
OCIS codes: 030.6600, 120.5240, 120.3940, 300.3700.

Preview

The problem of estimating the parameters of a
Gaussian-shaped function in a noisy data set (such as
the one shown in Fig. 1) is a common one in many
fields. For one-dimensional (1D) data sets, the mostprominent use of Gaussian profile estimation is probably in spectroscopy, where it is used to accurately
locate the wavelength of a Doppler-broadened emission or absorption line and to estimate the line
strength [1– 4]. The technique can also be found in
applications of position sensing for line-objects whichhave been blurred by the imaging process [5,6]. For
two-dimensional data sets, the two-dimensional (2D)
version of the technique is even more widely used, infields such as in astrometry, wave front sensing,
Gaussian beam characterization, bio-imaging, andalso in calibrations of volumetric sensors. Achieving
the most accurate possible estimate of the desired
parameters is of course an important goal in all ofthese applications.
In this paper we analyze the use of maximumlikelihood estimation (MLE) for a 1D Gaussian profile
corrupted by additive Gaussian-distributed noise anddiscuss the benefits and drawbacks of this approach.
The MLE equations for this problem, while not solvable analytically, can nevertheless be solved exactly
0003-6935/07/220001-10$15.00/0© 2007 Optical Society of America
5374

APPLIED OPTICS 兾 Vol. 46, No. 22 兾 1 August 2007

by a straightforward nonlinear optimization procedure. Previous efforts have all used approximations
in order to solve the MLE equations, but keeping toexact solutions allows us to use the Cramer–Rao
Bound (CRB) to obtain the parameter variances. The
approach is also sufficiently general to allow nonuniform sampling of the profile, thus allowing for bad
pixels or missing data. Although the exact form of theparameter covariance matrix is complicated, applying some reasonable approximations can reduce themto surprisingly simple analytic expressions. The various models for the background underneath theGaussian profile can have an appreciable effect here,
and we extend the result to more advanced background models.
We also show that the typical formulation of theGaussian profile estimation problem produces a biased statistical model, and that the bias error can belarge for systems which insufficiently sample the profile. However, because the expression for the bias
term is simple, we show a straightforward techniquefor removing the bias error. In the conclusion, we
combine our results for the variances and the bias todiscuss the problem of estimating the Gaussian profile position. This has been treated in numerical simulations of previous publications, but our results here
do not agree with previously published research.
This paper discusses the Gaussian profile estimation problem only for the 1D case. The extension to
2D is given in a subsequent publication.

Here Ᏼ is an operator which maps a continuous function onto a discrete one. Typically, Ᏼ is chosen such
that

冕

⬁

g៮m ⫽ ␦xf(x)␦(x ⫺ xm)dx⫺⬁⫺(xm ⫺x̄)2兾2w2⫽ ␦xAeFig. 1. Example data for a 1D Gaussian profile, generated fromthe signal (A ⫽ 100.0, x៮ ⫽ 23.5, w ⫽ 8.0, ␴ ⫽ 10.0, ␦x ⫽ 1, andM ⫽ 101). The solid line represents the object f共x兲 (true underlying
Gaussian profile); the filled circles represent the measured data
values, g.

1. Introduction

One of the great advantages of MLE, as opposed toleast-squares minimization, is that for a large number of measurements it provides the minimumvariance unbiased estimate of the parameter set–it is“efficient”. If there is any prior information on thevarious parameters (such as the positivity of theGaussian peak’s amplitude and width) then it iswell-known that Bayesian estimation techniques canprovide better estimates. The disadvantage with
Bayesian methods is that they result in more complex
expressions both for the parameter estimators andfor their variances, such that numerical rather thananalytic solution techniques must be used.
An additional advantage of MLE is that the
Cramer–Rao bound provides (again, asymptoticallyfor the case of a large number M of measurements)
not just a lower bound on the variance but the variance itself. Thus, it is unnecessary to locate the parameter variances by some numerical or Monte Carlo
method; a straightforward technique is available
which can provide analytical expressions for the variances of the fitted parameters.
We start with a zero-background model for the estimation problem (see Appendix A for the nonzero
background model). The deterministic function wewish to measure as accurately as possible is what wewill refer to as the “object” f共x兲, the Gaussian profile
f(x) ⫽ Ae⫺(x⫺x̄) 兾2w ,
2

2

(1)

where A is the peak amplitude, x៮ is the position of thepeak (the center), and w the Gaussian width (w ⫽FWHM兾冑8 ln 2 ⬇ FWHM兾2.35). The quantity wewish to measure is a continuous function, but thedetection process samples it to produce a discrete
data vector ḡ:ḡ ⫽ Ᏼf(x).

,(2)

where ␦共·兲 is the Dirac ␦-function, m is the pixel index(1 ⱕ m ⱕ M), ␦x is the pixel width, and xm gives theabscissa position for the center of pixel m. (In thispaper, we assume that the pixel width coincides withthe spacing between samples—i.e., 100% fill factor.)
While we are not compelled to make this choice for Ᏼ,
it produces a model which is simple both analyticallyand conceptually. The disadvantage of the choice isthat it oversimplifies the measurement process andintroduces a bias into the result. In Section 6 below
we present a more realistic operator for Ᏼ and discuss
how to modify the problem to incorporate the necessary changes.
When the object f共x兲 is detected, such as by an array
detector, the measurement result g can be modeled
as the sum of a noiseless discrete data vector, ḡ, anda random perturbation ng ⫽ ḡ ⫹ n ⫽ Ᏼf(x) ⫹ n,(3)

where the noise vector n is a Gaussian-distributed
random variable, such thatpr(gm) ⫽1

冑2␲␴m2e⫺(gm⫺g៮m) 兾2␴m ,
2

2

(4)

and ␴m2 is the variance of the noise at pixel m. Thisapproach and nomenclature, adopted from the field ofimage science, [7] may appear as excessive mathematical machinery for the problem but greatly simplifies the bias discussion (Section 6).
There are several reasons why we choose to base ourderivation on Gaussian statistics rather than Poisson
statistics. While one reason is that the Gaussian distribution closely approximates the Poisson case for signals with more than a few digital counts, the primaryreason is that the Gaussian model is more flexible,
allowing the ability to work with postprocessed data.
For example, applying pixel gain nonuniformity correction (“flat-fielding”) to the data results in pixel values which are real-valued and not integer-valued,
violating the assumption required for proper Poisson
statistics. Applying background subtraction invariably
produces a large number of negative-valued pixels inthe data, which is of course problematic in Poisson
statistics. In a naïve approach it is common to truncate
the negative values by setting them to zero and thengoing on to perform the estimation. This obviously introduces a positive bias into the result. Another approach to the background problem is to incorporateparameters in the fit which model the background signal as well, making direct subtraction unnecessary.
1 August 2007 兾 Vol. 46, No. 22 兾 APPLIED OPTICS

5375

This may or may not be a good choice, because theintroduction of extra parameters can significantly
worsen the accuracy of estimation for the remaining
parameters. Also problematic is the fact that thisapproach reduces the stability of the optimization
routine. It is usually best to avoid these issues byperforming an independent measurement in which
the background estimation can be done accurately.
The ML estimation process steps are:
1. Construct the likelihood function L ⫽ pr共g|␪兲,
the probability measuring the data given a set ofobject parameters ␪ ⫽ 共A, x៮ , w兲.
2. Calculate the gradient of the log-likelihood
function, often referred to as the “score”: ⵱ᐉ ⫽
⭸
ln L.
⭸␪i
3. Solve for the parameter set which maximizes
the score.
The probability of measuring a single data value gm
is given by pr共gm|␪兲. Assuming independence among
all measured data, the likelihood of the entire data
set is the product of the individual measurements,M

L ⫽ pr共gⱍ␪兲 ⫽ 兿

m⫽1

冉 冊 冋
1

1兾2

2␲␴m2

exp ⫺

冉

(gm ⫺ g៮m)2
2␴m2

册

M

(gm ⫺ g៮m)2

m⫽1

2␴m2

⫽ 关兿共2␲␴m2兲⫺1兾2兴exp ⫺ 兺

冊

.

(5)
Taking the logarithm to obtain the log-likelihood,冋ᐉ ⫽ ln L ⫽ ⫺册M (gm ⫺ g៮m)2
1
ln共2␲␴m2兲 ⫺ 兺
.
兺
2
m⫽1
2␴m2

Everything we want to estimate uses partial derivatives of ᐉ with respect to the various parameters, sothat we are free to drop any terms independent of theparameters, giving
M

1

m⫽1

2␴m2

ᐉ⬃⫺ 兺
⬃兺

1
␴m

2

(gm ⫺ g៮m)2(gmg៮m ⫺

1 2
g៮ ).
2 m

(6)

To obtain the ML estimators Â, x៮ˆ , and ŵ for theparameters A, x៮ , and w, we can set the appropriate components of the gradient to zero and solve. We
quickly find, however, that we obtain a set of three
coupled nonlinear equations, foiling access to an analytic solution. Previous approaches to this problem
APPLIED OPTICS 兾 Vol. 46, No. 22 兾 1 August 2007

␪(k⫹1) ⫽ ␪(k) ⫺ (H(k))⫺1⵱ᐉ(k).
If it is necessary to avoid calculating the matrix
inverse directly, then the update vector may also be
obtained by an iterative solution technique such as
QR-decomposition followed by back-substitution. Thisis a numerically stable method for obtaining a vector⌬␪ solving the matrix equation H共k兲⌬␪ ⫽ ⫺ⵜᐉ共k兲.
Calculating the elements of the gradient gives⭸ᐉ

2. Nonlinear Optimization by Newton’s Method

5376

1. Obtain initial guesses for the parameters ␪共0兲
⫽ 共A共0兲, x៮共0兲, w共0兲兲
2. Calculate the gradient vector ⵱ᐉ共k兲
3. Calculate the Hessian matrix H共k兲4. Check the algorithm stopping conditions
5. Update the parameter estimate with the equation⭸ᐉ ␦ A⵱ᐉ⫽E␰␥⫽关 兴 ⭸x៮兺,w⭸ᐉ ␦ A关⵱ᐉ兴 ⫽ ⭸w ⫽
兺E ␰ ␥w关⵱ᐉ兴1 ⫽ ⭸A ⫽ ␦x 兺 Em␥x

The elements of the gradient then have the form
⭸g៮m
⭸ᐉ
1
⫽ 兺 2 (gm ⫺ g៮m)
.
⭸␪i
⭸␪i
␴m

resorted to various approximation techniques to overcome this difficulty, but using a simple nonlinear
optimization routine allows access to the exact solution, an advantage which afterwards makes using the
Cramer–Rao bound credible. A disadvantage of theuse of a nonlinear optimization routine is that it introduces the real possibility of finding false solutions(such as local rather than global maxima) or failing tofind a solution at all when the algorithm diverges.
Luckily, the solution space for maximum likelihood
problems, the log-likelihood function, is generally
rather well-behaved when the number of parameters
is small and the number of data large, such that
pathological behavior in the optimization algorithm
is rare.
The most obvious choice for a nonlinear optimization method is Newton’s method, both because it ispossible to calculate the gradient vector and Hessian
matrix for this problem and also for the reason thatthe Hessian matrix is itself intimately related to alater step in the analysis: calculating the Fisher information matrix. An additional advantage of selecting Newton’s method is that it is very fast. A
drawback is that it requires a reasonably close starting guess.
Locating a solution via Newton’s method involves
an iterative approach with the following procedure[8] (where a superscript indicates the iteration index):

x

2

2

x

3

(7)

m m x

2

3

m m

x

where we define ␰m ⫽ 共xm ⫺ x៮ 兲, Em ⫽ e⫺␰m 兾2w , andgm⫺g៮m␥m ⫽ ␴ 2 . All of the sums are taken over the pixelmindex, 1 ⱕ m ⱕ M. The Hessian matrix for this problem has the form
2

2

⭸2ᐉ

an equality: K ⫽ F⫺1. This is a convention we willcontinue to use in the expressions below, remembering that it holds only asymptotically.
The first step toward obtaining the parameter variances is to calculate the Fisher information matrix,
which is conveniently just the element-by-element
expectation of the Hessian matrix [9], with the expectation taken over the data: ⭸A ⭸ᐉ ⭸ᐉ ⭸ᐉ H⬅⵱ ᐉ⫽.
⭸x៮ ⭸A⭸x៮ ⭸w⭸x៮ ⭸ᐉ ⭸ᐉ ⭸ᐉ  ⭸w⭸A ⭸w⭸x៮ ⭸w 
2

⭸2ᐉ
⭸2ᐉ
⭸A⭸x៮ ⭸A⭸w
2

2

2
␪

2

2

2

2

2

2

Since H is a symmetric matrix, it has only six independent elements, and so it is only necessary to calculate the upper triangular portion of the matrix.
Performing the partial derivative calculations gives
for the elements of the matrix
H11 ⫽ ⫺␦x2兺
H12 ⫽
H13 ⫽
H22 ⫽
H23 ⫽
H33 ⫽

␦x
w2
␦x

1
␴m2

Em2

兺 Em␰m⌫m
兺 Em␰m2⌫m

w3
␦xA

w4
␦xA

兺 Em关␰m2⌫m ⫺ w2␥m兴兺 Em关␰m ⌫m ⫺ 2␰mw ␥m兴
w
3

5

␦xA

2







,

(8)

兺 Em关␰m4⌫m ⫺ 3␰m2w2␥m兴
w6

冕 冉 冊

Fij ⫽ ⫺Ᏹ兵Hij其 ⫽ ⫺ ᐉ 共␪兲

(9)

Ᏹ兵·其 indicates taking the expectation value of the argument, and ᐉ is the log-likelihood of the data. Luckily, many of the Hessian matrix elements contain thefactor 共gm ⫺ g៮m兲, the expectation of which (taken overthe data vector g) is zero. This greatly simplifies theexpressions needed to form F:
1

F11 ⫽ ␦x2兺

␴m2

␦x2A

1

F12 ⫽
F13 ⫽
F22 ⫽

兺␴

␦x2A

1

2
2 Em ␰mm

兺 ␴ 2 Em2␰m2

w3

m

␦x A

2

w4

1

兺 ␴ 2 Em ␰m
m

␦x A

2

2

F23 ⫽

Em2

w2

2

gm⫺2g៮mwhere ⌫m ⫽ ␴ 2 .
mAssuming that the algorithm converges onto theglobal solution, inserting these formulas into Newton’s method solves for the maximum likelihood estimates Â, x៮ˆ , and ŵ to desired accuracy.

⭸2 ᐉ
dMg.
⭸␪i⭸␪j

1

.

(10)

2

兺 ␴ 2 Em2␰m3
w
5

m

F33 ⫽

2







␦x2A2
w6

1

兺 ␴ 2 Em2␰m4
m

3. Calculating Parameter Variances

Once we have our ML estimates of the Gaussianfunction parameters, we would like to know how accurate they are. Previous approaches to this problem
have used numerical methods to determine the parameter variances [1,4]. We show that not only is it
possible to obtain the exact expression for the boundon the parameter covariance matrix, but with somereasonable approximations the variances can be reduced to simple analytical formulas. With these formulas it is straightforward to analyze the generalbehavior of the system.
Cramer–Rao bound theory states that the covariance matrix K of the parameters is bounded by the
inverse of the Fisher information matrix F (Barrett
and Myers p. 897):
K ⱖ F⫺1.
For the majority of problems, the number of measurements M is sufficient that the bound is approximately
reached, allowing the relation above to be replaced by

The covariance matrix K is given by the inverse of
the above matrix, easily obtained by Cramer’s rule.
First calculating the determinant D ⫽ det共F兲, we canobtain the parameter variances themselves from the
diagonal elements of K:var(Â) ⫽ K11 ⫽F22F33 ⫺ F232
,
D

(11)

var(xˆ៮ ) ⫽ K22 ⫽F11F33 ⫺ F132
,
D

(12)

F11F22 ⫺ F122.
var(ŵ) ⫽ K33 ⫽D(13)(Note that we have used the symmetry of the matrix
F to replace F21, F31 and F32 with their counterparts
F12, F13 and F23.) The numerical values of the variance
1 August 2007 兾 Vol. 46, No. 22 兾 APPLIED OPTICS

5377

3

2

4. Parameter Variances in the Uniform Noise Case

Unfortunately, in the form given above the variances
say little about the general behavior of the estimationand how the parameter variances might depend on
system specifications. For this purpose we introduce
some approximations which allow the exact variances to be reduced to a simple form. The necessary
approximations we will use are:
1. Uniform noise: ␴m ⫽ ␴ (a good approximation
for detector-noise-limited systems).
2. Uniform sampling: the spacing between sample
locations xm is constant.
3. Complete sampling: the data region xmin ⫺ xmax
is sufficient such that the Gaussian peak f共x兲 is approximately zero outside of the sampled region. Thisrequirement, together with approximation #2, allows
us to write 兺Em2␰mp ⬇ 0 when p is odd.
4. The profile is well-sampled: since each of the Fij
are sums of sampled Gaussian functions, we can approximate the sums with integrals and then replace
the integrals with known analytic results.
With these approximations, we can write1

兺 Em2 ⫽ ␦x 兺 e⫺␰ 兾w ␦x
1
⬇
␦x

m2

冕

2

⬁

e⫺␰ 兾w d␰ ⫽
2

2

冑␲w
␦x

⫺⬁

.

2

so thatvar(Â) ⬇

兺 Em ␰m ⬇ 0,

兺 Em ␰m ⬇ 2␦x ,兺 Em2␰m3 ⬇ 0,

兺 Em4␰m4 ⬇ 4␦x ,

冕

.

(16)

Ae⫺(x⫺x̄) 兾2w dx ⫽ Aw冑2␲.
2

2

⭸U
⭸␪

T

Kˆ␪

⭸U
,
⭸␪

where ⭸␪ ⫽ 共 ⭸A , ⭸x៮ , ⭸w 兲 . Therefore,

A

2

(14)

2

The corresponding covariance matrix is
5378

A2␦x冑␲

⬁var(Û) ⫽
⭸U

2

2␴2w

Some remarkable facts about these results are thatthe variances of xˆ៮ and of ŵ are equal and that thevariance of Â is independent of A. It is also satisfying
to see that the parameter variances (16) have nodependence on the position of the Gaussian, x៮ . Figure
2 shows comparisons between the empirical parameter variances as determined by a Monte Carlo simulation and the parameter variances as determinedby the approximate formulas (16).
We can also see that the variances of the parameters increase as the pixel widths shrink. This is a
result of the fact that while decreasing the pixel
widths increases the sampling rate (and thus improves the potential resolution of the estimation), it
also divides a constant number of photons into
smaller bins. These bins have uniform noise, so thatby decreasing pixel size we are adding noise fasterthan we are improving the resolution.
Often we are interested in estimating the area U
under the Gaussian curve, given byso that F can be greatly simplified to

x

,

(17)

If Â and ŵ are unbiased estimators of the respective
parameters in the Gaussian function, then the unbiased estimator of U is formed by Û ⫽ Âŵ冑2␲. As alinear function of unbiased parameters, the variance ofÛ is given by the covariances of Â and ŵ together withthe gradient of U with respect to the parameters [10],3冑␲w5

w 0 2 
 A

␦ 冑␲
F⬇
 0 2w 0  .
␴
A
3A 
0
2
4w 

2w␦x冑␲

⫺⬁

冑␲w3

2

3␴2

var共xˆ៮ 兲 ⬇ var共ŵ兲 ⬇

U⫽
2

(15)

2

x

Similarly, we find that2⫺1

 2w 0 A  2w ␴
K⬇
 0 A 0 ,
冑␲␦
 ⫺1
2w 
0
A
A 

expressions can be obtained by substituting Eq. (10)into Eqs. (11)–(13).

APPLIED OPTICS 兾 Vol. 46, No. 22 兾 1 August 2007(18)

⭸U ⭸U ⭸U T

冢 冣冉

冊冢 冣

w冑2␲ T K11 K12 K130
K21 K22 K23
var(Û) ⫽K31 K32 K33
A冑2␲

w冑2␲
0
A冑2␲

⫽2␲(w2K11 ⫹ 2AwK13 ⫹ A2K33)⫽3冑␲␴2w兾␦x,
where we have used K31 ⫽ K13.

(19)冑2␲w3Em2␰m2兺 ␴ 2 ⬇ A␦ 2 ,m

x

Em2␰m3兺 ␴ 2 ⬇ 0,m

Em2␰m4兺 ␴2⬇

3冑2␲w5mA␦x2,

so thatA 0 1  A 
0 ,
F ⬇ 冑2␲  0
w

3A 
1 0

w

(20)

 2w 0 ⫺ 2 


w
1
0
0
K⬇

,
A
冑2␲
 1
w
0
⫺
 2
2A 

(21)

w

3A

1

var(Û) ⬇ 冑2␲Aw.
Fig. 2. Comparison between the parameter variances as measured by a Monte Carlo simulation (filled circles) and the approximate parameter variance formulas (16) (curves). The model
function used for the Monte Carlo simulation is that of Fig. 1, with
one of the model parameters allowed to vary, as shown in the figure
abscissa. The data for var共Û兲 has been rescaled to fit within theplotting range. (Note that the data here was simulated with
␦-sampling, so that the bias discussed in Section 6 below is notpresent in this simulation.)(22)

Note that in the Poisson-noise approximation var共xˆ៮ 兲⫽ var共ŵ兲 and var共Â兲 is a function of parameter A.
Unlike the uniform-noise case, we find that the parameter variances are now independent of the pixel
size ␦x. This reflects the fact that, in Poisson noise, as
the pixel sizes shrink, the increased sampling rate isexactly offset by a decrease in both signal and noise ineach pixel.
6. Bias in the Standard Model

5. Parameter Variances in the
Poisson-Noise-Dominated Case

If Gaussian statistics are being used for a Poissonnoise-dominated measurement, then it helps to avoid
the approximation of uniform noise variance (i.e.,␴m2 ⬇ ␴2) and instead equate the variance to theintegrated value at a given pixel: ␴m2 ⬇ ␦xAEm. The
integrals approximating the sums (10) still have analytic solutions for this case, giving
Em2

冑2␲w

兺 ␴ 2 ⬇ A␦ 2 ,
m

Em2␰m

兺 ␴ 2 ⬇ 0,m

x

Gaussian profile estimation is a classic problem, yet
to our knowledge there has been no mention of bias in
the model. The estimation process outlined so far inthis paper minimizes the error between a discrete
data set g and a discrete object model f even though
the object is known to be continuous. Ignoring thisdifference is the source of the bias.
The approach up to this point has been to model themeasured data g as coming from a system operator Ᏼwhich performs ␦-sampling:g៮m ⫽ ␦x

冕⬁f(x)␦(x ⫺ xm)dx.

⫺⬁

A more realistic model for an array detector, however,accounts for the physical dimension of the pixel overwhich the signal is integrated:
1 August 2007 兾 Vol. 46, No. 22 兾 APPLIED OPTICS

5379

g៮m ⫽

冕冉

⬁

f(x)rect⫺⬁冊x ⫺ xm
dx,
␦x

(23)

where the rectangle function is defined asrect(x兾L) ⫽再

1 : ⱍxⱍ ⬍ L兾2.
0 : ⱍxⱍ ⬎ L兾2(24)

The discussion of polynomial functions suggests
using a Taylor expansion of the object f共x兲 across thedimensions of each pixel, providing an explicit formula for the bias component. The Taylor series expansion of f共x兲 at location xm (the midpoint of pixel m)is given byf共x兲 ⫽ 关f共x兲兴 x⫽xm ⫹ 共x ⫺ xm兲关f⬘共 x兲兴x⫽xm

Ç
model

For pixel responses modeled as rect-functions (asshown in Fig. 3), the pixel value coincides with
␦-sampling only if the continuous function being integrated is linear across the dimensions of each pixel.
If the object contains quadratic or higher terms, rectsampling and ␦-sampling no longer coincide. The difference is what we refer to here loosely as the bias,
b共x兲. Ideally, we would like to know the parameter
biases directly, and perform a correction by␪ˆ unbiased ⫽ ␪ˆ biased ⫺ ␪ˆ correction.
Unfortunately, the parameter correction terms appear to be unavailable in analytic form. Instead, we
can provide a correction directly to fˆ共x兲 in postprocessing,fˆunbiased(x) ⫽ fˆbiased(x) ⫺ b(x).

2Fig. 3. Illustration of constant, linear, and quadratic functions
across an idealized pixel. A ␦-sampling measurement model selects
the function value at the midpoint of the pixel whereas a rectmodel selects the mean of the function across the pixel as theresult. (The dotted line represents the mean of the parabolic function.)
APPLIED OPTICS 兾 Vol. 46, No. 22 兾 1 August 20072

关f⬘共 x兲兴x⫽xm ⫹ · · ·(26)

Çbias error term, b共x兲While it is possible to calculate the Taylor series
terms to any desired order, the dominant error termis the leading (parabolic) term of the series. If weapproximate the bias component using only this parabolic term, we need only calculate its mean value
across the pixel dimensions to give b共x兲. Using thesecond derivative of a Gaussian function,f⬘(x) ⫽冋A (x ⫺ x៮)2
w2

w2

册

⫺ 1 e⫺(x⫺x̄) 兾2w ,
2

2

(27)

the mean of the parabolic term at pixel m gives thebias as(25)

The choice of using rect-functions to model the
pixel responses is not compulsory. Gaussian functions are also a common choice, and will produce a
different estimate of the bias. Without knowing the
precise form of the pixel response function, however,the bias can only be approximated. Rect-functions areused here to illustrate the principle.

5380

⫹

共x ⫺ xm兲

1b(xm) ⬇ f⬘(xm)
2␦x
⫽

冋

冕

xm⫹

1
␦x
2

1
xm⫺ ␦x2
2

A␦x2 (xm ⫺ x៮)
24w2

w2

(x ⫺ xm)2dx

册⫺ 1 e⫺(xm⫺x̄) 兾2w .
2

2

(28)

The structure of this result appears reasonable: thebias increases with A and inversely with w, which
both have the effect of narrowing the Gaussian. The
bias also increases quadratically with the pixel width␦x. Figure 4(b) gives an illustration of the bias curvefor the Gaussian object of Fig. 1 but now sampled byonly 13 pixels, so that ␦x ⬇ w. For the example shown,
the bias in Â is more than 3% and the bias in ŵ is
4.7%, whereas the x៮ˆ estimator is unbiased in thequadratic approximation. Comparing this to thestandard deviations of Â and ŵ shows that the bias
term is important when the sampling becomes thissmall (␦x ⱗ w).
In order to use bias-corrected parameter estimates,
we can first perform a biased estimation step. If weassume that the estimate produces a good approximation to f共x兲 then we can calculate b共x兲 and performthe estimation once more using “corrected” data,
gcorr ⫽ g ⫺ b, where b is the ␦-sampled bias. The
usefulness of performing bias-correction is probably
limited, however. Most data used for Gaussianprofile estimation is only approximately Gaussian, sothat in most cases the modeling error is more likely to
be larger than the bias discussed here.

Fig. 5. Comparison of var共x៮ˆ 兲 determined from Eq. (7) with aMonte Carlo simulation of a Poisson-noise-dominated measurement. In descending order, from top to bottom, the lines and data
points correspond to U ⫽ 20 000, 10 000, 5000, and 2000 (inunits of mean number of photoelectrons).

Fig. 4. (a) Illustration of the three functions involved in the bias
discussion: (1) the underlying continuous function f共x兲 (grey fill),
rescaled by 1兾␦x to match the dimensions of g; (2) the piecewise
continuous pixel model (stairstep), and (3) the discrete ␦-sampled
model (spikes). (b) The bias in f̂共x兲 resulting from using data measured from an array detector while estimating with a ␦-sampled
model. The curve is calculated from Eq. (28); the filled circles
represent data obtained by measuring the bias directly in a simulation. The model f共x兲 used for both (a) and (b) is that of Fig. 1 but
measured with 13 pixels rather than 101.

7. Conclusions

Two previous publications [1,4] which use MLE for
1D Gaussian profile estimation state an importantresult that there is an optimal choice for magnification of the Gaussian profile prior to sampling by agiven detector array. If the optical system is modified
such that the Gaussian width w is approximately
1 ⬃ 2 times the pixel spacing, then Refs. 2, 5 bothconclude that the variance of the position estimate xˆ៮
will be minimized.
For the above situation, where the detector is unmodified while the optics are allowed to change magnification, A and w are not free parameters but are
constrained to enforce a constant energy U via therelation U ⫽ Aw冑2␲. Plugging this into the equationfor var共x៮ˆ 兲 in Eq. (16) for uniform-noise and Eq. (21) for
Poisson noise gives

Uniform:var(xˆ៮ ) ⫽Poisson:var(xˆ៮ ) ⫽

冑2U␴2
␲A3␦x
U
2␲A2

,

.

Combining these results with the model of the bias,
we find that we can first decrease the pixel size ␦x to
render the bias negligible, after which further decreasing the pixel size increases the variance in theuniform-noise case, or has no effect on the variance in
the Poisson-noise case. (Figure 5 shows the simulation results of varying ␦x and measuring var共Û兲 under
the case of Poisson noise.) All other things being
equal, it is usually desirable to allow the pixel widths
to be as large as possible, which implies that w
ⱗ ␦x ⱗ 2w is a reasonable choice, even for the Poissonnoise case. For this pixel size, the bias is kept small,but not so small as to be wasteful. Thus, our result isin qualitative agreement with that of Ireland [1] andWinick [4], but is in disagreement about the asymptotic behavior of the estimator as the pixel width isdecreased in the presence of pure Poisson noise.
A reason why we find it implausible, for the case of
Poisson noise, that there should be an increase inestimator variance with finer sampling is the following argument. If we compare two systems, one ofwhich samples the Gaussian profile at twice the rateof the other (i.e., contains twice as many pixels, eachof which are half the width of the other system), then
we know that we can always bin the data from thefinely-sampled system, by adding pairs of pixels together, to replicate data from the coarsely-sampled
system. If the data noise is purely driven by Poisson
statistics, then the coarsely-sampled data and thebinned finely-sampled data will have identical statistics. A properly-designed estimation algorithm
should be able to take advantage of this, and thus,
1 August 2007 兾 Vol. 46, No. 22 兾 APPLIED OPTICS

5381

there should not be a penalty for increasing the sampling rate under Poisson noise.
The authors have posted a set of programs written in the IDL [11] language located at the url
http://www.ittvis.com/codebank/index.asp (the fileis gauss_mle1.pro under the heading Statistics)
which perform the estimation algorithm discussed
in this paper. Interested readers are encouraged todownload, distribute, and use the code freely.
Appendix A.
Backgroundwhere, again, ␰m ⫽ 共xm ⫺ x៮ 兲, Em ⫽ e⫺␰m兾2w , and ␥mgm⫺g៮m⫽ ␴ 2 . The Hessian matrix becomes2m

Extension to a Model with Nonzero

If the Gaussian profile is placed on a quadratic background, then we can model the measured profile(shown in Fig. 6) asf(x) ⫽ ␦x[b0 ⫹ b1x ⫹ b2x2 ⫹ Ae⫺(x⫺x̄) 兾2w ].
2

2

The gradient of the log likelihood equation becomes⭸ᐉ关⵱ᐉ兴1 ⫽ ⭸b ⫽ ␦x兺 ␥m,

H⫽







⭸2ᐉ
⭸b02

⭸2ᐉ
⭸2ᐉ
⭸2ᐉ
⭸2ᐉ
⭸2ᐉ
⭸b1⭸b0 ⭸b2⭸b0 ⭸A⭸b0 ⭸x៮ ⭸b0 ⭸w⭸b0
⭸2ᐉ

·

⭸b12

⭸2ᐉ
⭸2ᐉ
⭸2ᐉ
⭸2ᐉ
⭸b2⭸b1 ⭸A⭸b1 ⭸x៮ ⭸b1 ⭸w⭸b1
⭸2ᐉ

⭸2ᐉ
⭸2ᐉ
⭸2ᐉ
⭸A⭸b2 ⭸x៮ ⭸b2 ⭸w⭸b2

·

·

·

·

·

·

·

·

·

·

·

·

·

⭸b

2
2

⭸2ᐉ
⭸A

2

⭸2ᐉ
⭸x៮ ⭸A

⭸2ᐉ
⭸w⭸A

⭸2ᐉ

⭸2ᐉ
⭸w⭸x៮

⭸x៮

2

·

with elements0⭸ᐉ关⵱ᐉ兴2 ⫽ ⭸b ⫽ ␦x兺 xm␥m,1⭸ᐉ关⵱ᐉ兴3 ⫽ ⭸b ⫽ ␦x兺 xm ␥m,22⭸ᐉ关⵱ᐉ兴4 ⫽ ⭸A ⫽ ␦x兺 Em␥m,⭸ᐉ ␦xA关⵱ᐉ兴5 ⫽ ⭸x៮ ⫽ 2 兺 Em␰m␥m,w⭸ᐉ关⵱ᐉ兴6 ⫽ ⭸w ⫽

␦xA
w3

H11 ⫽ ⫺␦x2兺 Sm,
H12 ⫽ ⫺␦x2兺 Smxm,
H13 ⫽ ⫺␦x2兺 Smxm2,
H14 ⫽ ⫺␦x2兺 SmEm,
H15 ⫽ ⫺

H16 ⫽ ⫺

兺 Em␰m2␥m,

␦x2A
w2
␦x2A
w3

兺 SmEm␰m,
兺 SmEm␰m2,

H22 ⫽ ⫺␦x2兺 Smxm2,
H23 ⫽ ⫺␦x2兺 Smxm3,
H24 ⫽ ⫺␦x2兺 SmEmxm,
H25 ⫽ ⫺

H26 ⫽ ⫺

␦x2A
w2
␦x2A
w3

兺 SmEmxm␰m,兺 SmEmxm␰m2,

H33 ⫽ ⫺␦x2兺 Smxm4,
Fig. 6. Example data for a 1D Gaussian profile with a quadratic
background, generated from the signal (b0 ⫽ 25.0, b1 ⫽ 0.15, b2⫽ ⫺0.0008, A ⫽ 100.0, x៮ ⫽ 23.5, w ⫽ 4.0, ␴ ⫽ 5.0, ␦x ⫽ 1, andM ⫽ 101). The solid line represents the object f共x兲 (true underlying
Gaussian profile), the filled circles represent the measured data
values, g.
5382

APPLIED OPTICS 兾 Vol. 46, No. 22 兾 1 August 2007

H34 ⫽ ⫺␦x2兺 SmEmxm2,
H35 ⫽ ⫺

␦x2A
w2

兺 SmEmxm2␰m,

⭸2ᐉ
⭸w2







H36 ⫽ ⫺␦x2A
w3

兺 SmEmxm ␰m ,2

2

H44 ⫽ ⫺␦x2兺 SmEm2,
H45 ⫽

H46 ⫽

H55 ⫽

H56 ⫽

H66 ⫽

F66 ⫽ ⫺Ᏹ兵H66其 ⫽␦x兺 Em␰m⌫m,
w
2

␦x
w3

兺 Em␰m2⌫m,

␦xA
w4
␦xA
w5

兺 Em关␰m2⌫m ⫺ w2␥m兴,兺 Em关␰m3⌫m ⫺ 2␰mw2␥m兴,

␦xA

兺 Em关␰m4⌫m ⫺ 3␰m2w2␥m兴,
w6

1gm⫺2g៮mwhere Sm ⫽ ␴ 2 and ⌫m ⫽ ␴ 2 .
mm
Since the first three rows of the Hessian matrix
elements have no dependence on the measurements
g, the values of the Fisher information matrix F and
H are the same for those elements. (The “:” notation
below indicates all elements in the row.)
F1,: ⫽ ⫺H1,:
F2,: ⫽ ⫺H2,:
F3,: ⫽ ⫺H3,:
F44 ⫽ ⫺H44

F45 ⫽ ⫺Ᏹ兵H45其 ⫽F46 ⫽ ⫺Ᏹ兵H46其 ⫽F55 ⫽ ⫺Ᏹ兵H55其 ⫽

␦x2A
w2
␦x2A
w3

兺 SmEm2␰m,
兺 SmEm2␰m2,

␦x2A2
w4

F56 ⫽ ⫺Ᏹ兵H56其 ⫽兺 SmEm2␰m2,

␦x2A2
w5
␦x2A2
w6

兺 SmEm2␰m3,
兺 SmEm2␰m4.

From here it is possible to obtain K ⫽ F⫺1 analytically, and to use approximations to simplify the resulting expressions, but the product of the effort is athicket of lengthy equations. The parameter variances for the zero background case were shown in(16) to have no dependence on the position of theGaussian, x៮ . While this remains true when we add a
constant background term (b0) to the model, it is nottrue if the background is allowed to be a linear or
quadratic function of the abscissa. The expressions
involve a frustrating dependence on x៮ , xmin, and xmax
that greatly complicate their behavior. Thus, in thenonzero background case, it is probably best to obtain
K from F numerically and select the resulting diagonal elements as the parameter variances.
This work was supported in part by DOD contract
DAAE07-02-C-L011. The authors thank Peter Jansson and William Dallas of the University of Arizona
for helpful discussions.
References
1. J. Ireland, “Precision limits to emission-line profile measuring experiments,” Astrophys. J. 620, 1132–1139 (2005).
2. D. A. Landman, R. Roussel-Dupré, and G. Tanigawa, “On thestatistical uncertainties associated with line profile fitting,”
Astrophys. J. 261, 732–735 (1982).
3. D. D. Lenz and T. R. Ayres, “Errors associated with fitting
Gaussian profiles to noisy emission-line spectra,” Publ. Astron.
Soc. Pac. 104, 1104 –1106 (1992).
4. K. A. Winick, “Cramér-Rao lower bounds on the performance of
charge-coupled-device optical position estimators,” J. Opt. Soc.
Am. A 3, 1809 –1815 (1986).
5. N. Bobroff, “Position measurement with a resolution andnoise-limited instrument,” Rev. Sci. Instrum. 57, 1152–1157(1986).
6. B. E. A. Saleh, “Estimation of the location of an optical object
with photodetectors limited by quantum noise,” Appl. Opt. 13,
1824 –1827 (1974).
7. H. H. Barrett and K. Myers, Foundations of Image Science(Wiley, 2004).
8. J. Dennis, J. E. and R. B. Schnabel, Numerical Methods for
Unconstrained Optimization and Nonlinear Equations (SIAM),
1996).
9. H. L. Van Trees, Detection, Estimation and Modulation Theory, Part 1: Detection, Estimation, and Linear Modulation Theory (Wiley, 1968).
10. M. G. Kendall and A. Stuart, The Advanced Theory of Statistics, 3rd ed. (Charles Griffin, 1969), Vol. 2.
11. The Interactive Data Language (IDL) software is developedby ITT Visual Information Systems, http://www.ittvis.com/
index.asp.

1 August 2007 兾 Vol. 46, No. 22 兾 APPLIED OPTICS

5383