STOCHASTIC PROCESSES,

DETECTION AND ESTIMATION

6.432 Course Notes

Alan S. Willsky, Gregory W. Wornell, and Jeffrey H. Shapiro

Department of Electrical Engineering and Computer Science
Massachusetts Institute of Technology
Cambridge, MA 02139

Fall 2003
 3

Estimation Theory

3.1 INTRODUCTION

This chapter of the notes provides a fairly self-contained introduction to the fun-
damental concepts and results in estimation theory. The prototype problem we
will consider is that of estimating the value of a vector x based on observations of
a related vector y. As an example, x might be a vector of the position and velocity
of an aircraft, and y might be a vector of radar return measurements from several
sensors.
As in our treatment of hypothesis testing and detection theory, there are two
fundamentally rather different approaches to these kinds of estimation problems.
In the first, we view the quantity to be estimated as a random vector x. In this case,
the conditional density py|x (y|x) fully characterizes the relationship between x and
the observation y. In the second case, we view the quantity as a nonrandom but
unknown quantity x. In this case, we express the relationship between x and the
observed data y by writing x as a parameter of the density for y, i.e., py (y; x). We
emphasize that for nonrandom parameter estimation, a probability density is not
defined for x, and as a consequence, we will not need to distinguish between x and
x in this case. Note, however, that in both the random and nonrandom parameter
cases the observations y have some inherent randomness, and hence y is always
specified probabilistically.
Before we begin, it is worth commenting that the hypothesis testing prob-
lems we considered in the last chapter of the notes can, at least in principle, be
viewed as a special case of the more general estimation problem. In particular, we
can view the M-ary hypothesis testing problem as one of estimating the value of a

107
 108 Estimation Theory Chap. 3

quantity x that takes on one of M distinct values, each of which corresponds to one
of the hypotheses H0 , H1 , . . . , HM −1 . From this perspective, at least conceptually
we can view the problem of estimation of a vector x as one of making a decision
among a continuum of candidate hypotheses. However, in practice this perspec-
tive turns out to be a better way to interpret our estimation theory results than to
first derive them. As a result, we will develop estimation theory independently.

3.2 ESTIMATION OF RANDOM VECTORS: A BAYESIAN FORMULATION

A natural framework for the estimation of random vectors arises out of what is
referred to as “Bayesian estimation theory.” This Bayesian framework will be the
subject of this section. As will become apparent, there is a close connection be-
tween the Bayesian estimation problem we consider here and the Bayesian hy-
pothesis testing problem we discussed in Chapter 2.
In the Bayesian framework, we refer to the density px (x) for the vector x ∈ Rn
of quantities to be estimated as the prior density. This is because this density fully
specifies our knowledge about x prior to any observation of the measurement y.
The conditional density py|x (y|x), which fully specifies the way in which y
contains information about x, is often not specified directly but is inferred from a
measurement model.

Example 3.1
Suppose that y is a noise-corrupted measurement of some function of x, viz.,
y = h(x) + w (3.1)
where w is a random noise vector that is independent of x and has density pw (w).
Then
py|x (y|x) = pw (y − h(x)). (3.2)
Suppose in addition, h(x) = Ax and w ∼ N (0, Λ) where the matrix A and
covariance matrix Λ are arbitrary. Then
py|x (y|x) = N (y; Ax, Λ).

Note that the measurement model py|x (y|x) and prior density px (x) together
constitute a fully statistical characterization of x and y. In particular, the joint
density is given by their product, i.e.,
py,x (y, x) = py|x (y|x) px (x) (3.3)
from which we can get all other statistical information. As an example, we can get
the marginal density py (y) for the observed data via
Z +∞
py (y) = py|x (y|x) px(x) dx.
−∞
Sec. 3.2 Estimation of Random Vectors: a Bayesian Formulation 109

In turn, we can also get the posterior density for x, i.e., the density for x given that
y = y has been observed, via
py|x (y|x) px(x)
px|y (x|y) = . (3.4)
py (y)
In our treatment, we will use x̂(y) to denote our estimate of x based on ob-
serving that the measurement y = y. Note that what we are estimating is actually
an entire vector function x̂(·), not just an individual vector. In particular, for each
possible observed value y, the quantity x̂(y) represents the estimate of the corre-
sponding value of x. We call this function the “estimator.”
In the Bayesian framework, we choose the estimator to optimize a suitable
performance criterion. In particular, we begin by choosing a deterministic scalar-
valued function C(a, â) that specifies the cost of estimating an arbitrary vector a
as â. Then, we choose our estimator x̂(·) as that function which minimizes the
average cost, i.e.,
x̂(·) = arg min E [C(x, f(y))] . (3.5)
f (·)

Note that the expectation in (3.5) is over x and y jointly, and hence x̂(·) is that
function which minimizes the cost averaged over all possible (x, y) pairs.
Solving for the optimum function x̂(·) in (3.5) can, in fact, be accomplished
on a pointwise basis, i.e., for each particular value y that is observed, we find the
best possible choice (in the sense of (3.5)) for the corresponding estimate x̂(y). To
see this, using (3.3) we first rewrite our objective function in (3.5) in the form
Z +∞ Z +∞
E [C(x, f(y))] = C(x, f(y)) px,y (x, y) dx dy
−∞ −∞
Z +∞ Z +∞ 
= C(x, f(y)) px|y (x|y) dx py (y) dy. (3.6)
−∞ −∞

Then, since py (y) ≥ 0, we clearly will minimize (3.6) if we choose x̂(y) to minimize
the term in brackets for each individual value of y, i.e.,
Z +∞
x̂(y) = arg min C(x, a) px|y (x|y) dx. (3.7)
a −∞

As (3.7) indicates, the posterior density px|y (x|y) summarizes everything we

need to know about the x and y to construct the optimal Bayesian estimators for
any given cost criterion. From this perspective, we see that the posterior density
plays a role analogous to that played by the likelihood ratio in hypothesis testing
problems. As discussed earlier, computation of the posterior density is generally
accomplished via (3.4). However, since the denominator is simply a normalization
factor (independent of x), it is worth emphasizing that we can rewrite (3.7) more
directly in terms of the measurement model and prior density as
Z +∞
x̂(y) = arg min C(x, a) py|x (y|x) px (x) dx. (3.8)
a −∞
110 Estimation Theory Chap. 3

As an additional remark, we note that the result (3.7) is actually a direct ex-
tension of the corresponding M-ary Bayesian hypothesis testing result developed
in the preceding chapter of the course notes. Specifically, if x takes on one of only
M values—which, for convenience, we label H0 , H1 , . . . , HM −1 —then the pdf for
x consists of M impulses and the integral in (3.7) becomes a summation, i.e.,
M
X −1
x̂(y) = Ĥ(y) = arg min C(Hi , a) Pr [x = Hi | y = y] . (3.9)
a∈{H0 ,H1 ,...,HM −1 } i=0

Choosing a suitable cost criterion for a particular problem depends on a vari-

ety of factors. For example, the cost criterion should reflect the relative importance
of various kinds of errors in the application of interest. However, from a practi-
cal standpoint, if we choose extremely complicated cost criteria, solving for the
optimal estimator may be intractable. As a result, selecting a good cost criterion
involves a tradeoff between capturing the aspects of interest in the error behavior
and obtaining a framework that lends itself to analysis.
In the remainder of this section, we focus on some examples of practical cost
criteria. For simplicity, we will generally restrict our attention to the case of es-
timating scalar variables x from vector observations y. Keep in mind, however,
that the more general case of estimating vector variables x can be handled in a
component-wise manner.

3.2.1 Minimum Absolute-Error Estimation

One possible choice for the cost function is based on a minimum absolute-error
(MAE) criterion. The cost function of interest in this case is
C(a, â) = |a − â|. (3.10)
Substituting (3.10) into (3.7) we obtain
Z +∞
x̂MAE (y) = arg min |x − a| px|y (x|y) dx
a −∞
Z a Z +∞ 
= arg min (a − x)px|y (x|y) dx + (x − a)px|y (x|y) dx .
a −∞ a (3.11)
Differentiating the quantity inside braces in (3.11) with respect to a gives, via Leib-
nitz’ rule, the condition
Z a Z +∞ 
px|y (x|y) dx − px|y (x|y) dx = 0. (3.12)
−∞ a a=x̂MAE (y)

Rewriting (3.12) we obtain

Z x̂MAE (y) +∞
1
Z
px|y (x|y) dx = px|y (x|y) dx = . (3.13)
−∞ x̂MAE (y) 2
Sec. 3.2 Estimation of Random Vectors: a Bayesian Formulation 111

From (3.13) we see that the x̂MAE (y) is the threshold in x of the posterior density
px|y (x|y) for which half the probability is located above the threshold and, hence,
half is also below the threshold. This quantity is more generally known as the
median of a probability density. Hence, the MAE estimator for x given y = y is the
median of the posterior density.
Note that, in general, there is no explicit formula for the median of a density,
but rather it is specified implicitly as a solution to (3.13). As a result, the median is
often calculated through an iterative, numerical optimization procedure.

Example 3.2
Suppose we have the posterior density

1/(3y)
 0<x<y
px|y (x|y) = 2/(3y) y < x < 2y

0 otherwise


Then
x̂MAE (y) = (1 + ∆)y
for an appropriate choice of ∆ > 0. To solve for ∆, we use (3.13) to obtain
1 2
·y+ · y∆ = 1/2
3y 3y
from which we deduce that ∆ = 1/4.

We also note that the median of a density is not necessarily unique.

Example 3.3
Suppose
(
1/2y 0 < x < y and 2y < x < 3y
px|y (x|y) = (3.14)
0 otherwise
Then the median of (3.14) is any number between y and 2y; hence, the MAE estima-
tors for x given y = y are all of the form
x̂MAE (y) = α
where α is any constant satisfying y ≤ α ≤ 2y (assuming y ≥ 0) or 2y ≤ α ≤ y
(assuming y < 0).

3.2.2 Maximum A Posteriori Estimation

As an alternative to that considered in the previous section, consider the cost func-
tion (
1 |a − â| > 
C(a, â) = (3.15)
0 otherwise
112 Estimation Theory Chap. 3

which uniformly penalizes all estimation errors with magnitude bigger than .
This time, substituting (3.15) into (3.7) we obtain that the minimum uniform cost
(MUC) estimator satisfies
 Z a+ 
x̂MUC (y) = arg min 1 − px|y (x|y) dx
a a−
Z a+
= arg max px|y (x|y) dx. (3.16)
a a−

Note that via (3.16) we see that x̂MUC (y) corresponds to the value of a that makes
Pr [|x − x̂MUC (y)| <  | y = y] as large as possible. This means finding the interval
of length 2 where the posterior density px|y (x|y) is most concentrated.
If we carry this perspective a little further, we see that if we let  get suffi-
ciently small then the x̂MUC (y) approaches the point corresponding to the peak of
the posterior density. For this reason, this limiting case estimator is referred to as
the “maximum a posteriori” (MAP) estimator, which we denote using
x̂MAP (y) = arg max px|y (a|y) = lim x̂MUC (y). (3.17)
a →0

The peak value of a density is referred to as its mode. Hence, we see that the
MAP estimate of x based on observing y = y is the mode of the posterior den-
sity px|y (x|y). From our limiting argument, we see that the MAP estimator can be
viewed as resulting from a Bayes’ cost formulation in which all errors are, in the
appropriate sense, equally bad.
As a final remark, we note that the vector form of the MAP is a straightfor-
ward generalization of (3.17); specifically,
x̂MAP(y) = arg max px|y (a|y). (3.18)
a

We conclude our development of MAP estimation with a brief discussion of

some computational issues. For this discussion, let us restrict our attention to the
fairly typical case in which the posterior density is differentiable in x. In this case,
we first look for the MAP estimate among the stationary points of the posterior
density—the values of x for which the Jacobian is zero, i.e.,
∂
px|y (x|y) = 0. (3.19)
∂x
Recall that in this most general vector case, (3.19) is a set of equations resulting
from differentiation with respect to each component of x. A solution of (3.19)
is a local maximum of the posterior density if the corresponding Hessian matrix
satisfies
∂2
px|y (x|y) < 0 (3.20)
∂x2
where, as discussed in Appendix 1.A, the matrix inequality in (3.20) is to be in-
terpreted in the sense of negative definiteness. If there are several local maxima,
Sec. 3.2 Estimation of Random Vectors: a Bayesian Formulation 113

the relative sizes of the posterior density px|y (x|y) at each of these must be deter-
mined. In addition, if x takes on values only in some restricted range, then the
MAP estimate may be on the boundary of this set even if (3.19) is not satisfied at
such a point. Consequently, solving for the MAP estimator in general involves
finding all values of x corresponding to local maxima of px|y (x|y) as well as all
boundary points corresponding to the range of x, and taking as x̂MAP (y) the value
that maximizes px|y (x|y) over all these points.
As a final remark, it is worth pointing out that in many problems it is more
convenient to maximize other monotonic functions of the posterior density. For
example, maximizing ln px|y (x|y) with respect to x is sometimes easier than maxi-
mizing the posterior density directly. In this example, using
py|x (y|x)px (x)
px|y (x|y) = ,
py (y)
taking logarithms, and then differentiating with respect to x we obtain the MAP
equations
∂ ∂
ln py|x (y|x) + ln px (x) = 0. (3.21)
∂x ∂x
Solutions of (3.21) that also satisfy (3.20) are again the local maxima of the posterior
density px|y (x|y).

3.2.3 Bias and Variance

Let us briefly discuss some general and useful measures of performance for esti-
mators x̂(·) regardless of the cost criterion we choose. One very important quantity
is the estimate bias. Specifically, if we define the estimation error via
e(x, y) = x̂(y) − x, (3.22)
then the bias is the average value of this error, i.e.,
Z +∞ Z +∞
b = E [e(x, y)] = [x̂(y) − x] px,y (x, y) dx dy. (3.23)
−∞ −∞

The second performance measure is the error covariance

Λe = E (e(x, y) − b)(e(x, y) − b)T .
 
(3.24)
T
 
Note that using (3.24) the associated error correlation matrix E ee can be ex-
pressed in the form
E eeT = Λe + bbT ,
 
(3.25)
so that both the bias and covariance contribute to the error correlation and, in turn,
mean-square estimation error.
Since b is a deterministic vector it is, in principle, straightforward to correct
for the bias: take as the estimate x̂(y)−b. From this perspective, we see that adding
114 Estimation Theory Chap. 3

the constraint that our estimator be unbiased need not be a serious restriction.
It should be pointed out, however, that in some problems, it may be difficult to
compute b and therefore compensate for it. In such cases, there may be a tradeoff
between choosing an estimator with a small covariance or one with a small bias.

As a final remark, if error covariance is the performance metric of primary

interest for our estimator, then the natural cost criterion is in fact the least-squares
one, which we’ll now develop in detail.

3.2.4 Bayes’ Least-Squares Estimation

In this section we consider the mean-square error (MSE) cost criterion

N
X
2 T
C(a, â) = ka − âk = (a − â) (a − â) = (ai − âi )2 (3.26)
i=1

In this case, substituting (3.26) into (3.7) yields

Z +∞
x̂BLS (y) = arg min (x − a)T (x − a) px|y (x|y) dx (3.27)
a −∞

where we have used x̂BLS (·) to specifically denote the Bayes least-squares (BLS)
estimator. Since this estimator minimizes the mean-square estimation error, it is
often alternatively referred to as the minimum mean-square error (MMSE) estima-
tor and denoted using x̂MMSE (·).

Let us begin with the simpler case of scalar estimation, for which (3.27) be-
comes
Z +∞
x̂BLS (y) = arg min (x − a)2 px|y (x|y) dx. (3.28)
a −∞

As we did in the case of MAE estimation, we can perform the minimization

in (3.28) by differentiating with respect to a and setting the result to zero to find
the local extrema. Differentiating the integral in (3.28) we obtain
Z +∞  Z +∞
∂ 2 ∂
(x − a) px|y (x|y) dx = (x − a)2 px|y (x|y) dx
∂a −∞ −∞ ∂a
Z +∞
= −2 (x − a)px|y (x|y) dx. (3.29)
−∞
Sec. 3.2 Estimation of Random Vectors: a Bayesian Formulation 115

Setting (3.29) to zero at a = x̂BLS (y) we see that

Z +∞ 
(x − a)px|y (x|y) dx
−∞ a=x̂BLS (y)
Z +∞ Z +∞
= x px|y (x|y) dx − x̂BLS (y) px|y (x|y) dx
−∞ −∞
Z +∞
= E [x|y = y] − x̂BLS (y) px|y (x|y) dx
−∞
= E [x|y = y] − x̂BLS (y) = 0. (3.30)
A simple rearrangement of (3.30) then yields our final result
x̂BLS (y) = E [x|y = y] , (3.31)
i.e., that the BLS or MMSE estimate of x given y = y is the mean of the posterior
density px|y (x|y).
Note that since the quantity being minimized in (3.28) is nonnegative, and
since our derivation above concluded that there exists a single local extremum,
this extremum—and hence our estimate (3.31)—must correspond to a global min-
imum.1
The preceding derivation generalizes rather easily when x is a vector. Specifi-
cally, since the cost criterion (3.26) is a sum of individual squared estimation errors
for the components of x, the minimum is achieved by minimizing the mean-square
estimation error in each scalar component. Hence, we obtain
x̂BLS (y) = E [x|y] , (3.33)
from which we see that in the vector case as well the BLS estimate of x given y = y
is the mean of the posterior density px|y (x|y).
Let us next turn our attention to the performance characteristics of the BLS
estimator, in particular its bias and error covariance. First, we note that the BLS
estimator is always unbiased: using (3.23) we have
bBLS = E [e(x, y)] = E [x̂BLS (y) − x] = E [E [x|y]] − E [x] = 0, (3.34)
where the last equality follows from a simple application of the law of iterated
expectation.
Next, using (3.22), (3.24), and (3.34) we obtain that the associated error co-
variance is given by
h i
 T T
ΛBLS , Λe = E ee = E (x − E [x|y]) (x − E [x|y]) (3.35)

1
We can also verify this independently by taking a second derivative of the integral in (3.28),
i.e.,
+∞ +∞
∂2
Z  Z
(x − a)2 px|y (x|y) dx = 2 px|y (x|y) dx = 2 > 0, (3.32)
∂a2 −∞ −∞
which also establishes that the objective function is convex.
116 Estimation Theory Chap. 3

where we emphasize that the notation ΛBLS is used to refer to the error covariance
of the BLS estimator. Applying iterated expectation to (3.35) we see that the error
covariance can be written as
h h ii
ΛBLS = E E (x − E [x|y]) (x − E [x|y])T | y . (3.36)

However, the inner expectation in (3.36) is simply the covariance of the posterior
density, i.e., Λx|y , which in general depends on y.2 Hence, the error covariance of
the BLS estimator is simply the average of the covariance of the posterior density,
where this averaging is over all possible values of y, i.e.,
 
ΛBLS = E Λx|y (y) . (3.37)
As a final remark before we proceed to an example, note that using the iden-
tity (1.204) from Appendix 1.A of Chapter 1, we have that at its minimum value
the expected cost objective function in (3.27) can be expressed as
h i
T
E [C(x, x̂BLS (y))] = E (E [x|y] − x) (E [x|y] − x)
h n oi
= E tr (E [x|y] − x) (E [x|y] − x)T
 h i
= tr E (E [x|y] − x) (E [x|y] − x)T
= tr (ΛBLS ) . (3.38)

Example 3.4
Suppose x and w are independent random variables that are both uniformly dis-
tributed over the range [−1, 1], and let
y = sgn x + w .
Let’s determine the BLS estimate of x given y . First we construct the joint density.
Note that for x > 0, we have
(
1/2 0 < y < 2
py |x (y|x) =
0 otherwise
while for x < 0, we have
(
1/2 −2 < y < 0
py |x (y|x) = .
0 otherwise
Hence, the joint density is

1/4
 0 < x < 1 and 0 < y < 2
px,y (x, y) = py |x (y|x) px (x) = 1/4 −1 < x < 0 and −2 < y < 0

0 otherwise


2
Note that given an observed value of y, this posterior covariance Λx|y=y is in general a
function of y. To emphasize this dependence, and for future convenience, we’ll frequently use the
alternative notation Λx|y (y) for this covariance.
Sec. 3.2 Estimation of Random Vectors: a Bayesian Formulation 117

and so for y > 0 we have

(
1 0<x<1
px|y (x|y) = (3.39a)
0 otherwise

and for y < 0 we have (

1 −1 < x < 0
px|y (x|y) = . (3.39b)
0 otherwise
Thus from (3.39) we conclude that
(
1 1/2 y>0
x̂BLS (y) = E [x|y = y] = sgn y = . (3.40)
2 −1/2 y<0

Since
λx|y (y) = 1/12
is independent of y in this example, we have that the corresponding error variance
is simply  
λBLS = E λx|y (y ) = 1/12. (3.41)

Additional Properties of BLS Estimators

Let us briefly consider some additional important properties and an alternate char-
acterization of the Bayes’ least-squares estimate x̂BLS (y). Recall that we have al-
ready shown (see (3.34)) that the Bayes’ least-squares estimate is unbiased.
Next we show that Bayes’ least-squares estimates are unique in having an
important orthogonality property. Specifically, we have the following theorem.

Theorem 3.1 An estimator x̂(·) is the Bayes’ least-squares estimator, i.e., x̂(·) = x̂BLS (·),
if and only if the associated estimation error e(x, y) = x̂(y) − x is orthogonal to any
(vector-valued) function g(·) of the data, i.e.,
E [x̂(y) − x] gT (y) = 0.
 
(3.42)

In establishing this result, it will be convenient to first rewrite the condition

(3.42) as
E xgT (y) = E x̂(y)gT (y) .
   
(3.43)
and note, using the law of iterated expectation, that the left-hand side of (3.43) can
in turn be expressed in the form
E xgT (y) = E E xgT (y) | y = E E [x|y] gT (y) .
      
(3.44)
Then, to prove the “only if” statement, simply let x̂(·) = x̂BLS (·) in (3.43), and note
that in this case the right-hand expressions in both (3.44) and (3.43) are identical,
verifying (3.42).
118 Estimation Theory Chap. 3

To prove the converse, let us rewrite (3.42) using (3.43) and (3.44) as

0 = E xgT (y) − E x̂(y)gT (y)

   

= E E [x|y] gT (y) − E x̂(y)gT (y)

   

= E [E [x|y] − x̂(y)] gT (y) .

 
(3.45)
Then, since (3.45) must hold for all g(·), let us choose g(y) = E [x|y] − x̂(y) where
x̂(·) is our estimator. In this case (3.45) becomes
h i
E [E [x|y] − x̂(y)] [E [x|y] − x̂(y)]T = 0,

from which we can immediately conclude that x̂(y) = E [x|y].3

It is worth emphasizing that Theorem 3.1 ensures what we would expect
of an estimator that yields the minimum mean-square error: that since the error
e(x, y) = x̂(y) − x is uncorrelated with any function of the data we might construct,
there is no further processing that can be done on the data to further reduce the
error covariance in the estimate.
One final property of the BLS estimator is given in terms of the following
matrix inequality (again to be interpreted in the sense of positive semidefiniteness
as discussed in Appendix 1.A). Let Λe be the error covariance of any estimator
x̂(·). Then the error covariance of the BLS estimator, i.e., ΛBLS , satisfies
ΛBLS ≤ Λe (3.46)
with equality if and only if
x̂(y) − E [x̂(y) − x] = x̂BLS (y) = E [x|y] . (3.47)
In essence, this states that Bayes’ least-squares estimator is guaranteed to yield
less uncertainty in the value of x (as measured by the covariance) than any other
estimator—biased or unbiased.
Before proving this result, we point out that a useful corollary results from
the special case corresponding to choosing x̂(y) = mx . In this case the correspond-
ing error covariance is Λe = Λx , so we have that
ΛBLS ≤ Λx (3.48)
with equality if and only if
x̂BLS (y) = E [x|y] = E [x] . (3.49)
We stress that as discussed in Chapter 1, (3.49) is not equivalent to x and y being ei-
ther statistically independent or uncorrelated. While x and y being uncorrelated is
a necessary condition for (3.49) to hold, it is not sufficient one. On the other hand,
for (3.49) to hold it is sufficient but not necessary that x and y be independent.

3
 T
Here we are using a straightforward consequence of the Chebyshev inequality—that if
E zz = 0 then z = 0, or more precisely, Pr [z = 0] = 1.
Sec. 3.2 Estimation of Random Vectors: a Bayesian Formulation 119

To prove (3.46), let b denote the bias in our estimator x̂(·), let
g(y) = x̂(y) − x̂BLS (y) − b, (3.50)
and let us begin by noting that
h i
T
Λe = E (x̂(y) − x − b) (x̂(y) − x − b)
h i
= E [g(y) + (x̂BLS (y) − x)] [g(y) + (x̂BLS (y) − x)]T
h i
= E g(y)gT (y) + E (x̂BLS (y) − x) (x̂BLS (y) − x)T
 

T
+ E (x̂BLS (y) − x) gT (y) + E (x̂BLS (y) − x) gT (y) .
  
(3.51)
From Theorem 3.1 we get that the last two terms in (3.51) are zero. Using this
together with the definition of ΛBLS we get
Λe − ΛBLS = E g(y)gT (y) .
 
(3.52)
The right-hand side of (3.52) is in general positive semidefinite, which verifies
(3.46), and equal to zero if and only if g(y) = 0, which using (3.50) yields (3.47).

Weighted Least-Squares Estimators

As a final comment on Bayes’ least-squares estimation, in this section we show

that the posterior or conditional mean E [x|y] is also the optimal estimator for a
more general weighted least-squares cost criterion.
In particular, consider the weighted least-squares cost criterion
C(a, â) = (a − â)T M(a − â), (3.53)
where M is an arbitrary positive definite matrix. In this case (3.7) becomes
Z +∞
x̂WLS (y) = arg min (x − a)T M(x − a)px|y (x|y) dx. (3.54)
a −∞

The objective function in (3.54) is a function of the vector of variables com-

prising a. To find the minimum of this function we set the Jacobian of the function
to zero and look for local extrema. Applying this procedure to (3.54) we obtain
 Z +∞ 
∂ T
0= (x − a) M(x − a) px|y (x|y) dx
∂a −∞ a=x̂WLS (y)
Z +∞ 
∂  T

= (x − a) M(x − a) px|y (x|y) dx
−∞ ∂a a=x̂WLS (y)
 Z +∞ 
= 2 (x − a)T M px|y (x|y) dx . (3.55)
−∞ a=x̂WLS (y)
120 Estimation Theory Chap. 3

Minor rearrangement of (3.55) then yields

ME [x|y] = Mx̂WLS (y), (3.56)
which, since M is invertible, implies that (3.56) has a unique solution, and hence
the objective function in (3.54) has a unique local extremum. But since M is posi-
tive definite, the cost function C(·, ·) in (3.53) is non-negative for every value of a.
Hence, the unique local extremum must be a global minimum;4 hence from (3.56)
we obtain
x̂WLS (y) = E [x|y] = x̂BLS (y). (3.57)

3.2.5 Linear Least-Squares Estimation

Two important observations about the Bayes’ least-squares estimator x̂BLS (y) =
E [x|y] should be made. First, this estimator is in general a nonlinear (and often
highly nonlinear) function of the data y. Second, computing this estimator requires
that we have access to a complete statistical characterization of the relationship be-
tween x and y. In particular, we need full knowledge of px|y (x|y) or, equivalently,
py|x (y|x) and px (x).
However, there are many application scenarios when even though the least-
squares cost criterion is appropriate, the resulting Bayes’ least-squares estimator
is not practicable either because implementing the nonlinear estimator is compu-
tationally too expensive, or because a complete statistical characterization of the
relationship between x and y is not available from which to compute the estimator.
In these situations, we often must settle for a suboptimal estimator. One way
to obtain such an estimator is to add a constraint on the form of the estimator. As
an important example, in this section of the notes we’ll develop in detail estima-
tors that minimize the average Bayes’ least-squares cost (3.26), but subject to the
additional constraint that the estimator be a linear5 function of the data. Specifi-
cally, we let x̂LLS (·) denote this linear least-squares (LLS) estimator, and define it
as
x̂LLS (·) = arg min E kx − f(y)k2
 
(3.58a)
f (·)∈B

where
B = {f(·) | f(y) = Ay + d for some A and d} . (3.58b)

4
In fact, the objective function not only has a unique global minimum, but is convex as well:
Z +∞ Z +∞
∂2 T
(x − a) M(x − a) px|y (x|y) dx = M px|y (x|y) dx = M > 0.
∂a2 −∞ −∞

5
Throughout this course, we’ll use the term “linear” to refer to estimators of the form Ay + d
since this has become standard practice in estimation theory. More precise terminology would have
us refer to such estimators as “affine” and estimators of the form Ay as “linear.”
Sec. 3.2 Estimation of Random Vectors: a Bayesian Formulation 121

As we’ll see, this estimator is not only particularly efficient to implement, but
we’ll need access to only the joint second-order statistics of x and y in order to
compute it.
Although there are a variety of ways to derive the optimum estimator, we’ll
follow a powerful approach based on the abstract vector space concepts we de-
veloped in Section 1.7. With this formulation, several important perspectives and
properties of linear estimators will become apparent.
The key to exploiting vector space concepts in this problem lies in reinter-
preting (3.58) as a problem in linear approximation. For simplicity, we’ll begin by
examining the case in which we wish to estimate a scalar x based on observation
of a vector  T
y = y1 y2 · · · yM . (3.59)
In particular, let h·, ·i and k · k denote the inner-product and associated norm,
respectively, of the inner product space V = L2 (Ω) of finite mean-square random
variables with
hx, y i = E [xy ] (3.60)
and thus
kxk2 = E x 2 .
 
(3.61)
Then we can rewrite (3.58) as
x̂LLS (y) = arg min kw − xk2 (3.62a)
w ∈Y

where
Y = span(1, y1, y2 , . . . , yM )
( M
)
X
= w ∈V|w =d+ ai yi for some a1 , a2 , . . . , aM . (3.62b)
i=1

Note that by rephrasing our problem in terms of (3.62), we have abstracted

our estimation problem into one of approximation in an arbitrary inner product
space. That is, we can now develop a solution to (3.62) without having to take
into account the particular inner product that is relevant to this problem. As a
result, we’ll obtain a solution that can simultaneously solve a rich class of related
approximation problems, one of which is the particular one we are interested in at
present.
The problem that (3.62) poses in abstract vector space is the following. Given
a particular element x in a vector space, how do we optimally approximate this
element with an arbitrary linear combination of other elements y1 , y2 , . . . , yM in
this space. Note that all linear combinations of these approximating elements con-
stitute a subspace of the vector space, which can be thought of as a hyperplane.
Thus, our problem is to choose the element in this hyperplane that best approxi-
mates our element x, where x itself generally does not lie in this hyperplane.
122 Estimation Theory Chap. 3

( x-x )

Figure 3.1. Approximation by projec-

tion in R2 .

Based on our intuition about such approximations in RN , we would expect

that the optimum approximation x̂ would correspond to a projection of x onto the
hyperplane, which would make the approximation error orthogonal to the hy-
perplane. Fig. 3.1 depicts this result in R2 when there is M = 1 approximating
element.
In fact, this intuition and reasoning can be used to construct optimum ap-
proximations in arbitrary abstract vector spaces. The generalization and optimal-
ity of this result is given in terms of the celebrated Orthogonal Projection Theorem,
which we now develop.

Theorem 3.2 (Orthogonal Projection) Let Y be a subspace of a (complete) inner prod-

uct space V, and let x ∈ V be an arbitrary element. Then
x̂ = arg min kw − xk, (3.63)
w∈Y

i.e., the approximation x̂ ∈ Y minimizes kx̂ − xk, if

e = x̂ − x ∈ Y⊥ ,
i.e., if
(x̂ − x) ⊥ w for all w ∈ Y. (3.64)

In order to simplify our proof of Theorem 3.2, let us first establish Pythagoras’
Theorem: two elements v1 and v2 in an inner product space V are orthogonal, i.e.,
v1 ⊥ v2 , if and only if
kv1 + v2 k2 = kv1 k2 + kv2 k2 . (3.65)
Sec. 3.2 Estimation of Random Vectors: a Bayesian Formulation 123

To derive this intermediate result, it suffices to note

kv1 + v2 k2 = hv1 + v2 , v1 + v2 i
= hv1 , v1 i + hv1 , v2 i + hv2 , v1 i + hv2 , v2 i
= kv1 k2 + kv2 k2 + 2 hv1 , v2 i . (3.66)
Clearly, the right-hand sides of (3.66) and (3.65) are equal if and only if v1 and v2
are orthogonal.
Returning now to our proof of Theorem 3.2, since x̂ − x is orthogonal to
every element of Y, it is orthogonal to the particular element w − x̂ for w ∈ Y, i.e.,
(x̂ − x) ⊥ (w − x̂). Hence, by Pythagoras’ Theorem we have
kx − wk2 = k(x̂ − x) + (w − x̂)k2 = kx̂ − xk2 + kw − x̂k2
from which we can conclude
kx − wk2 ≥ kx̂ − xk2
with equality if and only if w = x̂. Thus x̂ is the solution to (3.63).
Note that the theorem yields a remarkably general result. Not only does it
not depend on the choice of inner product, but it doesn’t depend on the dimen-
sion of Y either. Indeed, Y could have infinite dimension, and in fact the infinite
dimensional case will be an important focus later in the course. However, in the
estimation problem we consider in this section, i.e., (3.62), we’ll assume the corre-
sponding subspace has finite dimension, i.e., M < ∞.
In the finite dimensional case, the Orthogonal Projection Theorem can be
used to obtain a matrix representation for the optimum approximation x̂. To see
this, when Y is M-dimensional and spanned by the elements
 T
y = y1 y2 · · · yM ,
any element of Y can be expressed as a linear combination of the yi ’s, and in par-
ticular, since x̂ ∈ Y, we can write
M
X
T
x̂ = â y = âi yi (3.67)
i=1

for an appropriate choice of weights

 T
â = â1 â2 · · · âM .
Let’s now use the orthogonality condition (3.64) to determine what these weights
must be. In particular, since x̂ must satisfy
hx − x̂, yj i = 0, j = 1, 2, . . . , M, (3.68)
we can substitute (3.67) into (3.68) to obtain
* M ! + M
X X
hx, yj i = âi yi , yj = âi hyi, yj i j = 1, 2, . . . , M. (3.69)
i=1 i=1
124 Estimation Theory Chap. 3

Collecting the M equations of (3.69) into vector form we obtain the normal equations
RT
xy = Ryy â, (3.70)
where
 T
â = â1 â2 · · · âM (3.71)
 
Rxy = hx, y1 i hx, y2 i · · · hx, yM i (3.72)
 
hy1 , y1 i hy2 , y1 i ··· hyM , y1 i
 hy1 , y2 i hy2 , y2 i ··· hyM , y2 i 
Ryy =  . (3.73)
 
.. .. .. ..
 . . . . 
hy1 , yM i hy2 , yM i · · · hyM , yM i

The matrix Ryy in (3.73) is referred to a the Grammian matrix associated with
the approximation problem. It is straightforward to verify that the Grammian is
a positive semidefinite matrix.6 Furthermore, if the y1 , y2 , . . . , yM are a basis for Y
(i.e., are a linearly independent set), then the Grammian is strictly positive definite
and hence invertible. In this case, the optimal weights are given by
â = R−1 T
yy Rxy .

When the yi ’s are not linearly independent, there is still a solution, but optimal
weights â are no longer unique. In particular, any solution of (3.70) will be optimal.
Again, we stress that this framework is remarkably general and can be used
to solve a host of approximation problems. For example, when we let V = L2 (R)
or V = `2 (Z), this framework solves the continuous- or discrete-time deterministic
linear least-squares approximation problem.7 Estimation problems can be viewed
as a specific class of approximation problems involving random variables. In par-
ticular, when V = L2 (Ω), the framework solves the linear least-squares estimation
problem that is of primary interest in this section, which we now develop.

6
Let z = aT y where  T
a = a1 a2 ··· aM
and note that
0 ≤ kzk2 = hz, zi = aT Ryy a.

7
A prototypical deterministic least-squares problem involves approximating some deter-
ministic known but arbitrary function x(t) as a linear combination of other deterministic known
but fixed functions y1 (t), y2 (t), . . . , yM (t), i.e.,
M
X
x̂(t) = ai yi (t)
i=1

so as to minimize the energy in the approximation error, i.e.,

Z +∞
[x̂(t) − x(t)]2 dt.
−∞

As an example, the approximating functions could be polynomials, e.g., yi (t) = ti .

Sec. 3.2 Estimation of Random Vectors: a Bayesian Formulation 125

In this case, we choose the particular inner product and associated norm
given by (3.60) and (3.61), respectively, and let our approximating elements be the
random variables (3.59). Now any affine estimator for x based on y can, without
loss of generality, be expressed in the form
x̂(y) = d + aT (y − my ). (3.74)
In order to accommodate the additional constant term, we need to augment our
observed data with one additional (deterministic) random variable that is the con-
stant 1, so our data is now  
1
ỹ+ = (3.75)
ỹ
with
ỹ = y − my . (3.76)
With this notation our estimator (3.74) takes the form
x̂(y) = aT
+ ỹ+ , (3.77)
where  
d
a+ = . (3.78)
a
Proceeding, the corresponding normal equations (3.70), i.e.,
T
 
E [x ỹ+ ] = E ỹ+ ỹ+ â+ ,
become
dˆ
    
mx 1 0
T =
Λxy 0 Λy â
from which we get the two equations
dˆ = mx (3.79)
and
ΛT
xy = Λy â. (3.80)
Substituting (3.79) and (3.80) into (3.74) we obtain, when Λy is nonsingular,
x̂LLS (y) = mx + âT ỹ = mx + Λxy Λ−1
y ỹ. (3.81)
Finally, substituting (3.76) into (3.81) gives
x̂LLS (y) = mx + Λxy Λ−1
y (y − my ). (3.82)

Turning next to the performance of the resulting LLS estimator, we first note
from (3.82) that the estimator is unbiased. This follows immediately from the fact
that the estimation error must be orthogonal to the constant 1, which was an ele-
ment of the data vector (3.75), i.e.,
E [(x̂LLS (y) − x) · 1] = E [x̂LLS (y) − x] = 0, (3.83)
126 Estimation Theory Chap. 3

which in turn means of course that the mean-square estimation error is the same
as the error variance.
Second, we note that the error variance associated with this estimator can
also be obtained in a particularly straightforward manner using the orthogonality
condition. In particular, via Pythagoras’ Theorem we have that

λLLS = λe = E (x − x̂LLS (y))2

 

= kx − x̂k2 = kxk2 − kx̂k2

= E x 2 − kmx + âT ỹk2
 

= λx − âT Λy â
= λx − Λxy Λ−1 T
y Λxy , (3.84)
where to obtain the last equality in (3.84) we have used (3.80). This completes our
derivation of the scalar LLS estimator.
Our results are easily extended to the case in which we want to construct the
linear least-squares estimate of an N-vector x based on observations y. To do this,
we solve the problem in a component-wise manner, constructing the optimum
linear least-squares solution for the estimation of each component of x. The result
is a set of estimates
x̂i,LLS (y) = mxi + Λxi y Λ−1
y (y − my ) i = 1, 2, . . . , N, (3.85)
which we collect into vector form as
x̂LLS (y) = mx + Λxy Λ−1
y (y − my ). (3.86)

Exploiting that this estimator is unbiased, the error covariance can be readily
calculated. In particular, with the error written as
e = x̂(y) − x = Λxy Λ−1
y ỹ − x̃

where x̃ = x − mx , we obtain
ΛLLS = E eeT
 
h
T i
= E Λxy Λ−1 −1


y ỹ − x̃ Λ Λ
xy y ỹ − x̃
= Λx − Λxy Λ−1 T
y Λxy . (3.87)

As a final comment before we proceed to an example, note that as we an-

ticipated both the LLS estimator (3.86) and its performance (3.87) depend only on
the joint second-order statistics of x and y. As a consequence, this means that we
do not need to know the complete statistical characterization of x, y and their re-
lationship in order to construct and evaluate this estimator. On the other hand,
we must also keep in mind that this implies the LLS estimator cannot exploit the
higher-order statistical dependencies among x and y to better estimate x. As a re-
sult, our LLS estimators can never give us a lower mean-square estimation error
than our BLS estimators, which do fully exploit this additional information.
Sec. 3.2 Estimation of Random Vectors: a Bayesian Formulation 127

Example 3.5
Let’s consider the random variables x and y from Example 3.4 again, but now find
the LLS estimator for x based on y . First we note that by symmetry
mx = my = 0. (3.88a)
Furthermore, since x and w are independent we have
λxy = E [xy ] = E [x(sgn x + w )] = E [|x|] = 1/2, (3.88b)
and
λy = E y 2 = E (sgn x + w )2 = E sgn2 x + E w 2 = 1 + 1/3 = 4/3.
       
(3.88c)

Substituting (3.88) into (3.86) and (3.87) we obtain

λxy 3
x̂LLS (y) = y= y (3.89)
λy 8
and
λ2xy 1 (1/2)2
λLLS = λx − = − = 7/48. (3.90)
λy 3 (4/3)
Comparing (3.90) with (3.41) we see that, as expected, constraining our estimator to
be linear leads to a larger mean-square estimation error.

Additional Properties of LLS Estimators

We’ve already established some important properties of LLS estimators. For in-
stance we’ve shown that the LLS estimator is always unbiased [see (3.83)]. In this
section, we develop several additional special properties of LLS estimators.
In Section 3.2.4 we showed in Theorem 3.1 that an estimator x̂(y) is the Bayes’
least-squares estimator if and only if the corresponding estimation error x̂(y) −
x is orthogonal to any function of the observed data y. From the orthogonality
principle we used to derive the linear least-squares estimator, we immediately
obtain the counterpart to Theorem 3.1 for linear estimators.

Theorem 3.3 A linear estimator x̂L (·) is the linear least-squares estimator, i.e., x̂L (·) =
x̂LLS (·), if and only if the associated estimation error e(x, y) = x̂L (y) − x is orthogonal to
any vector-valued linear (i.e., affine) function of the data, i.e.,
h i
T
E [x̂L (y) − x] [Fy + g] = 0 (3.91)

for any constant matrix F and any constant vector g.

The proof follows immediately from the fact that the estimation error is orthogonal
to any of the elements of (3.75), and thus any linear combination of these elements
as well.
128 Estimation Theory Chap. 3

As another property of LLS estimators, we have that the covariance of the

BLS and LLS estimators are related according to the matrix inequalities
0 ≤ ΛBLS ≤ ΛLLS ≤ ΛL (3.92)
where ΛL is the error covariance of any linear (i.e., affine) estimator x̂L (·), and
where the rightmost inequality is satisfied with equality if and only if
x̂L (y) − E [x̂L (y) − x] = x̂LLS (y). (3.93)

Before deriving (3.92), we make several observations. First, the leftmost in-
equality in (3.92) is merely a restatement of the fact that covariance matrices are
positive semidefinite, and equality holds when x can be determined with certainty
from the data—this is what is called the “singular estimation” scenario.
Second, the middle inequality in (3.92) is an immediate consequence of (3.46)
which holds for any estimator x̂(·) and therefore any linear estimator x̂L (·). Fur-
thermore, as we will see shortly, this middle inequality in (3.92) is satisfied with
equality when x and y are jointly Gaussian. However, the converse is not true:
there do exist non-Gaussian examples where the BLS estimator turns out to be a
linear estimator.
Finally, since x̂L (y) = mx is a valid linear estimator, and has an associated
error covariance of ΛL = Λx , we have as a special case of (3.92) the statement
ΛLLS ≤ Λx . (3.94)
From (3.86), we see that (3.94) is satisfied with equality if and only if x and y are
uncorrelated.
We derive the rightmost inequality in (3.92) by following an approach analo-
gous to that used to derive the corresponding result for BLS estimators, i.e., (3.46).
In particular, let bL denote the bias in our estimator x̂L (·), let
h(y) = x̂L (y) − x̂LLS (y) − bL , (3.95)
and let us begin by noting that
h i
ΛL = E (x̂L (y) − x − bL ) (x̂L (y) − x − bL )T
h i
= E [h(y) + (x̂LLS (y) − x)] [h(y) + (x̂LLS (y) − x)]T
h i
 T
 T
= E h(y)h (y) + E (x̂LLS (y) − x) (x̂LLS (y) − x)
T
+ E (x̂LLS (y) − x) hT (y) + E (x̂LLS (y) − x) hT (y) .
  
(3.96)
Since h(y) is a linear (i.e., affine) function of y, from Theorem 3.3 we get that the
last two terms in (3.96) are zero. Using this together with the definition of ΛLLS we
get
ΛL − ΛLLS = E h(y)hT (y) .
 
(3.97)
Sec. 3.2 Estimation of Random Vectors: a Bayesian Formulation 129

The right-hand side of (3.97) is in general positive semidefinite, which verifies the
rightmost inequality in (3.92), and equal to zero if and only if h(y) = 0, which,
using (3.95), yields (3.93).
We finish this section with two examples of linear least-squares estimation.

Example 3.6
Consider the scalar problem of estimating a random variable x whose mean is mx
and whose variance is σx2 based on observations of the form

y = hx + w (3.98)

where h is a known deterministic constant, and where w has zero mean, variance
σw2 , and is independent of x. In this case, to construct the LLS estimator (3.86) we
need only determine the appropriate statistics. In particular, we have

λxy = hσx2 (3.99)

σy2 = h2 σx2 + σw2 (3.100)
my = hmx (3.101)

so that the LLS estimator is

hσ 2
x̂LLS (y) = mx + 2 2 x 2 (y − hmx )
h σx + σw
σw2 h2 σx2
   
y
= 2 2 2
mx + 2 2 2
(3.102)
h σx + σw h σx + σw h

and the associated error covariance is

h2 σx4 σx2 σw2
λLLS = σx2 − = . (3.103)
h2 σx2 + σw2 h2 σx2 + σw2

Example 3.7
Next let’s consider the vector generalization of Example 3.6, which arises in a host
of practical problems. Specifically, suppose that x has mean mx and covariance Λx
and that our observations are noisy measurements of linear functions of x, i.e.,

y = Hx + w (3.104)

where H is a given matrix and w has, for convenience, zero-mean, covariance Λw ,

and is uncorrelated with x.8

8
There is no loss of generality in assuming that w has zero-mean; if the measurement noise
w were nonzero mean, we could simply subtract this mean from y to obtain an equivalent problem
with zero-mean. Also, it is certainly reasonable to consider a scenario in which w is correlated with
x, though slightly more complicated expressions result in this case.
130 Estimation Theory Chap. 3

To construct the linear least-squares estimator for this problem simply requires
that we determine the appropriate statistics in (3.86). In particular, we have
my = E [y] = HE [x] + E [w] = Hmx (3.105)
Λxy = E (x − mx )(y − my )T
 

= E (x − mx )(H(x − mx ) + w)T
 

= E (x − mx )(x − mx )T HT + E (x − mx )wT
   

= Λx HT , (3.106)
= E (y − my )(y − my )T
 
Λy
= E [H(x − mx ) + w] (x − mx )T HT + wT
  

= HΛx HT + Λw . (3.107)
Then, from (3.86), (3.87) we obtain that
x̂LLS (y) = mx + K (y − Hmx ) (3.108)
T T


ΛLLS = Λx − K HΛx H + Λw K , (3.109)
where K is a gain matrix defined as
K = Λx HT (HΛx HT + Λw )−1 . (3.110)
Note the intuitively appealing structure of (3.108)—our posterior estimate x̂LLS (y)
equals our prior estimate mx plus a correction term that is proportional to the dif-
ference between the observation y and our best prediction of the observation based
on the prior information, i.e., my = Hmx .
Note that the gain matrix K controls the relative weight placed on our prior
information versus the observation. In particular, as Λx increases (e.g., in the sense
of its trace), our prior information degrades in quality, and one would therefore
want to place more weight on y. If Λw increases (again in the sense of its trace), the
quality of the measurement decreases and we would want less weight placed on the
measurement. The gain matrix makes the tradeoff in a statistically optimal manner.
The gain matrix K also optimally captures interdependencies among the com-
ponents of the vector to be estimated. This is especially important in applications
where it is only possible to obtain measurements of some of the variables of inter-
est. For example, consider the estimation of vehicle position x1 and velocity x2 in a
tracking problem, and let
 
x
x= 1 .
x2
Furthermore, suppose we are only able to obtain measurements of the position, so
that
y = x1 + w .
In this case, the only way in which y helps us to estimate the velocity x2 is through its
correlation with position x1 . This dependency is exploited as efficiently as possible
via the gain matrix K.
As a final comment, an alternate expression for ΛLLS in (3.109) that is derived
in Appendix 3.A is (see (3.343))
Λ−1 −1 T −1
LLS = Λx + H Λw H. (3.111)
Sec. 3.2 Estimation of Random Vectors: a Bayesian Formulation 131

Which form is more useful for practical computations of ΛLLS depends on a num-
ber of factors. Note for example, that (3.111) involves inversions of a matrix of size
N , the dimension of x, while (3.109) involves inversion of a matrix of size M , the
dimensions of y. Depending on the relative sizes of M and N , one form may be
preferable to the other. Other considerations that influence the choice involve nu-
merical stability issues, which we won’t develop here.
In any case, the form (3.111) provides us with some valuable intuition. As
we’ll see shortly, the inverse of a covariance has a useful interpretation as a mea-
sure of information. What (3.111) states is that the information Λ−1LLS about x after
the measurement equals the prior information Λx plus the information HT Λ−1
−1
w H
contained in the measurement.

Estimation in the Jointly Gaussian Case

In this section, we establish yet another very special property of jointly Gaussian
random variables. In particular we have the remarkable result that if x and y are
jointly Gaussian random vectors, then
x̂BLS (y) = x̂LLS (y). (3.112)
To see this, let
eLLS = x̂LLS (y) − x, (3.113)
and note that by Theorem 3.3 we have that eLLS must be orthogonal to every linear
function of y and hence y itself. But since x and y are jointly Gaussian, this means
that e is actually statistically independent of y. This implies, for example, that
E [eLLS |y] = E [eLLS ] = 0 (3.114)
where the last equality follows from the fact that the LLS estimate is unbiased. But
we also have directly from (3.113) and from (3.33) that
E [eLLS |y] = E [x̂LLS (y)|y] − E [x|y] = x̂LLS (y) − x̂BLS (y). (3.115)
Comparing (3.114) and (3.115) completes our derivation.
Note that this result provides a convenient derivation of the mean and co-
variance associated with the posterior density px|y (x|y) in the jointly Gaussian
case, i.e., (1.149) and (1.150), respectively, from Chapter 1. To establish the mean
expression (1.149) it suffices to combine (3.112) with (3.33) and (3.86). To establish
the covariance expression (1.150) we first note, using (3.112) and the fact that eLLS
and y are jointly Gaussian, that
h i
T
Λx|y = E (x − E [x|y]) (x − E [x|y]) | y
= E eLLS eT
 
LLS | y
= E eLLS eT
 
LLS = ΛLLS . (3.116)
Combining (3.116) with (3.87) we get our desired posterior covariance expression.
Again we emphasize that the resulting posterior covariance (1.150) is not a func-
tion of y in this jointly Gaussian case.
132 Estimation Theory Chap. 3

Finally, to verify that the posterior density for x is actually Gaussian we need
only recognize that since
     
x mx Λx Λxy
∼N , T (3.117)
y my Λxy Λy
and since
px,y (x, y)
px|y (x|y) = ,
py (y)
we have that when viewed as a function of x alone (with y fixed), the posterior
density satisfies
 T  −1  !
1 x − mx Λx Λxy x − mx
px|y (x|y) ∝ exp − (3.118)
2 y − my ΛT
xy Λy y − my

which is an exponentiated quadratic in x. But the exponentiated quadratic family

of densities are precisely the Gaussian densities, and thus the posterior density is
Gaussian with the mean and covariance determined above. Obviously, results for
scalar x and y correspond to a special case of what we’ve obtained above.
Note that as a consequence of this result we have that although the BLS esti-
mator is in general a nonlinear function of the data, in the jointly Gaussian case it
turns out to be a linear function of the data. Furthermore, as we mentioned earlier
would be the case, from (3.112) we get immediately that the middle inequality in
(3.92), i.e.,
ΛBLS ≤ ΛLLS ,
is satisfied with equality in the jointly Gaussian case. In turn this implies that there
are no statistical dependencies among x and y beyond second-order ones that the
BLS can exploit in the jointly Gaussian case.
In summary then we have that when x and y are jointly Gaussian,
x̂BLS (y) = mx + Λxy Λ−1
y (y − my ) (3.119)
ΛBLS = E Λx|y = Λx|y = Λx − Λxy Λ−1 T
 
y Λxy . (3.120)
And we note that when x and y are uncorrelated (Λxy = 0), they are also inde-
pendent in this case, and the BLS estimator and its performance degenerate to the
prior statistics on x since y provides no information about x.
There are other consequences of the fact that the posterior density for x is
Gaussian. For example, since the Gaussian density is symmetric and unimodal, its
mean is also its mode. Since the posterior mean is the BLS estimator and the pos-
terior mode is the MAP estimator, we have that when x and y are jointly Gaussian
the two estimators coincide, i.e.,
x̂BLS (y) = x̂MAP(y). (3.121)
As a result, in the jointly Gaussian case, the MAP estimator inherits a variety of
useful properties. For example, in general MAP estimates are biased. However,
Sec. 3.2 Estimation of Random Vectors: a Bayesian Formulation 133

in the jointly Gaussian case they are unbiased. Note, however, that x and y need
not be jointly Gaussian for (3.121) to hold. Indeed, any px,y (x, y) such that the
corresponding posterior density px|y (x|y) is symmetric and unimodal, for instance,
will have this property.
We finish this section with a scalar example.
Example 3.8
Suppose that x and y are scalar, jointly Gaussian random variables with
     2 
x mx σx λxy
∼N , . (3.122)
y my λxy σy2
Then
px|y (x|y) = N (x; x̂BLS (y), λBLS ) , (3.123)
where
 
λxy σx
x̂BLS (y) = mx + (y − my ) = mx + ρxy (y − my ). (3.124)
σy2 σy
Furthermore,
λ2xy
λBLS = λx|y = σx2 − = σx2 (1 − ρ2xy ), (3.125)
σy2
where
λxy
ρxy =
σx σy
is the correlation coefficient.
Several observations should be re-emphasized. First, in the jointly Gaussian
case we were able to express the posterior variance as λBLS which does not depend
on y since λx|y (y) does not depend on the actual observed value of y, i.e.,
 
λBLS = E λx|y = λx|y . (3.126)
Second, note that because x and y are jointly Gaussian, when ρxy = 0 they are
also independent. In this case, y contains no information about x, which is reflected
in the fact that (3.124) and (3.125) reduce to the prior statistics on x. Conversely,
larger values of |ρxy | result in more weight being placed on information from the
measurement, and a reduction in the posterior uncertainty—i.e., the uncertainty in
x after incorporating knowledge of y .
Again we emphasize that because x and y are jointly Gaussian the resulting
estimator (3.124) is a linear (or more precisely affine) function of the data y, and that
in non-Gaussian cases BLS estimators are generally nonlinear functions of the data
y.
Finally, from our earlier comments, we note that the MAP estimator and BLS
estimators are identical in the jointly Gaussian case, so we immediately obtain
x̂MAP (y ) = x̂BLS (y ).
Note too that since the posterior density is symmetric and unimodal, its mean is also
its median. Since the median of the posterior density is the minimum absolute-error
estimator, we have as well
x̂MAE (y ) = x̂BLS (y ).
134 Estimation Theory Chap. 3

BLS Estimation with only Second-Moment Information

At the outset of Section 3.2.5, we showed that the construction of the BLS estimator
of a random variable x from a random vector y subject to the constraint that the es-
timator be linear requires only knowledge of the joint second-moment properties
of (x, y).
This observation raises an interesting related question. Suppose we only
have knowledge of the joint second-moment properties of a pair (x, y), then what
is the best possible estimator x̂(y) (in a BLS sense) we can construct, and how does
it perform? In particular, we might reasonably ask whether the LLS estimator is
also the solution to this problem.
To answer this question requires posing the problem as a game between two
adversaries: the system designer tries to find the best estimator, and nature tries to
find the model that makes the performance of the chosen estimator as bad as pos-
sible subject to the constraint that the (x, y) statistics match the prescribed second-
moment information.
For this game, we now determine the best estimator choice for the system
designer, the worst joint distribution we can encounter, and the resulting mean-
square estimator error. With the given moments being mx , my , σx2 , σy2 , and λxy , we
seek to evaluate
max min E (x − f (y))2 ,
 
(3.127)
pxy ∈M f (·)

where
          2 
x mx x x σx Λxy
M= pxy : E = , cov , = . (3.128)
y my y y ΛT
xy Λyy

We first note that, in general, a max-min can be upper bounded by a min-

max. For our problem, this means that
max min E (x − f (y))2 ≤ min max E (x − f (y))2 .
   
(3.129)
pxy ∈M f (·) f (·) pxy ∈M

Intuitively, the minimizer is more powerful on the left side of (3.129) while the
maximizer is more powerful on the right side. Indeed, the minimizer on the left
side of (3.129) gets to choose an estimating function f that depends on the distri-
bution chosen by the maximizer, while the opposite is true for the right side.
We can further upper bound the right side of (3.129) by substituting any func-
tion f0 . That is,
min max E (x − f (y))2 ≤ max E (x − f0 (y))2 .
   
(3.130)
f (·) pxy ∈M pxy ∈M

For any linear function f0 (y) = d+aT (y−my ) the right side of (3.130) only depends
on the second-moment statistics, which are fixed. Let us further choose d = mx and
aT = Λxy Λ−1
y , which results in
h
2 i
E (x − f0 (y))2 = E x − mx − Λxy Λ−1 = σx2 − Λxy Λ−1 T
 
y (y − m y ) y Λxy (3.131)
Sec. 3.2 Estimation of Random Vectors: a Bayesian Formulation 135

for any pxy ∈ M. Thus, the right side of (3.131) is an upper bound to (3.127).
We now proceed to show that the right side of (3.131) is also a lower bound
to (3.127). Similarly to our upper bound in (3.130), we can lower bound (3.127) by
choosing any distribution on x and y , i.e.,
max min E (x − f (y))2 ≥ min E (x − f (y))2
   
(3.132)
pxy ∈M f (·) f (·)

where the expectation on the righthand side of (3.132) is with respect to an arbi-
trary distribution p∗xy ∈ M. Let us choose p∗xy as that corresponding to x and y being
jointly Gaussian with the specified second-moment statistics. For jointly Gaussian
random variables, the BLS estimate is the linear estimate x̂ = mx +Λxy Λ−1 y (y−my ).
The resulting mean square error is that given on the right side of (3.131).
Since (3.127) is both upper and lower bounded by the right side of (3.131), we
conclude that 1) the system designer should choose the LLS estimator, 2) nature
should choose the jointly Gaussian model matching the second-moment statistics,
and 3) the resulting mean-square error performance will be σx2 − Λxy Λ−1 T
y Λxy .

Gram-Schmidt Orthogonalization

In many approximation problems in vector space, it turns out to be much more

convenient to work with approximating elements that are orthogonal to one an-
other, and normalized. Indeed, when the approximating elements form an ortho-
normal set, the normal equations take a particularly simple form since the Gram-
mian is then just the identity matrix!
Now any set of elements can be replaced with a corresponding set of ortho-
normal elements that span the same space or subspace. Furthermore, this new
orthonormal set can be generated from the original set of elements in a efficient
recursive manner. This is the essence of the Gram-Schmidt orthogonalization pro-
cess. You’ve probably seen this procedure used a linear algebra course for regular
vectors in RN . However, with our broader view of vector space in this course, we
now see that the appropriate abstraction of this procedure can be used with arbi-
trary vector spaces. In particular, if we apply this to some subspace of zero-mean
random variables, then the Gram-Schmidt procedure generates a set of uncorre-
lated random variables from a set of correlated ones. This notion will form the ba-
sis for some very efficient and powerful estimation algorithms we’ll discuss later
in the course.
To begin our discussion of the Gram-Schmidt process for arbitrary vector
spaces, we consider the general problem of approximating an element x ∈ V as a
linear combination of y1 , y2 , . . . , yM ∈ V, i.e., as in (3.67). As developed earlier, the
optimal values of the âi are obtained from the normal equations (3.70) by inverting
the Grammian, which is a computationally intensive task in general. In addition,
suppose that we add one more vector, yM +1 , to the set on which we will base our
approximation. In general, the optimum values of the coefficients â1 , â2 , . . . , âM
136 Estimation Theory Chap. 3

of the other yi ’s will all change—i.e., the procedure is not a simple recursive one.
Suppose, however, the yi are orthogonal, i.e., hyi , yj i = 0 if i 6= j. In this case, the
solution to (3.73) yields
hx, yii
âi = , i = 1, 2, . . . , M (3.133)
hyi , yi i
which implies that each âi can be calculated individually, in essence representing
the best approximation of x using that single element yi .
The preceding remarks suggest a powerful strategy for solving the normal
equations in general which is of particular importance in recursive approximation
in which the yi are received sequentially. The basic idea here is that if the yi ’s are
not orthogonal, we’ll transform them so that the resulting elements are orthogonal.
This procedure for accomplishing this, which we now describe, is referred to as
Gram-Schmidt orthogonalization.
First note that if we let ŷ[i|i − 1] denote the best linear approximation of
yi based on the preceding y’s, i.e., based on y1 , y2 , . . . , yi−1 , then the sequence
z1 , z2 , . . . obtained from the sequence y1 , y2 , . . . according to

z1 = y1 (3.134a)
zi = yi − ŷ[i|i − 1] i≥2 (3.134b)
has two important properties.
First, the z’s are orthogonal, i.e.,
hzi , zj i = 0 i 6= j (3.135)
To see this, assume, without loss of generality, that i > j. Then note that from
its definition (3.134b) as an approximation error and the Orthogonal Projection The-
orem, zi is orthogonal to y1 , y2 , . . . , yi−1 . However, zj is a linear combination of
y1 , y2 , . . . , yj with j ≤ i − 1. Therefore, (3.135) holds.
The second important property is
span(y1 , y2 , . . . , yk ) = span(z1 , z2 , . . . , zk ), k = 1, 2, . . . (3.136)
This implies that the best approximation that can be obtained in terms of the y’s
is the same as that which can be obtained in terms of the z’s. Since the latter are
orthogonal, the approximation in terms of the z’s is easier to compute.
To verify (3.136), we begin by noting that since zi is defined as a linear com-
bination of y1 , y2 , . . . , yi, we need only show the reverse, i.e., that each yi is also
a linear combination of z1 , z2 , . . . , zi . We’ll show this by mathematical induction.
First, from (3.134a) we see that (3.136) is trivially true for k = 1. Next, we assume
that (3.136) is true for k ≤ i − 1, and proceed to show that this implies it must be
true for k = i. In particular, using (3.134b) we see that
yi = zi + ŷ[i|i − 1] (3.137)
Sec. 3.2 Estimation of Random Vectors: a Bayesian Formulation 137

but ŷ[i|i − 1] is a linear combination of y1 , y2 , . . . , yi−1, which by the induction hy-

pothesis can also be written as a linear combination of z1 , z2 , . . . , zi−1 . This com-
pletes our proof.
The result just described suggests the following recursive approximation
procedure. Let x̂[i] denote the best approximation to x based on y1 , y2, . . . , yi or
equivalently, z1 , z2 , . . . , zi . We then begin by computing
hx, z1 i
x̂[1] = K1 z1 , K1 = (3.138)
hz1 , z1 i
At the next step, we receive y2 and first compute
hy2 , z1 i
z2 = y2 − γ21 z1 , γ21 = (3.139)
hz1 , z1 i
and then
hx, z2 i
x̂[2] = K1 z1 + K2 z2 , K2 = . (3.140)
hz2 , z2 i
We stress that once we have computed z2 , (3.140) tells us that it is very simple to
compute x̂[2] by updating our previous approximation, i.e., (3.140) states that
x̂[2] = x̂[1] + K2 z2 (3.141)
More generally, at step i we have
hx, zi i
x̂[i] = x̂[i − 1] + Ki zi , Ki = (3.142)
hzi , zi i
In general, of course, zi must be computed via
zi = yi − γi,1 z1 − γi,2z2 − · · · − γi,i−1zi−1 (3.143a)
with
hyi, zj i
γij = . (3.143b)
hzj , zj i
The preceding algorithm has a rather natural interpretation when our vector
space is the subspace of L2 (Ω) corresponding to zero-mean random variables. In
this case x and the y1 , y2 , . . . are all zero-mean random variables, and ŷ [i|i−1] is the
best estimate of yi given previous measurements of y1 , y2 , . . . , yi−1 . Moreover, zi is
the corresponding estimation error, i.e., we can think of zi as the one-step prediction
error in estimating yi. We can also think of zi as the new information in yi that is not
predictable based on observations of y1 , y2 , . . . , yi−1 . For this reason the sequence
z1 , z2 , . . . is often referred to as the innovations sequence. In the terminology of
stochastic processes to be developed in the next chapter of the notes, this innova-
tions sequence is an example of a white process, in that the sequence consists of
uncorrelated random variables, and the Gram-Schmidt procedure in this context
can be viewed as a whitening filter,—it takes the correlated sequence y1 , y2 , . . . as
input and produces the white sequence z1 , z2 , . . . as output.
138 Estimation Theory Chap. 3

We can also interpret the Gram-Schmidt procedure applied to random vari-

ables as another strategy for diagonalizing a covariance matrix. To see this, we
begin by observing that transformation from y ’s to z’s—and hence, via (3.136),
from z’s to y ’s—is linear. We can express this in matrix form as
   
y1 z1
 y2   z2 
 ..  = Γ  ..  . (3.144)
   
 .   . 
yM zM

for some matrix Γ. Note however, that because of the recursive nature of the algo-
rithm, the associated matrix Γ is lower triangular, and its coefficients in the lower
triangle are precisely the γij ’s we constructed via projections in (3.143b), i.e.,

1
 i=j
[Γ]ij = γij i > j

0 otherwise


Note that by substituting recursively for the z’s on the right-hand side of (3.143b),
we can also readily compute Γ−1 , which we see is also lower triangular. Specifi-
cally,
z1 = y1
z2 = y2 − γ21 z1 = y2 − γ21 y1
z3 = y3 − γ31 z1 − γ32 z2 (3.145)
= y3 − γ31 y1 − γ32 (y2 − γ21 y1 )
..
.

To complete our interpretation of Gram-Schmidt as a diagonalization, note

that using (3.144) we get that
Λy = E yyT = ΓΛz ΓT ,
 
(3.146)
and, in turn,
Λ−1
y = Γ
−T −1 −1
Λz Γ . (3.147)
Since the z’s are uncorrelated, this means that Λz is indeed diagonal. The factoriza-
tion (3.146) is generally referred to as an LDU (lower-triangular–diagonal–upper-
triangular) decomposition of a matrix, and as (3.147) indicates, once we have this
decomposition, it is particularly easy to compute Λ−1y . This interpretation provides
additional perspective on the the Gram-Schmidt procedure.
It is worth commenting that this is now the second method for diagonalizing
a covariance matrix that we’ve encountered. In Chapter 1 we developed a diago-
nalization based on the eigenvalue decomposition of the covariance matrix. There
are, in fact, others as well that are useful. What you should keep in mind is that
Sec. 3.2 Estimation of Random Vectors: a Bayesian Formulation 139

which diagonalization is most useful in practice depends greatly on the type of

problem being addressed.
As a final comment, we note that in many problems of interest, the computa-
tion of the z’s can be extremely simple. For example, in a large class of estimation
problems the observations yi are given explicitly as functions of the variable x to
be estimated. In such a case ŷ [i|i − 1] can be computed directly from x̂[i − 1]. We
illustrate this through the following example.

Example 3.9
Suppose that
yi = hi x + wi , i = 1, 2, . . . (3.148)
where the hi are known numbers, x has zero mean and variance σx2 , and the wi are
uncorrelated, zero-mean random variables with variances σi2 and are also uncorre-
lated with x.
We again let x̂[i] denote our optimum estimate of x based on observing y1 , y2 , . . . , yi ,
and denote the corresponding mean-square error in these estimates by σx2 [i]. With
this notation, we use σx2 [0] to denote the variance of x, i.e., the mean-square error
before any observation is made; hence σx2 [0] = σx2 .
To initiate the recursion, we consider the estimation of x based on the first
measurement y1 . This is just the scalar version of our LLS estimator problem, so
from (3.102) and (3.103) we obtain

x̂[1] = K1 y1 (3.149a)
σx2 [0]h1
K1 = , (3.149b)
h21 σx2 [0] + σ12
and the variance of the estimation error x̂[1] − x is
h21 λ2x [0] σ12 σx2 [0]
σx2 [1] = σx2 [0] − = . (3.150)
h21 σx2 [0] + σ12 h21 σx2 [0] + σ12
At the next step, we first need to compute the best estimate of y2 = h2 x + w2
based on y1 . However, since w2 is uncorrelated with y1 , we have that
ŷ [2|1] = h2 x̂[1], (3.151)
so
z2 = y2 − h2 x̂[1]. (3.152)
Note that no new estimates are needed to generate z2 . More generally,
zi = yi − hi x̂[i − 1], (3.153)
and furthermore, from (3.142) we have that

x̂[i] = x̂[i − 1] + Ki zi
= x̂[i − 1] + Ki [yi − hi x̂[i − 1]], (3.154)
where
λxzi
Ki = . (3.155)
λzi
 140 Estimation Theory Chap. 3

To calculate the statistics in (3.155), we obtain zi from (3.148) and (3.153) as

zi = hi (x − x̂[i − 1]) + wi . (3.156)
Then, since wi is uncorrelated with both x and y1 , y2 , . . . we see that
λzi = h2i σx2 [i − 1] + σi2 . (3.157)
Similarly, we obtain that

λxzi = hi E [x(x − x̂[i − 1])] + E [xwi ]

= hi E (x − x̂[i − 1])2 + hi E [x̂[i − 1](x − x̂[i − 1])]
 

= hi σx2 [i − 1], (3.158)

where we have used the fact that x and wi are uncorrelated, and the fact that the
error x̂[i − 1] − x is uncorrelated with any linear combination of y1 , y2 , . . . , yi−1 and
thus in particular with x̂[i − 1]. Substituting (3.157) and (3.158) into (3.155) and
combining the result with (3.154) we obtain the recursion for our estimator, viz.,
x̂[i] = x̂[i − 1] + Ki (yi − hi x̂[i − 1]) (3.159a)
hi σx2 [i − 1]
Ki = . (3.159b)
h2i σx2 [i − 1] + σi2

Since the computation (3.159) involves the sequence of mean-square errors

σx2 [1], σx2 [2], . . .
we must also calculate these recursively. We obtain this recursion by
direct computation. In particular, using (3.148) and (3.159) we obtain

σx2 [i] = E (x − x̂[i])2

 
h i
= E [(1 − Ki hi )(x − x̂[i − 1]) − Ki wi ]2
= (1 − Ki hi )2 σx2 [i − 1] + Ki2 σi2
h2i σx2 [i − 1]2
= σx2 [i − 1] −
h2i σx2 [i − 1] + σi2
σi2 σx2 [i − 1]
= . (3.160)
h2i σx2 [i − 1] + σi2
What we have just derived is a simple example of a Kalman filter. The general
development and study of this and related topics forms the focus of a major portion
of the advanced graduate subject in recursive estimation, 6.433.

3.3 NONRANDOM PARAMETER ESTIMATION

In many types of applications, we are interested in estimating certain parameters

of the observed data. For example, given noisy measurements of a sinusoid, we
might be interested in estimating the frequency of the sinusoid. In such cases it is
often unnatural or inappropriate to view these parameters as random. Rather, it
makes more sense to view them as deterministic quantities, but quantities that are
nevertheless unknown.
Sec. 3.3 Nonrandom Parameter Estimation 141

As an example, suppose we have a sequence of independent identically dis-

tributed Gaussian random variables y1 , y2 , . . . , yN , where the mean m and variance
σ 2 that parameterize the density are unknown. We’ll explore some approaches
to the problem of developing good estimators for these kinds of nonrandom pa-
rameters, and in doing so, we’ll explore the kinds of performance criteria that are
typically used in evaluating such estimators.

We stress at the outset that our Bayesian framework can’t be adapted in any
straightforward way to handle nonrandom parameter estimators. To see this, con-
sider the scalar parameter case with a least-squares cost criterion. If we attempt to
construct an estimate x̂(y) via

x̂(·) = arg min E (x − f (y))2

 
(3.161)
f (·)

we see that we obtain a degenerate solution. In particular, noting that the expec-
tation in (3.161) is over y alone (since x is deterministic), we immediately obtain
that the right-hand side of (3.161) is minimized by choosing x̂(y) = x, and hence
the optimum estimator according to (3.161) depends on the very parameter we’re
trying to estimate!

Obviously, regardless of the performance criterion we choose we’ll want to

restrict our search to “valid” estimators, i.e., estimators that don’t depend explic-
itly on the parameters we’re trying to estimate. In this portion of the notes, we’ll
develop some ways of thinking about and approaching the problem of finding
valid estimators that yield good performance.

In our treatment, we will use x to denote the vector of parameters we wish

to estimate, and write the density for the vector of observations y as py (y; x) so as
to make the parameterization explicit. In addition, we will use my (x) and Λy (x)
to denote, respectively, the mean vector and covariance matrix of y, again to make
the parameterization explicit.

3.3.1 Bias and Error Covariance

As in the case of random parameters, two important measures of the performance

of an estimator for nonrandom parameters are the bias and error covariance. How-
ever, there are some important distinctions between these quantities in the nonran-
dom case, which will become apparent in this section.

Using
e(y) = x̂(y) − x = x̂ − x (3.162)
142 Estimation Theory Chap. 3

as our notation for the error, we define the bias in an estimator x̂(·) as
bx̂ (x) = E [e(y)] = E [x̂(y) − x]
Z +∞
= [x̂(y) − x] py (y; x) dy
−∞
Z +∞ 
= x̂(y) py (y; x) dy − x (3.163)
−∞

Likewise, we express the error covariance as

h i
T
Λe (x) = E [e(y) − bx̂ (x)] [e(y) − bx̂ (x)] , (3.164)
where, again, the expectation is with respect to y.
We stress that both the bias (3.163) and error covariance (3.164) are, in gen-
eral, functions of the parameter x. Moreover, since the parameter x is unknown,
removing a bias bx̂ (x) as we discussed in the case of random parameter estimators
is generally not feasible in the nonrandom case—what bias to subtract would be a
function of the very quantity we wish to estimate. In other words, modifying such
an estimator by subtracting the bias will render it an “invalid” estimator.
In general the error correlation (and, in turn, its trace—the mean-square esti-
mation error) depends on both bias and error covariance; specifically
E e(y)eT (y) = Λe (x) + bx̂ (x)bT
 
x̂ (x) (3.165)
From this expression we see that we may not want simply to minimize Λe (x) if
this leads to a large bias. To illustrate this point, suppose, for example, we take as
our estimate x̂(y) a constant vector independent of y. In this case, Λe (x) = 0, but
the bias could be arbitrarily large.
For these reasons, a reasonable approach to developing good estimators for a
nonrandom parameter is to explicitly restrict our search for estimators to those that
are valid and unbiased,9 and among this class, choose the one having the smallest
variance. This is the notion underlying minimum-variance unbiased estimators,
which we discuss next.
As one final comment before we explore this topic, note that in contrast to the
case of random parameters, for nonrandom parameter estimators we have that the
error covariance is the same as theh covariance of the estimator itself,i i.e.,
Λe (x) = Λx̂ (x) = E (x̂(y) − E [x̂(y)]) (x̂(y) − E [x̂(y)])T .
To see this in the scalar parameter case, simply note that using (3.162), we have
λe (x) = E e2 (y) = E ((x̂(y) − x) − (E [x̂(y)] − x))2
   

= E (x̂(y) − E [x̂(y)])2
 

= λx̂ (x).

9
We say an estimator x̂(·) for a nonrandom parameter x is unbiased if bx̂ (x) = 0 for all
possible values of x.
Sec. 3.3 Nonrandom Parameter Estimation 143

We begin by considering the case where the parameter to be estimated is

a scalar x, which simplifies our exposition much as it did for the case of random
parameter estimation. We again note in advance that what we develop will readily
generalize to the vector case, since the estimation of a vector of parameters can be
accomplished in a component-wise manner.

3.3.2 Minimum-Variance Unbiased Estimators

To begin, let A denote the set of all estimators that are valid (i.e., don’t depend on
x) and unbiased, i.e.,
A = {x̂(·) | x̂(·) is valid and bx̂ (x) = 0}
Then, when it exists, a minimum-variance unbiased (MVU) estimator for x is de-
fined to be the estimator in A with the smallest variance, i.e.,
x̂MVU (·) = arg min λx̂ (x) for all x (3.166)
x̂∈A

Several observations regarding (3.166) are worth emphasizing. The first is

that x̂MVU (·) may not exist! For example, for some problems the set A is empty—
there are no valid unbiased estimators. In other cases, A is not empty, but no
estimator in A has a smaller variance than all the others for all values of the parameter
x. Suppose for example that A consists of three estimators x̂1 (·), x̂2 (·), and x̂3 (·),
whose variances are plotted as a function of the unknown parameter x in Fig. 3.2.
In this case, there is no estimator having a smaller variance than all the others for
all values of x.
It should also be emphasized that even when x̂MVU (·) does exist, it may be
difficult to find. In fact in general there is no systematic procedure for either de-
termining whether an MVU estimator exists, or for computing it when it is does
exist. However, fortunately there are cases in which such estimators can be com-
puted, as we’ll discuss later. As an example, we can determine the MVU estimator
in the linear/Gaussian case. Also, when we further restrict our attention to linear
MVU estimators, we’ll see that these can often be computed as well—in particular,
whenever the mean of the observations is a linear (affine) function of the parame-
ter.
Sometimes, it is useful to exploit a bound on λx̂ (x) in our quest for MVU
estimators. A particular useful bound for this purpose is the Cramér-Rao bound,
which we explore next.

3.3.3 The Cramér-Rao Bound

Again we consider first the estimation of an unknown scalar parameter x given a

measurement vector y with density py (y; x). When it exists, the Cramér-Rao bound
144 Estimation Theory Chap. 3

λ x (x)

λx1(x)

λx2(x)

λx3(x)

x
Figure 3.2. The variances of three un-
biased estimators.

gives a lower bound on the variance of any valid unbiased estimator x̂(·) for x. In
particular, the Cramér-Rao bound for any x̂(·) ∈ A is
1
λx̂ (x) ≥ , (3.167)
Iy (x)
where the nonnegative quantity Iy (x) is referred to as the Fisher information in y
about x, which is defined by
" 2 #
∂
Iy (x) = E ln py (y; x) . (3.168)
∂x
Some preliminary remarks are worth making. First, we stress that the Fisher
information cannot be computed in all problems, in which case no Cramér-Rao
bound exists. For example, for densities such as
(
1 x<y < x+1
py (y; x) = ,
0 otherwise
which are not strictly positive for all x and y, the logarithm in (3.168) doesn’t exist
and hence Iy (x) can’t be calculated.
Second, the notion of referring to (3.168) as an information measure comes
from the fact that Iy (x) is both nonnegative and additive, i.e., whenever
 T
y = y1 y2 · · · yM
consists of mutually independent components we have
N
X
Iy (x) = Iyi (x).
i=1
Sec. 3.3 Nonrandom Parameter Estimation 145

Third, the Fisher information (3.168) can be interpreted as a measure of cur-

vature: it measures, on average, how “peaky” ln py (y; x) is as a function of x. As
such, the larger Iy (x), the better we expect to be able to resolve the value of x from
the observations, and hence the smaller we expect λx̂ (x) to be. We’ll develop this
interpretation further in Section 3.4.

Example 3.10
Consider the scalar Gaussian problem
y = x + w,
where w ∼ N (0, σ 2 ). Then
1 1
ln py (y; x) = − (x − y )2 − ln(2πσ 2 ) (3.169)
2σ 2 2
Here the Fisher information is
1
Iy (x) = ,
σ2
so the smaller the variance σ 2 the sharper the peak of (3.169) is as a function of x.

To derive the Cramér-Rao bound (3.167), we begin by recalling that for unbi-
ased estimators the error
e(y) = x̂(y) − x (3.170)
has zero mean, i.e.,
E [e(y)] = 0, (3.171)
and variance
var e(y) = E e2 (y) = λx̂ (x).
 
(3.172)
Next we define
∂
f (y) = ln py (y; x) (3.173)
∂x
and note that using the identity
∂ 1 ∂
ln py (y; x) = py (y; x), (3.174)
∂x py (y; x) ∂x
we get that f (y) has zero mean:
 
1 ∂
E [f (y)] = E py (y; x)
py (y; x) ∂x
Z +∞
∂
= py (y; x) dy
−∞ ∂x
Z +∞
∂ ∂
= py (y; x) dy = 1 = 0, (3.175)
∂x −∞ ∂x
and, in turn, variance
var f (y) = E f 2 (y) = Iy (x).
 
(3.176)
146 Estimation Theory Chap. 3

Finally, again using the identity (3.174), the covariance between e(y) and f (y) is
given by

cov (e(y), f (y)) = E [e(y)f (y)]

Z +∞
∂
= (x̂(y) − x) py (y; x) dy
−∞ ∂x
 Z +∞   Z +∞ 
∂ ∂
= x̂(y) py (y; x) dy − x py (y; x) dy
∂x −∞ ∂x −∞
= 1 − 0 = +1. (3.177)
Now recall from from Chapter 1 that the correlation coefficient associated with
e(y) and f (y) satisfies, via the Cauchy-Schwarz inequality, the bound
[cov (e(y), f (y))]2
ρ2ef = ≤ 1. (3.178)
var e(y) var f (y)
Finally, substituting (3.172), (3.176) and (3.177) into (3.178) we get the Cramér-Rao
bound (3.167).
Several additional comments regarding the bound provide important in-
sights. First, the bound (3.167) in general depends on x, which we don’t know.
However, by plotting it as a function of x, we get a sense for the relative difficulty
of estimating x as a function of its true value. In addition, we can extract best- and
worst-case scenarios.
Second, any estimator that satisfies the Cramér-Rao bound with equality
must be a MVU estimator. Note however, that the converse is not true: the Cramér-
Rao bound may not be tight. Sometimes no estimator can meet the bound for all
x, or even for any x! An estimator which achieves the bound, i.e., satisfies (3.167)
with equality is referred to as an efficient estimator. Hence, efficient estimators are
MVU estimators, but the converse need not be true. We’ll return to a discussion of
efficiency shortly.
Third, the Cramér-Rao bound can be generalized in a variety of ways. For
example, a Cramér-Rao bound can be constructed for biased estimates. However,
in practice this bound is not particularly useful. Likewise, there is an analogous
bound for random parameters (see, e.g., Van Trees). However, this bound is not
widely used, primarily because we always have a tight bound on the error vari-
ance of random parameter estimates, viz.,
 
var e(x, y) = var [x̂(y) − x] ≥ E λx|y (y)
with equality if and only if x̂(y) = x̂BLS (y) = E [x|y].
In practice, the Fisher information (3.168) is frequently more useful when
re-expressed in the following form
 2 
∂
Iy (x) = −E ln py (y; x) . (3.179)
∂x2
Sec. 3.3 Nonrandom Parameter Estimation 147

To verify (3.179), we begin by observing

Z +∞
py (y; x) dy = 1. (3.180)
−∞

Differentiating (3.180) with respect to x and using the identity (3.174) yields
Z +∞
∂
py (y; x) ln py (y; x) dy = 0. (3.181)
−∞ ∂x
Finally, differentiating (3.181) once more with respect to x and again using (3.174)
we obtain
Z +∞  2  Z +∞  2
∂ ∂
py (y; x) ln py (y; x) dy + py (y; x) ln py (y; x) dy = 0, (3.182)
−∞ ∂x2 −∞ ∂x
which verifies that (3.168) and (3.179) are consistent.

Efficiency and Consistency

From our derivation of the Cramér-Rao bound (3.167) and in particular from (3.178),
we note that the Cramér-Rao bound is satisfied with equality if and only if the
functions e(y) and f (y) defined in (3.170) and (3.173), respectively, are perfectly
positively correlated, i.e., if and only if there exists some constant k(x) > 0 (i.e.,
that can only depend on x) such that
e(y) = k(x)f (y) for all y. (3.183)
As we mentioned earlier, we refer to estimators that satisfy the Cramér-Rao bound
with equality as efficient estimators. Rearranging (3.183) using (3.170) and (3.173),
we obtain that an efficient estimator x̂(·) must take the form
∂
x̂(y) = x + k(x) ln py (y; x). (3.184)
∂x
Hence, an efficient estimator exists if and only if (3.184) is a valid estimator, i.e., if
and only if the right-hand side of (3.184) is independent of x for some k(x).
However, k(x) cannot, in fact, be arbitrary. To see this, let us suppose that an
efficient estimator exists, so that (3.167) is satisfied with equality. Then, via (3.172)
we must have
1
E e2 (y) = λx̂ (x) =
 
. (3.185)
Iy (x)
Next note that using (3.183), (3.176), and (3.177) we obtain
E e2 (y) = E [e(y) · k(x)f (y)] = k(x)E [e(y)f (y)] = k(x)
 
(3.186)
Comparing (3.185) and (3.186), we can then conclude that
1
k(x) = . (3.187)
Iy (x)
148 Estimation Theory Chap. 3

Thus, substituting (3.187) into (3.184), we obtain the following characteriza-

tion for efficient estimators: an estimator x̂(·) is efficient if and only if it can be
expressed in the form
1 ∂
x̂(y) = x + ln py (y; x) (3.188)
Iy (x) ∂x
where the right-hand side must be independent of x for the estimator to be valid.
Three final remarks are important. First, note that an efficient estimator, i.e.,
a valid estimator satisfying (3.188) is guaranteed to be unbiased: taking the expec-
tation of (3.188) we get, using (3.175),
1
E [x̂(y)] = x + E [f (y)] = x.
Iy (x)
Second, we note that (3.188) implies that when it exists, an efficient estimator is
also unique—clearly no two estimators could satisfy (3.188) and be distinct. Finally,
since it meets a lower bound on the estimator variance, when it exists, an efficient
estimator must be the unique MVU estimator for a problem.
Let us now consider another desirable property of estimators. Suppose we
have an estimator for x based on a sequence of observations y1 , y2 , . . . , and let us
specifically denote by x̂M the estimate of x based on y1 , y2 , . . . , yM , i.e.,
x̂M = x̂(y1 , y2 , . . . , yM ).
Then we will say the estimator x̂M is a consistent estimator for x if
x̂M → x as M → ∞. (3.189)
Note that since x̂M is a random variable for each M we have to say what we
mean by the convergence condition (3.189), so let’s briefly discuss the four main
notions of convergence that are associated with sequences of random variables.
To make the discussion as general as possible, let z be an arbitrary random vari-
able and let z1 , z2 , . . . denote a sequence of related random variables. We note in
advance that z being a deterministic constant corresponds to a special case of our
discussion; z then has the degenerate density consisting of a single impulse at its
actual value.
The weakest form of convergence is termed “convergence in distribution.”
The sequence z1 , z2 , . . . is said to converge in distribution to z if, when Pz (·) de-
notes the distribution function for z and Pzn (·) denotes that for zn , we have
lim Pzn (z) = Pz (z)
n→∞

for all z at which Pz (z) is continuous. The notation

d
zn −→ z
is often used to specifically denote convergence in distribution. Note that this
does not say that the values of the random variables zn are getting close to the
Sec. 3.3 Nonrandom Parameter Estimation 149

value of z for large n, merely that their statistics are. This is, of course, the kind of
convergence that the Central Limit Theorem we discussed in Chapter 1 involves.
We also emphasize that convergence in distribution does not ensure convergence
in density, as was apparent in our discussion of the Central Limit Theorem in
particular.
A second form of convergence is termed “convergence in probability” or “p-
convergence.” We say that z1 , z2 , . . . converges in probability to z if for every fixed
 > 0 we have
lim Pr [|zn − z| > ] = 0.
n→∞

The notation
p
zn −→ z
is sometimes used to denote convergence in probability. This kind of convergence
is much stronger than convergence in distribution, and says something about the
actual values of the zn ’s converging to z. Convergence in probability implies con-
vergence in distribution, then, but of course the converse is not true. As an ex-
ample, the weak law of large numbers is a statement about the convergence in
probability of certain averages.
A still stronger notion of convergence is termed “mean-square convergence”
or “convergence in the mean.” We say z1 , z2 , . . . converges in mean-square (or “in
the mean”) to z if
lim E (zn − z)2 = 0.
 
n→∞

This kind of convergence is usually denoted using

m.s.
zn −→ z
or sometimes
l. i. m. zn = z.
n→∞

Using the Chebyshev inequality we discussed in Chapter 1, we can readily estab-

lish that convergence in mean-square implies convergence in probability. How-
ever, it is important to note (although it may not be obvious at first glance) that the
converse isn’t true.
Another very strong notion of convergence is termed “almost-sure” or “pro-
bability-1” convergence. We say z1 , z2 , . . . converges almost-surely (or with probability-
1) to z if
h i
Pr lim zn (ω) = z(ω), for all ω ∈ A ⊂ Ω, Pr [A] = 1 = 1.
n→∞

This kind of convergence is often denoted using

a.s.
zn −→ z,
and, while a technically somewhat difficult definition to digest, effectively re-
quires that (almost) every realization z1 , z2 , . . . of the sequence of random vari-
ables z1 , z2 , . . . converges to the corresponding realization z of z. It is this kind
150 Estimation Theory Chap. 3

of convergence that is involved in the strong law of large numbers. Almost-sure

convergence also implies convergence in probability, but again the converse is not
true.
Almost-sure convergence is the most desirable form of convergence in many
problems. However, it is often difficult to establish. By contrast, mean-square
convergence is often comparatively easier to work with, and is well-suited to en-
gineering problems. For this reason, we’ll generally restrict our attention to the
latter form of convergence. Keep in mind, however, that mean-square conver-
gence is not a weaker form of convergence than almost-sure convergence. In par-
ticular, mean-square convergence neither implies nor is implied by almost-sure
convergence. A full discussion of these issues and counterexamples is beyond the
scope of our treatment. However, good discussions can be found in a variety of
advanced probability texts.
Let us now return to our development of the notion of consistency of an es-
timator. Based on our discussion above, we’ll restrict our attention to consistency
in the sense of mean-square convergence of (3.189), so x̂M will be a consistent esti-
mator when
E (x̂M − x)2 → 0
 
as M → ∞. (3.190)
As a final comment before we explore a couple of examples, note that in
general there is no relationship between efficiency and consistency: an efficient
estimator need not be consistent, and a consistent estimator need not be efficient.

Example 3.11
Let’s continue with the linear Gaussian problem we began in Example 3.10, i.e.,
y = x + w, (3.191)
where w ∼ N (0, σ 2 ). In this case
√ 1
ln py (y; x) = − ln( 2πσ 2 ) − 2 (y − x)2 , (3.192)
2σ
so that
∂2 1
ln py (y; x) = − 2 , (3.193)
∂x2 σ
and from (3.179)
1
Iy (x) = . (3.194)
σ2
From the Cramér-Rao bound (3.167), we get that variance of any unbiased estimator
satisfies
λx̂ (x) ≥ σ 2 . (3.195)
Constructing the right-hand side of (3.188) using (3.192) and (3.194) we obtain
x̂(y) = y (3.196)
which we note is not a function of x and is therefore valid. Hence, we can immedi-
ately conclude that x̂ = x̂(y ) defined via (3.196) is unbiased and has a variance equal
to the Cramér-Rao bound, i.e.,
λx̂ (x) = σ 2 . (3.197)
Sec. 3.3 Nonrandom Parameter Estimation 151

Hence, we can conclude that (3.196) is an efficient estimator, and hence the unique
MVU estimator for the problem.

Example 3.12
Let’s consider a generalization of Example 3.11. In particular, suppose that we now
have a set of observations of x of the form
yi = x + wi i = 1, 2, . . . , M (3.198)
where the wi are independent identically-distributed random variables with densi-
ties N (0, σ 2 ). In this case,
M
−M 2 1 X
ln py (y; x) = ln(2πσ ) − 2 (yi − x)2 , (3.199)
2 2σ
i=1

and hence
M
∂ 1 X
ln py (y; x) = 2 (yi − x). (3.200)
∂x σ
i=1
From (3.200) and (3.168) we then obtain
M
.
Iy (x) = (3.201)
σ2
If we again construct an estimator from the right-hand side of (3.188) using
(3.200) and (3.201), we obtain
M
1 X
x̂(y) = yi (3.202)
M
i=1

which a valid estimator. Hence, (3.202) is unbiased and also an efficient estimator
for the problem, so its variance is
σ2
λx̂ = 1/Iy (x) = . (3.203)
M
Note too that our estimator (3.202) also happens to be consistent in this exam-
ple, i.e., from (3.203) we have
σ2
λx̂ = →0 as M → ∞.
M

3.3.4 Maximum Likelihood Estimation

To develop the topic of maximum likelihood estimators, we begin with the fol-
lowing observation regarding efficient estimators. Specifically, suppose an effi-
cient estimator exists for a particular problem of interest, and let x̂eff (·) denote this
estimator. Hence, for any particular value of the data y we have, rewriting (3.188),
1 ∂
x̂eff (y) = x + ln py (y; x). (3.204)
Iy (x) ∂x
152 Estimation Theory Chap. 3

which we can compute directly. Now since the right-hand side of (3.204) is inde-
pendent of the value of x, we are free to choose any value of x in this expression,10
so let us judiciously choose x to be the number
x̂ML (y) = arg max py (y; x). (3.205)
x

Since py (y; x) is typically referred to as the likelihood function of the data y, (3.205)
is referred to as the maximum likelihood (ML) estimator for x based on y.
From (3.205) we see that provided the likelihood function is strictly positive
and differentiable, the ML estimator satisfies
 
∂
ln py (y; x) = 0. (3.206)
∂x x=x̂ML (y)

Thus, since Iy (x) > 0 for all x except in the trivial case, (3.204) becomes
x̂eff (y) = x̂ML (y). (3.207)
From this we can conclude that when it exists, the (unique) efficient estimator is
equivalent to the ML estimator for the problem. For future convenience, we’ll use
λML (x) to denote the variance (and hence error variance) of the estimator (3.205).
However, several points should be stressed. This does not mean the ML
estimators are always efficient! When an efficient estimator doesn’t exist for a
problem, then the ML estimator need not have any special properties. This means,
for example, that when an efficient estimator does not exist, the ML estimator may
not have good variance properties or even be unbiased.
Nevertheless, ML estimators are highly practical—in particular, there exists
a systematic procedure for obtaining them from data. In problems where the like-
lihood, for a particular observed value of the data y, is a sufficiently tractable and
differentiable function of the parameter x, we may compute the ML estimate for
that y as follows. First, we analytically determine local maxima of the likelihood
function, i.e., solutions to
∂
py (y; x) = 0, (3.208)
∂x
for which
∂2
py (y; x) < 0.
∂x2
Then, we search over these local maxima and any boundary values for the largest
value of the likelihood function. In some problems, it often turns out to be easier to
maximize some monotonic function of the likelihood rather than the likelihood it-
self. For example, in a variety of problems maximizing the log-likelihood function
ln py (y; x) simplifies computations significantly.
It is worth pointing out, however, that the number of problems for which
solutions to (3.208) can be obtained as closed-form expressions is relatively small.

10
In particular, we need not choose x to be its true value.
Sec. 3.3 Nonrandom Parameter Estimation 153

More typically, iterative numerical techniques such as gradient searches (Newton-

Raphson) are used to find the local maxima of the likelihood function, from which
the global maximum is selected. In addition to general-purpose iterative ascent
algorithms, there also exist iterative ascent algorithms that are specifically tailored
to the special characteristics of likelihood functions. One class of these algorithms
are the so-called Estimate-Maximize (EM) algorithms, which have proven useful
in a wide range of practical estimation problems.
There are additional reasons why ML estimators have proven popular in
many applications even when they aren’t efficient estimators. For example, in
many cases these estimators have good asymptotic properties, i.e., when the size of
the vector y gets sufficiently large.
To develop the necessary concepts, as we did in Section 3.3.3, we again con-
sider the estimation of a parameter x based on a sequence of related observations
y1 , y2 , . . . , and we let x̂M denote the estimate based on the first M observations, i.e.,
on y1 , y2, . . . , yM . We then say x̂M is an asymptotically unbiased estimator if
E [x̂M ] → x as M → ∞. (3.209)
In addition, we say that an at least asymptotically unbiased estimator x̂M is weakly
asymptotically efficient if
λx̂M − Iy−1
1 ,y2 ,...,yM
(x) → 0 as M → ∞. (3.210)
Furthermore, we say that an estimator is strongly asymptotically efficient if it is
weakly asymptotically efficient and
Iy−1
1 ,y2 ,...,yM
(x) − Iy−1
1 ,y2 ,...,y∞
(x)
→1 as M → ∞. (3.211)
λx̂M − λx̂∞
Also, note that in general the concepts of consistency and asymptotic efficiency
need not be related. However, if for all x
Iy1 ,y2 ,...,yM (x) → ∞ as M → ∞,
then an estimator that is even weakly asymptotically efficient is also consistent.
In many problems, the ML estimator is not efficient but has the property that
it is asymptotically efficient and often consistent. Moreover, in a substantial subset
of such problems it is not only asymptotically efficient but also asymptotically
Gaussian; specifically
x̂M ∼ N (x, 1/Iy1 ,y2 ,...,yM (x)) as M → ∞.
While these “folk theorems” are often used casually to justify the optimality
of ML estimators, it should be emphasized that while they are often true, it is easy
to construct counterexamples for which the ML estimator is neither asymptotically
efficient nor even asymptotically unbiased. For example, as will become apparent
later in this chapter, for observations of the form
y1 = h(x) + w1 (3.212)
yi = wi i≥2 (3.213)
154 Estimation Theory Chap. 3

where h(·) is an invertible nonlinear function and where the wi are independent
identically-distributed N(0, σ 2 ) random variables, the ML estimator for x based on
y1 , y2 , . . . , yM for any M is
x̂ML (y) = h−1 (y1 ) (3.214)
However, for almost any choice of h(·) the ML estimator (3.214) is neither efficient
nor unbiased. Thus, since (3.214) is also independent of M, it is neither asymptot-
ically efficient nor unbiased either.
One class of problems for which ML estimators are always efficient and
therefore MVU estimators are the linear/Gaussian problems. We consider the
canonical scalar version of this problem in the following example.

Example 3.13
Consider the scalar linear/Gaussian problem
y = hx + w (3.215)
where w ∼ N (0, σw2 ). Note that
 
1 1
py (y; x) = N (y; hx, σw2 ) =p exp − 2 (y − hx)2
(3.216)
2πσw2 2σw

so that the ML estimator simply inverts h and ignores the noise, i.e.,
y
x̂ML (y) = . (3.217)
h
It is straightforward to verify this estimator is unbiased, i.e.,
 
hx + w 1
E [x̂ML (y ) − x] = E − x = E [w ] = 0 (3.218)
h h
and that its variance is
w2 σ2
 
λML (x) = E 2
= w2 (3.219)
h h
Note that in this case the estimator variance turns out to be independent of x.
Furthermore, the estimator variance is equal to the reciprocal of the Fisher informa-
tion for the problem, i.e.,
σ2
λML (x) = w2 = 1/Iy (x),
h
and therefore the ML estimator is efficient.
It is interesting to compare the ML estimator in this example to the LLS esti-
mator for the closely related problem developed in Example 3.6. In both examples,
the measurement models (3.215) and (3.98) are identical, but in this example x is a
nonrandom parameter while in Example 3.6 we have a random parameter x with
zero-mean and variance σx2 .
If we add the Gaussian assumptions to Example 3.6, we can conclude that the
resulting LLS estimator (3.102) is also the BLS estimator, the MAP estimator, and
the MAE estimator for the problem. For this reason, we’ll simply use x̂B (y) to de-
note this estimator for the remainder of this example. Furthermore, since our ML
Sec. 3.3 Nonrandom Parameter Estimation 155

estimate is efficient, it is the MVU estimator for the nonrandom parameter estima-
tion problem, so we’ll use x̂MVU (y) to denote this estimator for the remainder of this
example.
First, comparing (3.102) and (3.217) we see that
lim x̂B (y) = x̂MVU (y) (3.220)
σx2 →∞

which indicates that as our prior knowledge about x in the random parameter case
deteriorates (so that px (x) becomes increasingly flat) the Bayesian estimate approaches
the MVU estimate. In fact, from (3.102) we can see that the Bayesian estimate is a
linear combination of the best prior estimate mx and the MVU estimate y/h, where
the weights are determined by the relative quality of the prior information and the
measurement. Indeed, if we define a signal-to-noise ratio (SNR) of the form
mean-square contribution of “signal” portion of y h2 σ 2
SNR = = 2x (3.221)
mean-square contribution of noise in y σw
we have that    
1 SNR
x̂B (y) = mx + x̂MVU (y) (3.222)
1 + SNR 1 + SNR
Similarly, we can relate the performance of these estimators according to
1 1 1
= + 2, (3.223)
λB λMVU σx
to which we can attach the interpretation that the information after the measurement
equals the sum of the information in the measurement plus the prior information.

Let’s consider a couple of other examples of ML estimators that happen to

be efficient.

Example 3.14
Suppose that the random variable y is exponentially-distributed with unknown
mean x ≥ 0, i.e.,
1
py (y; x) = e−y/x u(y). (3.224)
x
Since py (y; x) and ln py (y; x) have the same maximum, we obtain the ML estimate
as the solution of
∂ ∂ h yi
ln py (y; x) = − ln x −
∂x ∂x x
1 y
= − + 2 = 0. (3.225)
x x
In particular, from (3.225) we get
x̂ML (y) = y. (3.226)
Since the mean of y is x, this estimate is unbiased. Furthermore, using the fact
that
λML (x) = var y = x2
156 Estimation Theory Chap. 3

we obtain
(y − x)2
   
∂
Iy (x) = E ln py (y; x) =E
∂x x4
1 1 1
= 4 x2 = 2 = . (3.227)
x x λML (x)
Hence, the Cramér-Rao lower bound is tight and the ML estimate is efficient. Note
that in this case the variance of the estimator and thus the Cramér-Rao bound are
functions of x.

All of our results on nonrandom parameter estimation apply equally well to

the case in which y is discrete-valued, as we illustrate with the following example.

Example 3.15
Suppose we observe a vector
 T
y = y1 y2 · · · yM
of independent Poisson random variables with unknown mean x, i.e., for i = 1, 2, . . . , M
we have
xyi e−x
pyi [yi ; x] = Pr [yi = yi ; x] = . (3.228)
yi !
In this case,
M
X M
X M
X
ln py [y; x] = ln pyi [yi ; x] = (yi ln x − x) − ln(yi !) (3.229)
i=1 i=1 i=1

so that x̂ML (y) is the unique solution to

M
∂ ln py [y; x] X  yi 
= − 1 = 0. (3.230)
∂x x
i=1

In particular, we obtain
M
1 X
x̂ML (y) = yi , (3.231)
M
i=1
which again is then unbiased.
Since the variance of a Poisson random variable equals its mean, we have
M
1 X x
λML = 2
x= . (3.232)
M M
i=1

Using (3.179) with (3.229) we get that the Fisher information is

"M #
1 X M
Iy (x) = 2 E yi = , (3.233)
x x
i=1

so comparing (3.233) with (3.232) we get that the ML estimate is efficient. Further-
more, since λML → 0 as M → ∞, we see that the ML estimate is also consistent.
Sec. 3.3 Nonrandom Parameter Estimation 157

3.3.5 Estimation of Nonrandom Vectors

In this section, we explore some extensions of the preceding results to the problem
of estimating a vector of nonrandom parameters x. To begin, let’s briefly discuss
the extension of the Cramér-Rao bound to this case. In particular, we have that the
covariance matrix Λx̂ (x) of any unbiased estimator satisfies the matrix inequality

Λx̂ (x) ≥ I−1

y (x), (3.234)

where Iy (x) is now the Fisher Information matrix

" T  #
∂ ln py (y; x) ∂ ln py (y; x)
Iy (x) = E
∂x ∂x
 2 
∂ ln py (y; x)
= −E (3.235)
∂x2

Note that from the diagonal elements of (3.234) we obtain a set of scalar
Cramér-Rao bounds on the variances of individual components of x. Also an un-
biased efficient estimate x̂(y) exists if and only if
 T
∂ ln py (y; x)
x̂(y) = x + I−1
y (x) (3.236)
∂x

is a valid estimator, i.e., if and only if the right-hand side of (3.236) does not depend
on x. Also, if an efficient unbiased estimate exists, it is the ML estimate.
To derive the matrix Cramér-Rao bound (3.234), we follow an approach anal-
ogous to that used to obtain (3.167), but which requires some additional steps. In
particular, we begin by recalling that for unbiased estimators the error

e(y) = x̂(y) − x (3.237)

has zero mean, i.e.,

E [e(y)] = 0 (3.238)
and covariance
E e(y)eT (y) = Λx̂ (x).
 
(3.239)
Next we define
∂
f T (y) = ln py (y; x) (3.240)
∂x
and note that using the identity

∂ 1 ∂
ln py (y; x) = py (y; x), (3.241)
∂x py (y; x) ∂x
158 Estimation Theory Chap. 3

we get that f(y) has zero mean:

 
1 ∂
E f T (y) = E
 
py (y; x)
py (y; x) ∂x
+∞
∂
Z
= py (y; x) dy
−∞ ∂x
+∞
∂ ∂
Z
= py (y; x) dy = 1 = 0, (3.242)
∂x −∞ ∂x
and, in turn, covariance
Λf (x) = E f(y)f T (y) = Iy (x).
 
(3.243)
Finally, again using the identity (3.241), the covariance between e(y) and f(y) is
given by

cov (e(y), f(y)) = E e(y)f T (y)

 
Z +∞
∂
= (x̂(y) − x) py (y; x) dy
−∞ ∂x
 Z +∞   Z +∞ 
∂ ∂
= x̂(y) py (y; x) dy − x py (y; x) dy
∂x −∞ ∂x −∞
= I − 0 = I. (3.244)

Next, for an arbitrary choice of c we let

ẽ(y) = cT e(y) (3.245)
and
˜ = cT I−1 (x)f(y) = f T (y) I−1(x) c.
f(y) (3.246)
y y

Then both ẽ(y) and f˜(y) have zero-mean and, using (3.239), (3.243) and (3.244), we
have

var ẽ(y) = cT Λx̂ (x)c (3.247a)

˜ = cT I−1 (x)c
var f(y) (3.247b)
y
 
cov ẽ(y), f˜(y) = c Iy (x)c.
T −1
(3.247c)

Now since the covariance between ẽ(y) and f˜(y) satisfies the bound
h  i2
cov ẽ(y), f˜(y) ˜
≤ var ẽ(y) var f(y) (3.248)

we can substitute (3.247) into (3.248) to obtain, after some simple manipulation,
cT I−1
 T T −1

y (x)c c Λx̂ (x)c − c I y (x)c ≥ 0. (3.249)
Sec. 3.3 Nonrandom Parameter Estimation 159

However, since I−1y (x) is positive semidefinite, the term to the left of the brackets
in (3.249) is non-negative. Hence, the term in brackets must be non-negative. But
then since c is arbitrary this means Λx̂ (x) − I−1
y (x) must be positive semidefinite,
which establishes (3.234) as desired.
Finally, equality is satisfied in (3.248) (and therefore (3.249)) if and only if
˜ for some function k(x) that doesn’t depend on y, i.e., if and only
ẽ(y) = k(x)f(y)
if,
cT e(y) = cT k(x) I−1
y (x) f(y). (3.250)
However, since (3.250) holds for any choice of c we must have
e(y) = k(x) I−1
y (x) f(y). (3.251)

Again k(x) can’t be arbitrary. In particular, when the bound (3.234) is satis-
fied with equality we have
E e(y) eT (y) = Λx̂ (x) = I−1
 
y (x). (3.252)
However, using (3.251), (3.243), and (3.244) we have
E e(y) eT (y) = E e(y) f T (y) I−1 = I−1
   
y (x) k(x) y (x) k(x). (3.253)
Comparing (3.252) with (3.253) we obtain
k(x) = 1, (3.254)
which when substituted into (3.251) yields the following: x̂(y) is an efficient esti-
mator, i.e., satisfies the bound (3.234) with equality if an only if it can be expressed
in the form (3.236) where the right-hand side must be independent of x for the
estimator to be valid.
We can also readily verify that the second form of the Fisher information in
(3.235) is equivalent to the first. Analogous to our approach in the scalar case, we
begin by observing
Z +∞
py (y; x) dy = 1. (3.255)
−∞

Computing the Jacobian of (3.255) with respect to x and using the identity (3.241)
yields
Z +∞  T
∂
py (y; x) ln py (y; x) dy = 0. (3.256)
−∞ ∂x
Finally, computing the Hessian of (3.255) with respect to x and again using (3.241)
we obtain
Z +∞  2 
∂
py (y; x) ln py (y; x) dy
−∞ ∂x2
Z +∞  T  
∂ ∂
+ py (y; x) ln py (y; x) ln py (y; x) dy = 0. (3.257)
−∞ ∂x ∂x
160 Estimation Theory Chap. 3

which verifies that the two expressions in (3.235) are consistent.

It is also straightforward to verify that when an efficient estimator x̂eff (y)
exists, it must be the ML estimator. Again we follow an approach analogous to the
scalar case. Since (3.236) must not be a function of x when an efficient estimator
exists, we can then freely choose any value of x in this expression without effect.
If we choose the value x = x̂eff (y), we obtain
 T
−1 ∂
Iy (x) ln py (y; x) = 0. (3.258)
∂x
x=x̂eff (y)

But since Iy (x) is nonsingular except in the trivial case, we have that the term in
brackets in (3.258) must be zero, i.e.,
x̂eff (y) = x̂ML (y) = arg max py (y; x). (3.259)
x

Again we stress that one should not infer from these results that the ML
estimator is always efficient. When no efficient estimator exists, the ML estimate
can still be computed; however it need not have any special properties. As in
the scalar case, though, even when an efficient estimator doesn’t exist, the ML
estimator often has good asymptotic properties in several problems. One class of
problems in which the ML estimator is always efficient are the linear/Gaussian
problems. We conclude this section with the canonical example.
Example 3.16
Suppose we that our observed data y depends on our parameter vector x through
the linear model
y = Hx + w, (3.260)
where w ∼ N (0, Λw ). In this case
 
1 T −1
py (y; x) = N (y; Hx, Λw ) ∝ exp − (y − Hx) Λw (y − Hx) (3.261)
2
so that maximizing py (y; x) with respect to x is equivalent to minimizing
1
J(x) = (y − Hx)T Λ−1
w (y − Hx) (3.262)
2
with respect to x. Since (3.262) is a non-negative function, its unique stationary
point, which we obtain by setting the Jacobian of (3.262) to zero, is its global mini-
mum and thus gives the ML estimate
x̂ML (y) = (HT Λ−1 −1 T −1
w H) H Λw y (3.263)
This estimate is unbiased, since
E [x̂ML (y)] = (HT Λ−1 −1 T −1
w H) H Λw (Hx + E [w]) = x (3.264)
and its error covariance is
h   T −1 −1 T −1 T i
ΛML = E (HT Λ−1 w H)−1 T −1
H Λw w (H Λw H) H Λw w

= (HT Λ−1 −1 T −1 −1 T −1
w H) H Λw Λw Λw H(H Λw H)
−1

= (HT Λ−1 −1
w H) . (3.265)
 Sec. 3.4 Nonlinear Estimation 161

Note that for this estimate to make sense, HT Λ−1 w H must be invertible, and
this in turn requires that the dimension of y (or, more precisely, the rank of Λw ) be
at least as large as the dimension of x. Phrased differently, the number of degrees of
freedom in the measurements must equal or exceed the number of parameters to be
estimated.
The Fisher information matrix for this problem is obtained using the second
form of (3.235) and yields
d2
Iy (x) = − J(x) = HT Λ−1 w H (3.266)
dx2
which by comparison to (3.265) allows us to conclude that the ML estimate is, in fact,
efficient. Note as well that the estimator covariance (and thus the Fisher matrix) is
independent of x in this example.
As in the scalar case, it is again interesting to compare the ML estimator in this
example to the LLS estimator for the closely related problem developed in Exam-
ple 3.7. In both examples, the measurement models (3.260) and (3.104) are identical,
but in this example x is a nonrandom parameter vector while in Example 3.7 we
have a random parameter x with zero-mean and covariance Λx .
If we added the Gaussian assumptions to Example 3.7, we can again conclude
that the resulting LLS estimator (3.102) is also the BLS estimator and the MAP esti-
mator for the problem. For this reason, we’ll simply use x̂B (y) to denote this estima-
tor and ΛB to denote its error covariance for the remainder of this example. Further-
more, since our ML estimate is efficient, it is the MVU estimator for the nonrandom
parameter estimation problem,11 so we’ll use x̂MVU (y) to denote this estimator and
ΛMVU to denote its covariance for the remainder of this example.
In this case we have, using the alternative matrix forms developed in Ap-
pendix 3.A,

x̂B (y) = mx + ΛB HT Λ−1

w (y − Hmx ) (3.267)
Λ−1
B = Λ−1
x + Λ−1
MVU (3.268)
(3.269)
From these expressions we see that as Λx → ∞ (again in the sense of its trace),
ΛB → ΛMVU (3.270)
and, in turn,
x̂B (y) → x̂MVU (y). (3.271)

3.4 NONLINEAR ESTIMATION

Quite frequently in practice our observations y correspond to some noisy nonlin-

ear function of the parameters x. Let us explore this general problem in the con-
text of the nonrandom parameter estimation theory of Section 3.3. There are lots
of important examples of problems that fall into this category, and we’ll explore

11
This result is referred to as the Gauss-Markov theorem.
162 Estimation Theory Chap. 3

in detail one involving estimation of the parameters of a sinusoid in Section 3.4.1.

Before we do that, however, let us begin with some preliminary observations.

Example 3.17
Consider the following nonlinear measurement
y = h(x) + w , (3.272)
where w ∼ N (0, σ 2 ). In this case,
py (y; x) = N (y; h(x), σ 2 ), (3.273)
so that  
∂ ln py (y; x) y − h(x) dh(x)
= . (3.274)
∂x σ2 dx
Let’s compute the Cramér-Rao bound on the performance of arbitrary unbi-
ased estimates x̂(·) for x. Using (3.274) we obtain that
"  # 
y − h(x) dh(x) 2 dh(x) 2 1 dh(x) 2
      
w 2
Iy (x) = E = E = 2 ,
σ2 dx dx σ2 σ dx
(3.275)
so that
σ2
λx̂ (x) ≥ (3.276)
(dh(x)/dx)2
for any unbiased estimate. Now an efficient estimate exists if and only if (3.188) is a
valid estimator, i.e., if and only if
 
1 ∂ h(x) y
x+ ln py (y; x) = x − + (3.277)
Iy (x) ∂x dh(x)/dx dh(x)/dx
is a function only of y. However, since the right-hand term in (3.277) is the only one
that depends on y and since y can be arbitrary, we can conclude that no efficient
estimate can exist unless dh(x)/dx does not depend on x. However, this will only be
the case when h(·) is a linear (affine) function. Hence, efficient estimates fail to exist
in the strictly nonlinear case.
Consider, for example, h(x) = x3 . In this case, (3.277) becomes
y − x3 2 1 y
x+ = x+ , (3.278)
3x2 3 3 x2
from which we see that there is no efficient estimate.
Since an efficient estimate generally doesn’t exist, the ML estimate, which
we’ll now compute, needn’t have any special properties in the nonlinear case. When
h(·) is invertible, as we’ll assume in this example, we get immediately from (3.274)
that
x̂ML (y) = h−1 (y), (3.279)
where h−1 (·) is the inverse function of h(·), i.e., h−1 (h(x)) = x. Calculating the bias
bML (x) and variance λML (x) of this estimate is difficult in general, though in general
it will be biased. And when biased, this means we cannot even conclude that its
variance λML (x) satisfies (3.276) for even one value of x.
Sec. 3.4 Nonlinear Estimation 163

While ML estimates in nonlinear problems needn’t have good variance char-

acteristics, they do have some attractive features. For example, suppose that a
parameter θ is related to x via
θ = g(x),
where g(·) is a nonlinear but invertible transformation. Then it is a straightforward
exercise to show that the ML estimates are also related by
θ̂ML (y) = g(x̂ML(y)), (3.280)
i.e., ML estimates commute under nonlinear transformations. This is an extremely
convenient property, and one that is not shared by most of the other estimators
we’ve explored in this chapter of the notes. For example, Bayesian estimators
almost never commute with nonlinear transformations, i.e., if x is a random pa-
rameter, then
θ̂B (y) 6= g(x̂B (y)),
for almost any nontrivial cost criterion. We also remark that although (3.280)
doesn’t apply when g(·) is not invertible, straightforward extensions of this result
can be developed to handle the non-invertible case.
We remark, however, that even if x̂ML (y) has nice properties, these properties
are generally not preserved under the transformation (3.280). For example, since
typically h i
E θ̂ML (y) = E [g(x̂ML(y))] 6= g (E [x̂ML (y)])

we wouldn’t expect θ̂ML (y) to be unbiased even if x̂ML (y) were.

Nonlinear estimation problems are distinguished in other important ways
from inherently linear estimation problems. In the remainder of the chapter, we
explore such distinguishing characteristics. To illustrate the main ideas before we
develop them in detail, we first explore as a case study a particular nonlinear es-
timation problem that arises in an extraordinarily wide range of practical applica-
tions.

3.4.1 Sinusoid Estimation

In this section we explore a basic problem involving sinusoid estimation. In par-

ticular, given noisy observations of the form
y [n] = A cos(ω0 n + Θ) + w [n], n = 0, 1, . . . , N − 1 (3.281)
we wish to estimate one or more of the nonrandom parameters A, ω0 , or Θ, where
A > 0 and 0 ≤ ω0 < π. In (3.281), we’ll assume that the noise samples w[n] are
independent, identically distributed N(0, σ 2 ) random variables; this will facilitate
our analysis, and is often a good model in applications.
Among the enormous number of applications in which this model arises are
analog communications, Doppler radar, noise cancellation, interference suppres-
sion, radio astronomy, and sonar direction-finding.
164 Estimation Theory Chap. 3

Analog Communication

Amplitude Modulation (AM) In AM systems, the frequency ω0 is known,

but the amplitude varies with time and carries the information. In such
systems we can generally approximate the amplitude as constant over
the block of N samples, and consider the problem of recovering the am-
plitude for each block as an estimation problem. In such problems, the
phase Θ may be known, but more typically is an unknown parameter
that, while not of interest, must be simultaneously estimated. Note that
if ω0 and Θ are both known, then the resulting estimation problem is
inherently linear, and will exhibit certain associated characteristic be-
havior.
Phase Modulation (PM) In PM systems, it is the phase Θ that carries the
information while the frequency remains essentially fixed (and known).
In such problems, the amplitude is generally distorted by the channel
and while similarly not of interest, must be jointly estimated as well.
Frequency Modulation (FM) In FM systems, the frequency ω0 carries the
information and varies with time accordingly. In such systems, we then
wish to estimate ω0 , which is modeled as essentially constant over the
block of length N. Typically, the communication channel distorts both
the amplitude and phase, so these quantities must be simultaneously
(i.e., jointly) estimated.

Noise and Interference Cancellation A wide variety of noise and interference en-
countered in practice is inherently sinusoidal in nature. Examples include
60 Hz (line-frequency) interference in systems due to AC power supplies,
noise from rotating machinery, propeller noise in aircraft and on ships, and
narrowband jamming—hostile or inadvertent—in wireless communication
systems. In such cases, the sinusoidal term in (3.281) may be the unwanted
interference and w [n] may represent the (broadband) signal of interest. For
these scenarios, an effective interference suppression strategy involves es-
timating the parameters of the sinusoidal interferer, then subtracting it out
from the observations to recover the signal of interest.

Doppler Radar In radar systems, (3.281) can be used to model the radar return,
where the deviation of ω0 from some nominal value is a Doppler shift used
to measure the velocity of the target.

Radio Astronomy In radio astronomy applications, one is often interested in de-

tecting and locating spectral lines corresponding to emissions from distant
sources of radiation, and even most experimental apparatus designed for
searching for extraterrestrials uses (3.281) as the basic model for the signal
being sent by ET!
Sec. 3.4 Nonlinear Estimation 165

wavefronts
far field
source

sensor array

d φ …
0 1 2 3 N-1

d cos φ

Figure 3.3. Estimating the direction of arrival of a far field source using an
N -element linear array of sensors.

Sonar Direction-Finding Sonar systems are often used to locate the direction from
which an acoustic source is propagating. To illustrate this, suppose the source
is emitting a pure tone (sinusoid) of the form

x(t) = A cos(Ω0 t), (3.282)

and that this signal is being picked up at a linear, horizontal array of N sen-
sors (hydrophones), as depicted in Fig. 3.3. Let d denote the distance be-
tween sensors, and let us assume that the source is sufficiently distant to
allow a so-called “far-field” approximation: the signal arrives at the array as
a plane-wave, with the wavefronts consisting of straight lines (rather than
circles) as Fig. 3.3 reflects. Let φ denote the angle at which the plane wave
impinges on the array.

In this system, the propagation time to the nth sensor is

d
tn = t0 − n cos φ, n = 0, 1, . . . , N − 1 (3.283)
c

where c is the propagation speed (i.e., phase velocity), so that the signal ob-
served at the nth sensor is, for some Θ0 ,

yn (t) = A cos(Ω0 (t − tn ) + Θ0 ) + wn (t). (3.284)

166 Estimation Theory Chap. 3

If we take a snapshot of all the sensors at a particular time instant t = t∗ , we

obtain the vector of samples

y [n] = yn (t∗ )
  
d
= A cos Ω0 cos φ n + Θ + wn (t∗ )
c
= A cos(ω0 n + Θ) + w [n], (3.285)
where Θ = Θ0 + Ω0 (t∗ − t0 ), ω0 = (Ω0 d/c) cos φ, and w [n] = wn (t∗ ). Hence, by
estimating the spatial frequency ω0 , we can indirectly obtain an estimate of
the direction-of-arrival φ.

Performance Issues and Cramér-Rao Bounds

Let us begin by considering the most general problem, wherein the parameters
A, ω0 , and Θ in the model (3.281) are all unknown, and explore the form of the
associated Cramér-Rao bounds. These bounds will give us some insight into how
we can expect estimator performance to vary with the signal-to-noise ratio
1
γ = A2 /σ 2 , (3.286)
2
the data length N, and the actual values of the parameters A, ω0 , and Θ.
When we collect the unknown parameters into a vector x, i.e.,
 
A
x = ω0  ,
 (3.287)
Θ
and do the same for the data, i.e.,
 
y[0]
 y[1] 
y= , (3.288)
 
..
 . 
y[N − 1]
the elements of the Fisher information matrix then take the form
∂2
 
[Iy (x)]ij = −E `(y; x) , (3.289)
∂xi ∂xj
where
N −1
N 2 1 X 2
`(y; x) = ln py (y; x) = − ln(2πσ ) − 2 y[n] − A cos(ω0 n + Θ) . (3.290)
2 2σ n=0

The calculation of the quantities (3.289) is straightforward but somewhat

lengthy; the details are provided in Appendix 3.B. For arbitrary N, the results
Sec. 3.4 Nonlinear Estimation 167

are also somewhat cumbersome. However, as the Appendix shows, in the large
N regime, corresponding to at least moderately sized data sets, the Fisher infor-
mation can be expressed using order notation12 in the following comparatively
simple form
o(N 2 )
 
N/2 + o(N) o(N)
1
Iy (x) = 2  o(N 2 ) A2 TN /2 + o(N 3 ) A2 SN /2 + o(N 2 ) , (3.291)
σ 2 2 2
o(N) A SN /2 + o(N ) A N/2 + o(N)
where
N −1
X 1
SN = n = N(N − 1) (3.292)
n=0
2
N −1
X 1
TN = n2 = N(N − 1)(2N − 1). (3.293)
n=0
6

Computing the inverse of (3.291) we then obtain our Cramér-Rao bound on

unbiased estimates Â(y), ω̂0 (y), and Θ̂(y) of the parameters in the large N regime.
In particular, we obtain, using (3.286),
!
2σ 2
 
Â(y) 1  −1  1 1
var ≥ 2 Iy (x) 11 ≈ 2 = ∼O (3.294a)
A A AN γN γN
 
 −1  12 1
var ω̂0 (y) ≥ Iy (x) 22 ≈ ∼O (3.294b)
γN(N 2 − 1) γN 3
 
 −1  2(2N − 1) 1
var Θ̂(y) ≥ Iy (x) 33 ≈ ∼O (3.294c)
γN(N + 1) γN

It should be emphasized that the size of N necessary for the approximations

in (3.294) to be valid depends on the true value of ω0 —in general, the closer ω0 is to
0 or π, the larger N must be. More specifically, as the development in Appendix 3.B
reveals, the approximation is valid provided
 
π 1
 ω0  π 1 − . (3.295)
N N

12
For functions f (·) and g(·) we use the notation f (N ) ∼ o(g(N )) to indicate that f (N ) grows
strictly slower than g(N ), i.e.,
f (N )
lim = 0,
N →∞ g(N )

As related order notation, we write f (N ) ∼ O(g(N )) if f (N ) grows no faster than g(N ), i.e.,
f (N )
lim < ∞.
N →∞ g(N )

As examples, f (N ) ∼ o(N ) means that f (N ) grows slower than linearly with N , while f (N ) ∼
O(N ) means that f (N ) grows no faster than linearly with N .
168 Estimation Theory Chap. 3

The asymptotic Cramér-Rao lower bounds (3.294) reveal some key charac-
teristics of the estimation problem. As we would expect, all the bounds decrease
inversely with the SNR γ and the data length N. However, data length has the
most profound impact on the bound for the frequency estimate. This suggests that
it may be possible to estimate this parameter with very high accuracy at moderate
data lengths. Whether this is possible depends, of course, on whether estimators
can be developed whose performance comes close to the bound. We explore this
issue, among others, in the context of developing asymptotic ML estimates for the
parameters in the next section.
As a final remark, it is worth emphasizing that the Fisher information (3.291)
contains all the information necessary to asymptotically bound the performance of
related sinusoid estimation problems. In particular, when some of the parameters
A, ω0 , Θ are known, the associated Cramér-Rao bounds for the remaining parame-
ters are obtained by inverting a submatrix of (3.291) formed by discarding the rows
and columns corresponding to the known parameters. Using this approach, it can
be readily verified that, for example, when the frequency ω0 is known, the asymp-
totic Cramér-Rao bounds for  and Θ̂ are still O(1/γN) as in (3.294a) and (3.294c),
respectively. Likewise, when both the frequency ω0 and phase Θ are known, the
Cramér-Rao bound on  remains O(1/γN) as in (3.294a).

Maximum Likelihood Estimates

In this section, we obtain ML estimates for the sinusoid estimation problem, de-
velop their properties, and relate the performance of resulting estimators to the
corresponding Cramér-Rao bounds.
To begin, the ML parameter estimates
x̂(y) = arg max `(y; x)
x

are the solutions to the nonlinear least-squares problem

(Â, ω̂0 , Θ̂) = arg min J(A, ω0 , Θ) (3.296a)
(A,ω0 ,Θ)

where, via (3.290),

N −1
X  2
J(A, ω0 , Θ) = y[n] − A cos(ω0 n + Θ) . (3.296b)
n=0

In principle, the optimization (3.296) can be performed numerically without

exploiting any of the special structure in the problem. However, in the large N
regime, more direct expressions are possible, which have intuitively satisfying in-
terpretations and lead to efficient implementations. These expressions are in terms
of the normalized, length-N discrete-time Fourier transform of the data segment
Sec. 3.4 Nonlinear Estimation 169

comprising y, i.e.,13
N −1
jω 1 X
YN (e ) = √ y[n]e−jωn . (3.299)
N n=0
The magnitude-squared of (3.299), i.e., |YN (ejω )|2 , is referred to as the periodogram
of the data.
A natural periodogram-based estimator for the sinusoid estimation problem
is defined as follows.

Definition 3.1 (Periodogram-Based Estimator) The periodogram-based frequency es-

timate ω̂0 is the location of the peak of the periodogram, i.e.,
2
ω̂0 = arg max YN (ejω ) . (3.300)
ω

In turn, the magnitude of this peak yields the associated amplitude estimate, i.e.,
4 2
Â2 = YN (ej ω̂0 ) , (3.301)
N
and the associated phase estimate corresponds to the (negated) phase of YN (ejω ) at the
location of the peak, i.e.,
!
j ω̂0

Im Y N (e )
Θ̂ = −]YN (ej ω̂0 ) = − tan−1 (3.302)
Re {YN (ej ω̂0 )}

Note that the estimator in Definition 3.1 is both intuitively appealing and
highly practical. Indeed, to identify a sinusoid it is rather natural to compute the
Fourier transform of the noisy data segment and locate the amplitude, frequency,
and phase of its peak. Moreover, these estimators can be implemented very ef-
ficiently in practice. In particular, the frequency estimate can be computed by
taking a sufficiently large discrete Fourier transform (DFT) of the data—i.e., with
sufficient zero-padding of the data—and searching for index of the largest DFT co-
efficient. The computation of these DFT’s can be conveniently carried out using an
efficient fast Fourier transform (FFT) algorithm, which has O(N log N) complexity.
The estimator of Definition 3.1 also has some important optimality proper-
ties, and in fact is closely related to the ML estimator for the problem. In particular,

13
It is often convenient to view (3.299) as the Fourier transform of a windowed version of
the sequence y[n], i.e.,
YN (ejω ) = F {gN [n] y[n]} (3.297)
where g[n] is the unit-energy window
( √
1/ N 0≤n≤N −1
gN [n] = . (3.298)
0 otherwise
170 Estimation Theory Chap. 3

the periodogram-based and ML estimators are effectively equivalent when N is

large in the sense of (3.295); from this perspective, we can view the periodogram-
based estimator as the asymptotic ML estimator. A relatively straightforward
derivation of this result is developed in Appendix 3.C.
In order to analyze the performance of the periodogram-based estimator, it
is useful to decompose the periodogram into signal and noise components. In
particular, we write YN (ejω ) in the form
YN (ejω ) = XN (ejω ) + WN (ejω ) (3.303)
where
N −1
jω 1 X
WN (e ) = √ w [n] e−jωn (3.304)
N n=0
and
N −1
jω 1 X
XN (e ) = √ x[n] e−jωn , (3.305)
N n=0
with
x[n] = A cos(ω0 n + Θ). (3.306)
The noise component WN (ejω ) has mean
N −1
jω 1 X
E [w [n]] e−jωn = 0
 
E WN (e ) = √ (3.307)
N n=0
and variance
N −1 N −1
h
jω 2
i 1 XX
E WN (e ) = E [w [n] w [m]] ejω(n−m) = σ 2 . (3.308)
N n=0 m=0

With this decomposition, the “signal” and “noise” components of the peri-
odogram are naturally defined as, respectively,
2 2
E YN (ejω ) = XN (ejω )

(3.309)
and h i
2
var YN (ejω ) = E WN (ejω ) = σ2 . (3.310)
These two components are depicted in Fig. 3.4. In turn, we define the SNR at a
particular frequency as
2 2
|E [YN (ejω )]| |XN (ejω )|
γ(ω) = = . (3.311)
var YN (ejω ) σ2
In the large N regime (3.295) we have
( 2  2 )
2
2 A sin(ω − ω 0 )N/2 sin(ω + ω 0 )N/2
XN (ejω ) ≈ + . (3.312)
4N sin(ω − ω0 )/2 sin(ω + ω0 )/2
Sec. 3.4 Nonlinear Estimation 171

ω0 =π/2 N=32
SNR: −12 dB
0
Fourier Transform Power (dB)
−5 SNR: −6 dB

−10
SNR: 0 dB

−15
SNR: 6 dB Figure 3.4. Signal and noise com-
−20 ponents of the periodogram, when
ω0 = π/2 and N = 32. The
−25
solid curve depicts the signal compo-
nent, 10 log10 (|XN (ejω )|2 ), while the
dash lines depict the
 noise components
−30

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 σ 2 = 10 log10 (E |WN (ejω )|2 ) corre-
normalized frequency ω/2π sponding to various values of SNR.

and thus the SNR (3.311) effectively attains its peak at ω = ω0 . Since
2 A2 N
XN (ejω0 ) ≈ ,
4
the peak SNR is, using (3.286),
A2 N 1
2
= γN. γ(ω0 ) = max γ(ω) =
(3.313)
ω 4σ 2
Note that (3.313) is a factor of N/2 larger than γ, the SNR for the original data.
The performance characteristics of the periodogram-based parameter esti-
mators have some special features. To illustrate this, the variance in the amplitude
and frequency estimates are plotted as a function of SNR in Figs. 3.5 and 3.6, re-
spectively, along with the associated Cramér-Rao bounds. These figures reveal a
distinct threshold phenomenon: for a given data length, there exists a SNR thresh-
old above which the estimator variance closely tracks the Cramér-Rao bound, and
below which the estimator variance diverges sharply from the bound. The phe-
nomenon is particularly pronounced for the frequency estimator, but arises with
the amplitude estimator as well.
This threshold behavior, which is also referred to as the “capture” effect, can
be understood as follows. When the SNR is high enough that the peak in the
periodogram at the true frequency protrudes prominently above the noise, the
peak can be located quite accurately and the parameter estimation errors are due
to slight, noise-induced distortion of the true peak. This is the regime in which the
estimator performance tracks the Cramér-Rao bound. On the other hand, when
the SNR is low enough that the correct peak lies below the noise and is obscured
by other peaks, catastrophic estimation errors due to the estimator selecting the
wrong peak entirely, leading to anomalous parameter estimates. This is the regime
in which the estimator performance diverges from the Cramér-Rao bound. Sample
periodograms corresponding to the different regimes are depicted in Fig. 3.7.
172 Estimation Theory Chap. 3

2
10

1
10

ω0 = π/2
0
10
var A/A
ˆ

−1
10

N=16
−2
10
N=32
Figure 3.5. Variance of the ML esti-
N=64
10
−3
mate of the amplitude of an unknown
sinusoid in white Gaussian noise as
−4
a function of SNR γ for various data
10
−20 −15 −10 −5 0 5 10 15 20 lengths N . The dashed lines are the as-
SNR γ (dB) sociated Cramér-Rao bounds.

−1
10

−2
10

ω0 = π/2
−3
10

−4
ˆ0

10
var ω

−5
10
N=16

−6
10
N=32
Figure 3.6. Variance of the ML esti-
−7 mate of the frequency of a unknown
10 N=64
sinusoid in white Gaussian noise as
−8
a function of SNR γ for various data
10
−20 −15 −10 −5 0 5 10 15 20 lengths N . The dashed lines are the as-
SNR γ (dB) sociated Cramér-Rao bounds.
 Sec. 3.4 Nonlinear Estimation 173

20

15 SNR: 6 dB ω0 =π/2 N=64

10
Periodogram (dB)

−5

−10

−15

−20
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
normalized frequency ω/π

20

15 SNR: −1 dB ω0 =π/2 N=64

10
Periodogram (dB)

−5

−10

−15

−20
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
normalized frequency ω/π

20

15 SNR: −6 dB ω0 =π/2 N=64

10
Periodogram (dB)

Figure 3.7. Sample periodograms for

−5
a sinusoid of frequency ω0 = π/2 and
−10
length N = 64. The top figure corre-
sponds to SNR of 6 dB, which is above
−15 threshold. The middle figure corre-
sponds to the approximate threshold
−20
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
SNR of -1 dB. The bottom figure corre-
normalized frequency ω/π sponds to a sub-threshold SNR of -6 dB.
 174 Estimation Theory Chap. 3

1
10

ω0 = π/2
0
10

SNR γ = −9 dB
−1
10
var A/A
ˆ

SNR γ = −3 dB
−2
10

SNR γ = 3 dB

10
−3 Figure 3.8. Variance of the ML esti-
mate of the amplitude of an unknown
sinusoid in white Gaussian noise as a
10
−4 function of data length N for various
10
1
10
2 3
10 SNRs γ. The dashed lines are the as-
Data length N sociated Cramér-Rao bounds.

0
10

SNR γ = −9 dB
ω0 = π/2
−2
10

−4 SNR γ = −3 dB
10
ˆ0
var ω

−6
10

SNR γ = 3 dB
−8
10

Figure 3.9. Variance of the ML esti-

−10
10 mate of the frequency of a unknown
sinusoid in white Gaussian noise as a
−12
10
function of data length N for various
1
10 10
2
10
3
SNRs γ. The dashed lines are the as-
Data length N sociated Cramér-Rao bounds.

Among other features revealed by Figs. 3.5 and 3.6, we see that since the
peak SNR (3.313) is proportional to data length N, the threshold SNR decreases
as N increases. Also, the slope of the bounds in the two figures are the same, re-
flecting the same inverse dependence on SNR γ [cf. (3.294a) and (3.294b)], though
the offsets are quite different due to the different nature of the dependence on N.
Performance variations with block length N are more fully apparent in Figs. 3.8
and 3.9. These figures show the variance in the estimates of A and ω0 , respectively,
plotted as a function of N for several values of the SNR γ. Note that the slope
of the Cramér-Rao bounds is greater by a factor of 3 (on the log-log scale) for the
frequency estimate. This is because of the 1/N 3 vs. 1/N dependence in the bounds
apparent in comparisons of (3.294b) and (3.294a).
It is important to emphasize that, by contrast, linear estimation problems do
not exhibit the kind of threshold behavior observed above. In fact, for linear esti-
mation problems involving Gaussian data, we established that ML estimates are
Sec. 3.4 Nonlinear Estimation 175

efficient, so the associated Cramér-Rao bounds are accurate predictors of the per-
formance attainable in practice. This is the case in sinusoid estimation problems
where only the amplitude is unknown.
These distinctions underlie the familiar differences in the way signal quality
varies in, e.g., AM and FM radio reception. AM reception has the characteristic
that the quality degrades steadily with increasing distance from the source of the
transmission. On the other hand, FM systems have the characteristic that within
a certain radius of the source the quality of the reception is higher than corre-
sponding AM systems, but that that outside this service area reception deteriorates
sharply as the SNR drops below threshold.
More generally, the capture effect is a dominant feature of systems in many
applications where there are inherent nonlinearities. In the next section, we dis-
cuss how the effect arises in this more general setting, and view the sinusoid esti-
mation problem as a special instance of the phenomenon.

3.4.2 Threshold Behavior and the Capture Phenomenon

In this section, let’s consider a vector generalization of Example 3.17 in which the
measurements y depend on the parameter vector x via
y = h(x) + w (3.314)
with w ∼ N(0, Λw ), so that the measurements take the form of Gaussian random
vector.
Let us first determine the associated Cramér-Rao bound for the problem. To
begin, first note that, provided Λw > 0 and h(·) is differentiable,
M 1 1
ln py (y; x) = − ln(2π) − ln |Λw | − (y − h(x))T Λ−1
w (y − h(x)), (3.315)
2 2 2
so
∂ dh(x)
ln py (y; x) = (y − h(x))T Λ−1
w . (3.316)
∂x dx
In turn, using (3.316) in (3.235) we obtain the Fisher matrix
dh(x)T −1
 
T −1 dh(x)
Iy (x) = E Λw (y − h(x))(y − h(x)) Λw
dx dx
T
dh(x) −1  dh(x)
Λw E wwT Λ−1

= w
dx dx
T
dh(x) −1 dh(x)
= Λw . (3.317)
dx dx
From (3.317) we see that Iy (x) > 0 if and only if the Jacobian matrix14 dh(x)/dx
is nonsingular. In this case, the Cramér-Rao bound on the covariance of unbiased

14
This matrix was defined in Appendix 1.B of Chapter 1.
176 Estimation Theory Chap. 3

estimates x̂(y) of x is
−1
dh(x)T −1 dh(x)

Λx̂ (x) ≥ Λw . (3.318)
dx dx

As in Example 3.17, it is straightforward to show that an efficient estimate

fails to exist when h(·) is strictly nonlinear, i.e., unless h(·) is a linear (affine) func-
tion. Nevertheless, the Cramér-Rao bound (3.329) does have the useful interpreta-
tion as the performance of a closely related linear system. To see this, consider a
linearization of h(·) about a particular value x∗ that is near x, and let us use x̃ to
denote the deviation, i.e.,
x̃ = x − x∗ . (3.319)
This linearization is obtained from a Taylor series expansion of h(x) about x∗ ,
which, as described in Appendix 1.B of Chapter 1, takes the form
 
dh(x)
h(x) = h(x∗ ) + x̃ + · · · (3.320)
dx x=x∗

In particular, when the higher-order terms (· · · ) are neglected, the measurement

deviations
ỹ = y − h(x∗ ) (3.321)
are related to the parameter deviation x̃ according to the linear model
 
dh(x)
ỹ ≈ x̃ + w. (3.322)
dx x=x∗

To verify this it suffices to substitute (3.320) into (3.321). This linearization is de-
picted in Fig. 3.10 for the case in which both x and y are scalars.
When (3.322) holds with equality, the minimum variance unbiased estimate
of x̃ based on ỹ is the ML estimate, as is that for x based on y. In particular, via the
Gauss-Markov theorem (Example 3.16) we have
dh(x)T
x̂MVU (y) = x∗ + Φ−1 (x∗ ) Λ−1
w (y − h(x∗ )) (3.323)
dx x=x∗

where
dh(x)T −1 dh(x)
Φ(x) = Λw . (3.324)
dx dx
Moreover, via the Gauss-Markov theorem we also have that the covariance of this
estimate is
ΛMVU = Φ−1 (x∗ ), (3.325)
which when x∗ = x corresponds to the Cramér-Rao bound for the problem [cf.
(3.318)].
This analysis implies that if we know a priori that x lies in a neighborhood of
x∗ , and that the neighborhood is small enough that the Jacobian matrix dh(x)/dx
Sec. 3.4 Nonlinear Estimation 177

h(x* ) + dh (x-x*)
dx x=x*

h(x)

x Figure 3.10. Linearization of the mea-

x* surement function h(·) in the Gaussian
nonlinear estimation problem.

can be well approximated as essentially constant in that neighborhood, then the

Cramér-Rao bound will yield an accurate measure of achievable performance. For
this reason, the Cramér-Rao bound is frequently referred to as a local bound.
When, as is the case more generally, such a priori localization information
is not available, then ML estimates x̂ML (y) are generally not efficient. However,
these estimates are often asymptotically efficient at high SNR, where estimation
errors are small. To understand this property, we explore the ML estimate and its
relationship to the Cramér-Rao bound in more depth.
To begin, the ML estimate is obtained by maximizing (3.315) or, equivalently,
as the solution to the following nonlinear least-squares minimization problem
x̂ML (y) = arg min J(a) (3.326a)
a

where
J(a) = (y − h(a))T Λ−1
w (y − h(a)). (3.326b)
Because the resulting estimate will generally be biased, the bound (3.318) won’t
apply, but again may asymptotically.
Although the ML estimate may lack specific optimality properties, it is an at
least intuitively reasonable estimate, as the form of (3.326b) reveals. In particular,
the maximum likelihood estimate minimizes a weighted sum of squared errors,
where the weighting is determined by Λ−1 w so that more accurate observations are
weighted more heavily. For this reason, this type of “least-squares” estimate is
often used without any justification in terms of ML estimation.
In the sequel, it will be convenient to express the ML estimate in a different
form. In particular, expanding out the quadratic form (3.326b) and discarding
the quadratic term in y (since it doesn’t depend on a and hence won’t affect our
178 Estimation Theory Chap. 3

optimization over a) we see that the ML estimate can be rewritten as

x̂ML (y) = arg max r(y; a), (3.327a)
a

where
1 T
r(y; a) = hT (a)Λ−1 −1
w y − h (a)Λw h(a). (3.327b)
2

To simplify our exposition, in the sequel let us restrict our attention to the
case in which x is an unknown scalar x. In this case, the measurement model
(3.314) and Cramér-Rao bound (3.318) specialize to, respectively,
y = h(x) + w (3.328)
and −1
dhT (x) −1 dh(x)

Λx̂ (x) ≥ Λw . (3.329)
dx dx
Also, we rewrite (3.327b) in this case as
1 T
r(y; a) = hT (a)Λ−1 −1
w y − h (a)Λw h(a) (3.330)
2
and view the first term on the right-hand side of (3.330) as a (weighted) inner
product of the observed value y and its candidate values h(a). The second term in
(3.330) is an energy (norm) term, i.e., it is a measure of the corresponding SNR we
expect to see. As we’ll see, in the linear case the dominant effect a has is on signal
energy; this is the case, for example, when we are estimating the amplitude of a
sinusoid of known frequency (the AM problem). On the other hand, in many non-
linear problems a has little or no effect on SNR; this is the case, for example, when
we are estimating the frequency of a sinusoid (the FM problem). It’s these fun-
damental differences that generally leads to the Cramér-Rao bound being overly
optimistic in the nonlinear problem.
The statistics of the objective function (3.330), and hence overall system per-
formance, can be expressed completely in terms of “nonlinear inner products” of
the form
C(x1 , x2 ) = h(x1 )T Λ−1
w h(x2 ) (3.331)
which in essence measure how similar the observations y will be on average if a
is x1 versus x2 . In particular, since r(y; a) is a linear function of y, it is a Gaussian
random variable, and thus is fully described by its mean
1
mr (a) = E [r(y; a)] = C(a, x) − C(a, a) (3.332)
2
and variance
λr (a) = var r (y; a) = C(a, a) = hT (a)Λ−1
w h(a). (3.333)
Fig. 3.11 illustrates a typical example of what these statistics look like as a function
of a. This example corresponds to a special case of the sinusoid estimation problem
Sec. 3.4 Nonlinear Estimation 179

0 Figure 3.11. Statistics of an objective

function r(y; a). This example is gener-
ated from the sinusoid estimation prob-
lem involving N = 16 data samples
corrupted by white Gaussian noise,
−5 with the unknown parameter x be-
ing the normalized frequency of the si-
nusoid (1/2 in this case). The solid
curve depicts the mean mr (a), while
the dashed curve depicts a quadratic
function with the same curvaturep at its
−10
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 peak. The dotted curve depicts λr (a)
a corresponding to an SNR of 12 dB.

of Section 3.4.1 in which only the frequency parameter is unknown (and is denoted
using x). As this figure reflects, on average the peak of the objective function lies
at the true parameter value. Noise in y perturbs the values of r(y; x) away from
the solid curve in the figure, which in turn leads to the peak value shifting and,
hence, estimation error. As the figure also reflects, the standard deviation of the
associated perturbations does not vary strongly with the independent variable a.
A quadratic with the same curvature as mr (a) at its peak is also depicted in
Fig. 3.11. The curvature of mr (a) at its peak (i.e., a = x) is intimately related to
the Cramér-Rao bound for the problem. To see this, note that using (3.315) and
(3.327b) we can express r(y; a) as

r(y; a) = ln py (y; a) + f (y) (3.334)

where f (y) does not depend on a. Thus, the curvature is given by

d2 ∂2
 
mr (a) =E r(y; a)
da2 a=x ∂a2
a=x
 2 
∂
=E ln py (y; a)
∂a2
a=x
= −Iy (x). (3.335)

More generally, for the case of vector parameters x, it is straightforward to verify

that the Hessian matrix (as defined in Appendix 1.B) for mr (a) at its peak a = x is
the negative of the Fisher information, i.e.,

d2
mr (a) = −Iy (x). (3.336)
da2 a=x
180 Estimation Theory Chap. 3

Figure 3.12. Statistics of an objective

function r(y; a). This example is gener-
ated from the sinusoid estimation prob-
−5 lem involving N = 16 data samples
corrupted by white Gaussian noise,
with the unknown parameter x being
the amplitude of the sinusoid (unity in
this case). The solid curve depicts the
mean mr (a), while the dotted curve de-
−10
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 picts λr (a) corresponding to an SNR of
a 12 dB.

Note that when h(·) is a linear function, e.g., h(x) = cx, then mr (a) is quadratic
in a:
1 T −1
2
mr (a) = cT Λ−1


w c ax − c Λw c a . (3.337)
2
This is illustrated in Fig. 3.12. Since the Cramér-Rao bound is tight in this linear
case, the curvature of mr (a) at its peak fully characterizes the performance of the
ML estimator.
From this perspective, in the more general nonlinear case the Cramér-Rao
bound corresponds to is fitting a quadratic at the peak of mr (a) as depicted in
Fig. 3.11, and using the curvature as a measure of the performance. However,
the dashed line in Fig. 3.11 falls off sharply as a function of x, consistent with the
fact that in the linear case signal energy is a strong function of x. For this rea-
son, making very large errors in a linear problem is extremely unlikely because of
the enormous noise energy needed to cause such an error. However, in nonlinear
problems where signal energy is at most a weak function of x, behavior like that
depicted in Fig. 3.11 is more typical. In this case much smaller noise values are
needed to push a value of r(y; a) located far from a = x above the value of r(y; x),
and therefore the Cramér-Rao bound tends to grossly underestimate the probabil-
ity of large estimation errors. Whether this is ultimately significant or not depends
upon the size of the noise variance.
This type of behavior manifests itself as the capture phenomenon we saw
with the sinusoid estimation problem explored in Section 3.4.1. For small noise
variances, the Cramér-Rao bound is accurate, since large errors occur with neg-
ligible frequency. As the noise increases, however, a threshold effect occurs at
some value of the noise variance beyond which large errors become significant.
At this point achievable performance becomes considerably worse than the opti-
mistic Cramér-Rao bound prediction.
Sec. 3.4 Nonlinear Estimation 181

This interpretation provide additional insight into the qualitative differences

in behavior between AM and FM radio reception discussed at the end of Section
3.4.1. For FM systems, where the parameter of interest is frequency, the associated
mr (a) is as depicted in Fig. 3.11. For AM systems, where the parameter of interest
is amplitude and appears linearly in the observations, the associated mr (a) is a
quadratic as depicted in Fig. 3.12. As a comparison between these figures reflects,
the peak is much broader in the AM case than the FM case. Since the breadth of
the peak is determined by the curvature, this implies that at reasonably high SNR,
errors in FM demodulation are much smaller than those in AM demodulation.
As the SNR degrades, however, the AM reception degrades in a rather uniform
manner with SNR, while FM reception degrades abruptly when the SNR reaches
the threshold value.

3.4.3 Computation of ML Estimates

We close this chapter with a discussion of issues associated with the computation
of the ML estimate. As (3.327) reflects, x̂ML is in principle determined by evaluat-
ing r(y; a) for all values of a and choosing that for which the maximum is attained.
Obviously this isn’t viable in practice. One practical approach that can be used is
the following successive-linearization strategy:

1. Make an initial guess of the estimate, i.e., let i = 0 and let x̂(i) = x∗ for some
x∗ .
2. Linearize h(x) about x = x̂(i) , i.e., assume x = x̂(i) + ∆i and solve the lin-
earized estimation problem for ∆ ˆ i.
ˆ i , and increment i.
3. Generate a new estimate via x̂(i+1) = x̂(i) + ∆
4. Go to step 2.

This iterative algorithm typically converges to a local maximum of the objective

function r(y; x). Alternative methods for finding such local maxima can also be
used.
To find the global maximum, a two-stage procedure can be used in princi-
ple. First a coarse search is performed essentially as an M-ary hypothesis test. In
particular, we choose a discrete set of values x1 < x2 < · · · < xM and let
x̂(0) = xı̂ where ı̂ = arg max r(y; xi ). (3.338)
i

The estimate (3.338) can then be used as the initial estimate for a local search based
on, e.g., successive-linearization.
Note that for the coarse search to be useful, the grid points x1 , x2 , . . . , xM
need to be chosen so that r(y; xi ) is likely to be larger than r(y; xj ) for j 6= i if
the actual value of x is closest to xi . For example, in a scenario like that depicted
 182 Estimation Theory Chap. 3

in Fig. 3.11, we might partition the a-axis up into small intervals of width on the
order of the width of the main lobe of the solid curve and take as the xi the centers
of these intervals.
Finally, the two-stage algorithm leads to one more convenient interpretation
of the capture effect in a nonlinear estimation problem. Specifically, when we
choose the spacing between the xi to be sufficiently small that the E [r(y; xi)] are
roughly quadratic near xi the Cramér-Rao bound provides an accurate measure
of estimation error provided the correct xi is chosen in the first, coarse estimation
stage. As the noise variance increases, however, there is an increasing probabil-
ity of making an error in this first stage, and it is this behavior that leads to the
threshold phenomenon.

3.A ALTERNATE FORMULAS FOR LINEAR LEAST-SQUARES ESTIMATION

As mentioned in Section 3.2.5 there are a number of alternate expressions for the
quantities involved in linear-least squares estimation. Recall that the problem is
that of estimating x given the measurements
y = Hx + w (3.339)

For this problem, the LLS estimator (which also corresponds to the BLS es-
timator when x and w are independent Gaussian random vectors) and its perfor-
mance are given by, respectively, (3.108) and (3.109) with (3.110), which we repeat
here for convenience:
x̂LLS (y) = mx + K (y − Hmx ) (3.340)
ΛLLS = Λx − K HΛx HT + Λw KT ,


(3.341)
where
K = Λx HT (HΛx HT + Λw )−1 . (3.342)

Let us first derive the following alternative form for the error covariance
−1
ΛLLS = Λ−1 T −1

x + H Λ w H . (3.343)

Showing (3.343), i.e.,

ΛLLS Λ−1 T −1
 
x + H Λw H = I

is, using (3.341) with (3.342), equivalent to showing that

h −1 i
Λx − Λx HT HΛx HT + Λw HΛx [Λ−1 T −1

x + H Λw H] − I = 0. (3.344)

Eq. (3.344) can in fact be obtained from the expressions (1.231)–(1.235) in Ap-
pendix 1.A to Chapter 1 for inverting block matrices. Here, however, we verify
Sec. 3.A Alternate Formulas for Linear Least-Squares Estimation 183

this more directly. Expanding the expressions in (3.344) and rearranging terms we
find that the left-hand side of (3.344) is equivalent to
  −1 −1 
Λx HT − HΛx HT + Λw Λw + I − HΛx HT + Λw HΛx HT Λ−1

w H
 −1  
= Λx HT I − HΛx HT + Λw HΛx HT + Λw Λ−1

w H

= 0 (3.345)

so that (3.343) is verified.

Let us consider one additional expression for the error covariance. In prac-
tice, (3.343) is typically not well-suited for actual numerical computation of the
error covariance ΛLLS . Neither, however, is (3.341). Specifically, as a covariance
matrix ΛLLS is at least positive semidefinite. Eq. (3.341) expresses ΛLLS as the dif-
ference between two positive semidefinite matrices, and, in cases in which this dif-
ference involves subtracting large numbers, it is possible that numerical errors can
lead to the computed value of ΛLLS losing its definiteness. A detailed investigation
of numerical computation issues is beyond the scope of this course. However, we
point out that rewriting the error covariance in the form

ΛLLS = [I − KH] Λx [I − KH]T + KΛw KT , (3.346)

which involves the sum of positive definite matrices, is much preferred for numer-
ical computation. To derive (3.346), we use (3.339), (3.340) to write

e(x, y) = x − x̂(y) = [I − KH] (x − mx ) − Kw (3.347)

Then, since x and w are uncorrelated, we immediately obtain the expression (3.346)
for ΛLLS .

As a final comment, we note that the gain (3.342) can also be written in an
alternative form, viz.,
K = ΛLLS HT Λ−1
w . (3.348)

To verify this, we use (3.341) to write

h −1 i
ΛLLS HT Λ−1 Λx − Λx HT HΛx HT + Λw HΛx HT Λ−1

w = w
−1
= Λx HT HΛx HT + Λw HΛx HT + Λw − HΛx HT Λ−1
   
w
T
 T
−1
= Λx H HΛx H + Λw (3.349)

which agrees with (3.342).

 184 Estimation Theory Chap. 3

3.B FISHER INFORMATION CALCULATIONS FOR SINUSOID ESTIMATION

We compute the Fisher matrix entries one at a time. First,

 2 
∂
[Iy (x)]11 = −E `(y; x)
∂A2
N −1
1 X
= 2 cos2 (ω0 n + Θ)
σ n=0
N −1
1 X
= 2 [1 + cos(2ω0 n + 2Θ)]
2σ n=0
" N −1 #
N N 1 X
= 2+ 2 cos(2ω0 n + 2Θ)
2σ 2σ N n=0
N N
= 2
+ 2 Re {ξ(ω0)} , (3.350)
2σ 2σ
where we introduce the function ξ(·) defined via
N −1
1 X j(2ω0 n+2Θ)
ξ(ω0 ) = e . (3.351)
N n=0

As we’ll see, this function and its first and second derivatives, respectively
N −1
0 2 X j(2ω0 n+2Θ)
ξ (ω0 ) = ne (3.352)
N n=0

and
N −1
00 4 X 2 j(2ω0 n+2Θ)
ξ (ω0 ) = n e , (3.353)
N n=0

play a central role in the Fisher information for the problem.

Next, using (3.352),

∂2
 
[Iy (x)]12 = −E `(y; x)
∂A∂ω0
N −1
1 X
=− 2 An cos(ω0 n + Θ) sin(ω0 n + Θ)
σ n=0
" N −1 #
AN 1 X
=− 2 n sin(2ω0 n + 2Θ)
2σ N n=0
AN 1
=− Im {ξ 0(ω0 )} , (3.354)
2σ 2 2
Sec. 3.B Fisher Information Calculations for Sinusoid Estimation 185

and, using (3.351),

∂2
 
[Iy (x)]13 = −E `(y; x)
∂A∂Θ
N −1
1 X
=− 2 A cos(ω0 n + Θ) sin(ω0 n + Θ)
σ n=0
" N −1 #
AN 1 X
=− 2 sin(2ω0 n + 2Θ)
2σ N n=0
AN
=− Im {ξ(ω0 )} . (3.355)
2σ 2
Proceeding, using (3.293) and (3.353),
 2 
∂
[Iy (x)]22 = −E `(y; x)
∂ω0 2
N −1
1 X 2 2 2
=− 2 A n sin (ω0 n + Θ)
σ n=0
N −1
A2 X 2
= 2 n [1 − cos(2ω0 n + 2Θ)]
2σ n=0
N −1
A2 X 2
= 2 n [1 − cos(2ω0 n + 2Θ)]
2σ n=0
"N −1 # " N −1 #
A2 X 2 A2 N 1 X 2
= 2 n − n cos(2ω0 n + 2Θ)
2σ n=0 2σ 2 N n=0
A2 A2 N 1
= 2
TN − 2
Re {ξ 00 (ω0 )} , (3.356)
2σ 2σ 4
and, using (3.292) and (3.352),
∂2
 
[Iy (x)]23 = −E `(y; x)
∂ω0 ∂Θ
N −1
1 X 2
=− 2 A n sin2 (ω0 n + Θ)
σ n=0
N −1
A2 X
=− 2 n [1 − cos(2ω0 n + 2Θ)]
2σ n=0
"N −1 # "N −1 #
2 2
A X AN X
= 2 n − n cos(2ω0 n + 2Θ)
2σ n=0 2σ 2 n=0
A2 A2 N 1
= S N − Re {ξ 0 (ω0 )} . (3.357)
2σ 2 2σ 2 2
186 Estimation Theory Chap. 3

Finally, using (3.351),

∂2
 
[Iy (x)]33 = −E 2 `(y; x)
∂ Θ
N −1
1 X 2 2
=− 2 A sin (ω0 n + Θ)
σ n=0
N −1
A2 X
= 2 [1 − cos(2ω0 + 2Θ)]
2σ n=0
" N −1 #
A2 N A2 N 1 X
= − cos(2ω0 n + 2Θ)
2σ 2 2σ 2 N n=0
A2 N A2 N
= − Re {ξ(ω0 )} . (3.358)
2σ 2 2σ 2
To develop the asymptotic (large N) behavior of the Fisher information, we
explore the corresponding behavior of the function ξ(ω0 ) for 0 ≤ ω0 < π. To begin,
since (3.351) is a finite geometric sum, we readily obtain
1 j2Θ 1 − ej2ω0 N 1 sin ω0 N
|ξ(ω0)| = e · j2ω
= , (3.359)
N 1−e 0 N sin ω0
which is depicted in Fig. 3.13. As this figure reflects, and consistent with (3.359),
the function ξ(ω0 ) has a mainlobe at ω0 = 0 of unit height and a series of sidelobes
spaced apart by π/N; i.e.,
(
 nπ  1 n=0
ξ = . (3.360)
N 0 n = 1, 2, . . . , N − 1
The first sidelobe has its peak near ω0 = 3π/2N, and successive sidelobes get pro-
gressively smaller. In particular, when N is large, this first sidelobe has height
 
3π 1 sin(3π/2) 2
ξ = ≈ , (3.361)
2N N sin(3π/2N) 3π
while the height of the smallest sidelobe is no larger than 1/N. Indeed, using
(3.359), we see  
π  1 πN 1
ξ = sin ≤ . (3.362)
2 N 2 N
More generally, in the large N regime (3.295) we have15
ξ 0 (ω0 ) ξ 00 (ω0 )
|ξ(ω0)|  1,  1,  1, (3.363)
N N2
15
More precisely, for any fixed ω0 such that 0 < ω0 < π and given an ε > 0 that can be
arbitrarily small, there exists an N = N (ω0 ) such that
ξ 0 (ω0 ) ξ 00 (ω0 )
|ξ(ω0 )| < ε, < ε, < ε.
N N2
 Sec. 3.C Maximum Likelihood Sinusoid Estimator Derivation 187

N=9

0.8
|ξ(ω0)|

0.6

0.4

0.2

0
Figure 3.13. Plot of the magnitude of
−0.2 0 0.2 0.4 0.6 0.8 1 the function ξ(ω0 ) defined in (3.351)
normalized frequency ω0/π
when N = 9.

i.e.,
ξ(ω0 ) ∼ o(1), ξ 0(ω0 ) ∼ o(N), ξ 00 (ω0 ) ∼ o(N 2 ). (3.364)

Using (3.364) in (3.350)–(3.358) yields the desired (3.291).

3.C MAXIMUM LIKELIHOOD SINUSOID ESTIMATOR DERIVATION

To obtain our solution, we perform an invertible transformation of the parameter

set of the form
(A, ω0 , Θ) −→ (α1 , α2 , ω0 )
where
α1 = A cos Θ, α2 = A sin Θ (3.365)
and develop ML estimates (α̂1 , α̂2 , ω̂0 ) of the new parameters. Then via the invari-
ance property of ML estimation, the corresponding ML estimates (Â, ω̂0 , Θ̂) follow
immediately as
 
α̂2
q
2 2 −1
 = α̂1 + α̂2 Θ̂ = − tan . (3.366)
α̂1

A straightforward solution in terms of the new parameters involves casting

the problem as a nested optimization: we first find the optimum α̂1 and α̂2 in terms
of ω0 , then solve for the optimum ω̂0 , i.e.,
 
min J(α1 , α2 , ω0 ) = min min J(α1 , α2 , ω0 )
α1 ,α2 ,ω0 ω0 α1 ,α2

= min J(α̂1 (ω0 ), α̂2 (ω0 ), ω0 )

ω0

= J(α̂1 (ω̂0 ), α̂2 (ω̂0 ), ω̂0). (3.367)

188 Estimation Theory Chap. 3

The first minimization is straightforward because α1 and α2 appear linearly

in the objective function. In particular, using
   
1 0
 cos ω0   sin ω0 
c(ω0 ) =  , s(ω ) = (3.368)
   
..  0  .. 
 .   . 
cos ω0 (N − 1) sin ω0 (N − 1)
we can write
y = α1 c(ω0 ) + α2 s(ω0 ) + w
= H(ω0 )α + w (3.369)
where 
  α1
H(ω0) = c(ω0 ) s(ω0 ) and α= . (3.370)
α2
Hence, for a given ω0 , solving for α is a standard linear-Gaussian ML estimation
problem. In particular, we obtain,
α̂(ω0 ) = arg min J(α1 , α2 , ω0 )
α1 ,α2

= min ky − H(ω0 )αk2

α
−1
= H(ω0 )T H(ω0 ) H(ω0 )T y.

(3.371)
For the second stage of optimization, we then have
ω̂0 = arg min J1 (ω0 ) (3.372)
ω0

where
J1 (ω0 ) = J(α̂1 (ω0 ), α̂2 (ω0 ), ω0 )
= ky − H(ω0 )α̂(ω0 )k2
h −1 i
= yT I − H(ω0 ) H(ω0 )T H(ω0) H(ω0)T y.

(3.373)

Thus, from (3.373), we obtain

ω̂0 = arg min J1 (ω0 ) = arg max J2 (ω0 ) (3.374)
ω0 ω0

where −1
J2 (ω0 ) = yT H(ω0 ) H(ω0 )T H(ω0 ) H(ω0 )T y.

(3.375)
Now
c(ω0 )T c(ω0 ) c(ω0 )T s(ω0 )
 
T
H(ω0) H(ω0 ) = (3.376)
c(ω0 )T s(ω0 ) s(ω0 )T s(ω0 )
and
c(ω0 )T y
 
T
H(ω0 ) y = . (3.377)
s(ω0 )T y
Sec. 3.C Maximum Likelihood Sinusoid Estimator Derivation 189

But in the large N regime (3.295) we have, again using order notation,
N −1
X N
c(ω0 )T c(ω0 ) = cos2 ω0 n = + o(N) (3.378a)
n=0
2
N
X −1
c(ω0 )T s(ω0 ) = cos ω0 n sin ω0 n = o(N) (3.378b)
n=0
N −1
T
X N
s(ω0 ) s(ω0 ) = sin2 ω0 n = + o(N). (3.378c)
n=0
2

Hence, using (3.377) and (3.376) with (3.378) in (3.375) we obtain

−1 
c(ω0 )T y
 
 T T
 N/2 + o(N) o(N)
J2 (ω0 ) = c(ω0 ) y s(ω0 ) y
o(N) N/2 + o(N) s(ω0 )T y
2 
(c(ω0 )T y)2 + (s(ω0 )T y)2

≈
N
N −1
!2 N −1
!2 
2  X X
= y[n] cos ω0 n + y[n] sin(ω0 n) 
N n=0 n=0

= 2 |YN (ejω0 )|2 (3.379)

where YN (ejω ) is as given in (3.299). Thus, using (3.379) in (3.374), we conclude
that for N in the regime (3.295), the ML estimate ω̂0 effectively corresponds to the
location of the peak of the periodogram of the data, i.e., (3.300).
Given ω̂0 , the remaining parameters follow almost immediately. In particu-
lar, since

α̂(ω̂0 ) = (H(ω̂0 )T H(ω̂0 ))−1 H(ω̂0 )T y

2 c(ω̂0 )T y
 
≈
N s(ω̂0 )T y
2 Re YN (ej ω̂0 )
  
=√ j ω̂ , (3.380)
N Im YN (e 0 )
we have
4 2
Â2 = α̂12 (ω̂0 ) + α̂22 (ω̂0 ) = α̂(ω̂0 )T α̂(ω̂0 ) = YN (ej ω̂0 ) (3.381)
N
corresponding to (3.301), and
!
Im YN (ej ω̂0 )
  
−1 α̂2 (ω̂0 )
Θ̂ = − tan = − tan−1 (3.382)
α̂1 (ω̂0 ) Re {YN (ej ω̂0 )}

corresponding to (3.302).