Detection and Estimation
Detection and Estimation
Fall 2003
1
Probability, Random
Vectors, and Vector Spaces
1.1 INTRODUCTION
In these notes, we explore a powerful and remarkably broad framework for gener-
ating, modeling, and processing signals characterized by some degree of random-
ness or uncertainty. As we’ll see, this framework will be useful not only concep-
tually, but also practically, allowing us to develop a wide range of efficient algo-
rithms for various kinds of applications. In this first chapter of the course notes,
we develop important foundations for this framework, which we’ll build on in
subsequent chapters.
We build this foundation by combining the tools of probability theory with
concepts from the theory of vector spaces. We assume you’ve had a good deal of
exposure to basic probability concepts in your undergraduate curriculum. And
we assume you’ve also developed considerable experience with Euclidean vector
spaces and the associated tools of linear algebra from your undergraduate curricu-
lum. One role this chapter serves is to collect together and summarize those con-
cepts and techniques from this material that we will exploit extensively through-
out the course. However, the larger and more important purpose of this chapter is
to introduce and develop new ideas that arise from exploiting these ideas jointly.
As an example, we’ll develop the concept of a random vector, and explore some
important ways for characterizing such quantities. And we’ll introduce the no-
tion of abstract (non-Euclidean) vector spaces, which we’ll use in turn, to explore,
e.g., the notion of vector spaces of random variables. Some of these ideas will un-
doubtedly seem quite unusual at first, and will take some time and effort to digest.
5
6 Probability, Random Vectors, and Vector Spaces Chap. 1
However, as we’ll see they lead to some powerful geometric perspectives that will
play a key role in the course.
A detailed outline of the chapter is as follows. We begin with a compact sum-
mary of those probability concepts that will be of most use to us. Building on this
foundation, we then introduce random vectors using a vector-matrix notation, and
develop key concepts and properties. Finally, we introduce the concept of abstract
vector space, and develop several important examples of such spaces, including
those involving random variables. The accompanying appendices summarize im-
portant concepts and results from linear algebra and vector calculus that we rely
on in this and future chapters. Additional results from linear algebra and vector
space theory will be developed as we need them in subsequent chapters.
Pr [Ω] = 1 (1.2)
Pr [A ∪ B] = Pr [A] + Pr [B] if A ∩ B = ∅. (1.3)
Two of the many consequences of these axioms are
Pr [∅] = 0 (1.4)
and
Pr [A ∪ B] = Pr [A] + Pr [B] − Pr [A ∩ B] . (1.5)
Finally, (1.5) can be used with induction to establish the union bound: if the
Ai , i = 1, 2, . . . , n are an arbitrary collection of events, then
"n # n
[ X
Pr Ai ≤ Pr [Ai ] , (1.6)
i=1 i=1
where equality in (1.6) holds if and only if the Ai are a collection of mutually exclu-
sive events, i.e., if Ai ∩ Aj = ∅ for i 6= j.
1
In fact we cannot compute the probability of every subset of Ω. Those that we can we will
term valid subsets. In formal mathematical treatments a probability space is specified in terms
of a sample space, a probability measure, and a collection of valid sets. At our level of treatment,
however, you can assume that any subset we mention or construct—either explicitly or implicitly—
is valid.
Sec. 1.2 Axioms of Probability and Basic Concepts 7
1.2.2 Independence
For three events A, B, and C to be mutually independent, for example, this means
that we require that all the following hold:
Pr [A ∩ B ∩ C] = Pr [A] Pr [B] Pr [C] (1.14)
Pr [A ∩ B] = Pr [A] Pr [B] (1.15)
Pr [A ∩ C] = Pr [A] Pr [C] (1.16)
Pr [B ∩ C] = Pr [B] Pr [C] . (1.17)
In particular, (1.14) alone is not sufficient; (1.15)–(1.17) are also required.
In these notes, we adopt the useful convention of using fonts without serifs for
random variables, and the corresponding fonts with serifs for sample values and
dummy arguments. For example, x, y , z, and Θ will denote random variables, and
x, y, z, and Θ corresponding generic sample values.
Formally, a random variable x is a real-valued function on the sample space
Ω. The probability distribution function for x is defined by
Px (x) = Pr [ω | x(ω) ≤ x] = Pr [x ≤ x] , (1.18)
where the last expression we use for notational convenience. This distribution is a
complete characterization of the random variable. Likewise, the probability density
function (pdf) px (x), which is related to the distribution by
dPx (x)
px (x) = , (1.19)
dx
is also a complete characterization.2 This follows from the fact that for any valid
A, we can write Z
Pr [x ∈ A] = px (x) dx. (1.20)
A
If x takes on particular values with nonzero probability, then Px (x) will con-
tain step-discontinuities and px (x) will contain impulses. For example,
1 1
px (x) = δ(x + 1) + δ(x − 1) (1.21)
2 2
is the density of a random variable taking on the values ±1 each with the prob-
ability 1/2. To accommodate the possibility of px (x) having impulses and remain
consistent with (1.18), we write the inverse of (1.19) as
Z x+
Px (x) = px (u) du, (1.22)
−∞
2
We’ll assume in our treatment that densities always exist for the quantities of interest, at
least in this generalized sense (i.e., allowing impulses). However, it is worth keeping in mind
that there exist random variables whose probability distributions are not differentiable even in a
generalized sense.
Sec. 1.3 Random Variables 9
using x+ in the upper limit of (1.22) to indicate that the endpoint x is included in
the interval. Also, since [cf. (1.19)]
Pr [x0 < x ≤ x0 + δx]
px (x0 ) = lim ,
δx→0 δx
we have the frequently useful approximation valid for suitably small δx:
Pr [x0 < x ≤ x0 + δx] ≈ px (x0 ) δx. (1.23)
1.3.1 Expectations
Remember that since z = g(x) is itself a random variable we may also write (1.24)
in the form Z +∞
E [g(x)] = E [z] = z pz (z) dz, (1.25)
−∞
where pz (z) is the probability density for z = g(x). If g(·) is a one-to-one and
differentiable function, a simple expression for pz (z) can be derived, viz.,
px (g −1 (z))
pz (z) = (1.26)
|g ′ (g −1(z))|
where g ′ (·) denotes the first derivative of g(·). If g(·) is not invertible, the more
general method-of-events approach for deriving densities, which we briefly re-
view later in the multivariate case, can be employed to obtain pz (z). In terms of
the ultimate goal of evaluating E [g(x)], whether (1.24) or (1.25) turns out to be
more convenient depends on the problem at hand.
Several expectations that are important partial characterizations of a random
variable are the mean value (or first moment)
E [x] = x , mx , (1.27)
the mean-squared value (or second moment)
E x 2 = x 2,
(1.28)
and the variance (or second central-moment)
E (x − mx )2 = x 2 − m2x , var x , σx2 , λx .
(1.29)
In (1.27)–(1.29) we have introduced a variety of notation that will be convenient to
use in subsequent sections of these notes. The standard deviation is σx , the square
10 Probability, Random Vectors, and Vector Spaces Chap. 1
root of the variance. One important bound provided by these moments is the
Chebyshev inequality
σ2
Pr [|x − mx | ≥ ε] ≤ 2x (1.30)
ε
Observe that (1.30) implies that x is a constant (i.e., Pr [x = α] = 1 for some con-
stant α) if σx2 = 0. Note that the corresponding “only if” statement follows imme-
diately from (1.29).
The Chebyshev bound is a particularly convenient bound to use in practice
because its calculation involves only the mean and variance of the random vari-
able, i.e., it doesn’t depend on the detailed form of the density. However, for this
same reason the Chebyshev bound is not a particularly tight bound.3
and, as is apparent from the integral in (1.31), corresponds to the Fourier transform
of the density (to within a minor sign change). As a Fourier transform, we can
recover px (x) from Mx (jv) via the inverse formula
Z +∞
1
px (x) = e−jvx Mx (jv) dv. (1.32)
2π −∞
and hence the characteristic function is an equivalent complete characterization of
a random variable.
Characteristic functions are particularly useful in computing certain expec-
tations involving the random variable. For example, the moments of x can all be
efficiently recovered from Mx (jv) by differentiation, i.e.,
1 dn
n
E [x ] = n n Mx (jv) (1.33)
j dv v=0
Observe that (1.33) implies that the characteristic function can be expanded in
terms of the power series
+∞
X (jv)k
Mx (jv) = E xk (1.34)
k=0
k!
when all the moments of the form (1.33) exist. This result implies, in turn, that
knowledge of all moments is an equivalent characterization for such random vari-
ables: given these moments we can reconstruct Mx (jv) via (1.34).
3
As an aside, an alternative bound that is typically much tighter but which requires access
to more information about the random variable is the Chernoff bound.
Sec. 1.3 Random Variables 11
Random variables that take on only integer values can be fully developed within
the framework we’ve been describing. In particular, their probability densities
consist entirely of uniformly-spaced impulses with suitable weights. However, to
make manipulation of these quantities less cumbersome, it is sometimes conve-
nient to adopt some special notation for specifically discrete random variables. In
particular, we define the probability mass function (pmf) of an integer-valued ran-
dom variable k as
pk [k] = Pr [k = k] (1.37)
using square brackets to distinguish masses from densities, and to remind us that
the argument is integer-valued. The density can, of course, be derived from the
mass function via
+∞
X
pk (k) = pk [k] δ(i − k) (1.38)
i=−∞
We will frequently deal with several random variables, and we will find it conve-
nient to use vector notation in this case. Before we do that, however, let us recall
some complete joint characterizations of a pair of random variables x and y . One
such characterization is the joint distribution function for x and y , which is defined
by
Px,y (x, y) = Pr [x ≤ x and y ≤ y] . (1.40)
A second complete joint characterization is the joint density of x and y , i.e.,
∂ 2 Px,y (x, y)
px,y (x, y) = . (1.41)
∂x ∂y
If A ⊂ R2 is a valid set, where R2 denotes4 the plane of all pairs (x, y), then
Z Z
Pr [(x, y ) ∈ A] = px,y (x, y) dx dy. (1.42)
A
Again, from (1.41) we also have the following approximation valid for suitably
small δx and δy:
Pr [x0 < x ≤ x0 + δx and y0 < y ≤ y0 + δy] ≈ px,y (x0 , y0) δx δy. (1.43)
4
See the Appendix 1.A for a discussion of such spaces.
Sec. 1.4 Pairs of Random Variables 13
1.4.2 Independence
On many occasions we will exploit the fact that expectations are linear opera-
tions. For example, for arbitrary constants α and β we have
E [αx + βy ] = αE [x] + βE [y ] .
This means that in computations we can typically interchange expectations with
summations, integrations, and other linear operations.
In addition to (1.27)–(1.29) for x and their counterparts for y , some additional
expectations that constitute useful partial characterizations of the statistical rela-
tionship between x and y are
14 Probability, Random Vectors, and Vector Spaces Chap. 1
Correlation:
E [xy ] (1.51)
Covariance:
E [(x − mx )(y − my )] = E [xy ] − mx my
= λxy , cov (x, y ) . (1.52)
= E [x] E [y ] .
Z +∞
E [E [x|y]] = E [x|y = y] py (y) dy
−∞
Z +∞ Z +∞
= x px|y (x|y) dx py (y) dy
−∞ −∞
Z+∞ Z +∞
= x px,y (x, y) dx dy = E [x] (1.56)
−∞ −∞
The identity (1.56), which we’ll use on many occasions in this course, is called the
law of “iterated expectation.”
Z +∞ Z +∞
E [x|y = y] = x px|y (x|y) dx = x px (x) dx = E [x] ,
−∞ −∞
the converse is not true: x and y are not necessarily independent if (1.57) holds.
As a simple counterexample we have the joint density
1 1 1
px,y (x, y) = δ(x, y − 1) + δ(x + 1, y) + δ(x − 1, y). (1.58)
3 3 3
Likewise, it is true that x and y are uncorrelated if (1.57) holds since, using iterated
expectation, we have
however, the converse is again not true: if x and y are uncorrelated, we cannot
deduce that (1.57) holds. A simple counterexample is the density (1.58) with x and
y interchanged:
1 1 1
px,y (x, y) = δ(x − 1, y) + δ(x, y + 1) + δ(x, y − 1).
3 3 3
16 Probability, Random Vectors, and Vector Spaces Chap. 1
Rather than collecting all the random variables of interest into a single vector
x, in many problems it is often more natural and more convenient to divide them
among several random vectors of possibly different sizes.
In the case where we divide our random variables into two random vectors
x ∈ RN and y ∈ RM , we can define the joint distribution
Px,y (x, y) = Pr [x1 ≤ x1 , x2 ≤ x2 , . . . , xN ≤ xN , y1 ≤ y1 , y2 ≤ y2 , . . . , yM ≤ yM ]
(1.65)
5
Since we use bold face fonts for vectors, random vectors will be denoted using bold face
fonts without serifs, and sample values will be denoted using bold face fonts with serifs. For
example, x, y, z, and Θ will be random vectors, and x, y, z, and Θ will be associated sample
values.
Sec. 1.5 Random Vectors 17
Two random vectors x and y are independent (meaning that the two cor-
responding collections of random variables are mutually independent of one an-
other) if knowledge of any of the elements of y provides no information about any
of the elements of x (or vice versa), i.e., if
px|y (x|y) = px (x). (1.71)
Analogous to our earlier results, using (1.68) we find that (1.71) is equivalent to
the condition
px,y (x, y) = px (x) py (y). (1.72)
All of these formulas extend to more than two random vectors. For instance,
a collection of K random vectors x1 , x2 , . . . , xK are mutually independent if for
every i we have
pxi |{xj ,j∈J}(xi | {xj , j ∈ J}) = pxi (xi ), (1.73)
where J is any subset of indices between 1 through K but excluding i. The condi-
tion (1.73) is equivalent to the requirement that
n
Y
px1 ,x2 ,...,xN (x1 , x2 , . . . , xN ) = pxi (xi ). (1.74)
i=1
18 Probability, Random Vectors, and Vector Spaces Chap. 1
Note that by integrating out any subset of vectors in (1.74) we obtain lower-order
independence relations among arbitrary subsets of the random vectors as well,
i.e., if I is an arbitrary subset of distinct indices selected from 1 to K, then
Y
p{xi ,i∈I} ({xi , i ∈ I}) = pxi (xi ). (1.75)
i∈I
Suppose
y1 g1 (x)
y2 g2 (x)
y = .. = g(x) = ..
. .
yM gM (x)
is an M-dimensional random vector obtained as a function of the N-dimensional
random vector x. We can always in principle calculate the distribution for y from
the method of events:
Py (y) = Pr [g1 (x) ≤ y1 , g2 (x) ≤ y2 , . . . , gM (x) ≤ yM ]
Z
= px (x) dx (1.79)
A(y)
where
A(y) = {x | g1 (x) ≤ y1 , g2 (x) ≤ y2 , . . . , gM (x) ≤ yM } (1.80)
We can then obtain the density via
∂ M Py (y)
py (y) = . (1.81)
∂y1 ∂y2 · · · ∂yM
then
E [f1 (x)]
E [f2 (x)]
E [f(x)] = . (1.85)
..
.
E [fM (x)]
Some important expectations are:
Mean Vector:
E [x] = mx (1.86)
Correlation Matrix:
E xxT
(1.87)
Covariance Matrix:
cov (x, x) = Λxx = E (x − mx )(x − mx )T = E xxT − mx mT
x (1.88)
Cross-Covariance Matrix:
cov (x, y) = Λxy = E (x − mx )(y − my )T = E xyT − mx mT
y (1.89)
Conditional Mean:
Z +∞
mx|y (y) = mx|y=y = E [x|y = y] = x px|y (x|y) dx (1.90)
−∞
Conditional Covariance:
Λx|y (y) = Λx|y=y
Z +∞
= (x − E [x|y = y]) (x − E [x|y = y])T px|y (x|y) dx (1.91)
−∞
As before we can think of the conditional statistics mx|y and Λx|y in (1.90)
and (1.91), respectively, as deterministic quantities that are functions of a partic-
ular value y = y. Alternatively, mx|y and Λx|y can be viewed as functions of y
and therefore random variables in their own right. As before, the law of iterated
expectation applies, i.e.,
E [E [x|y]] = E [x] .
For notational convenience, we will often drop one of the subscripts in deal-
ing with the covariance of a random vector x, i.e., we will often write Λx in-
stead of Λxx . In terms of dimensions, note that if x is N-dimensional and y is
M-dimensional, then Λx is N × N, Λxy is N × M, and Λyx is M × N. Furthermore
the (i, j)th and (i, i)th elements of Λx are
[Λx ]ij = cov (xi , xj ) (1.92)
[Λx ]ii = σx2i (1.93)
Sec. 1.5 Random Vectors 21
In turn, (1.101) implies that Mx (jv) can be expanded in a power series of the form
+∞ X
+∞ +∞ k1 k2
X X h
k1 k2 kN (jv1 ) (jv2 )
i (jvN )kN
Mx (jv) = ··· E x1 x2 · · · xN ··· , (1.102)
k1 =0 k2 =0 kN =0
k1 ! k2 ! kN !
22 Probability, Random Vectors, and Vector Spaces Chap. 1
provided all the constituent moments exist. Hence, many classes of random vec-
tors are completely characterized by the complete set of moments of the form
(1.101).
In addition, note that the collection of random variables x1 , x2, . . . , xN are mu-
tually independent if and only if
Mx (jv) = Mx1 (jv1 ) Mx2 (jv2 ) · · · MxN (jvN ). (1.103)
To establish the “only if” part, it suffices to note that if the x1 , x2, . . . , xN are mutu-
ally independent then
h T i
Mx (jv) = E ejv x = E ej(v1 x1 +v2 x2 +···+vN xN )
As a special case, let a be a vector of numbers and consider the scalar random
variable
XN
T
z =a x= ai xi . (1.111)
i=1
Example 1.1
Let x and y be two scalar random variables, and consider the random vector
x
w= . (1.113)
y
Then
σx2 σx2
λxy ρxy σx σy
Λw = 2 = . (1.114)
λxy σy ρxy σx σy σy2
For Λw to be positive definite, it must be true that the determinant of Λw is positive:
det(Λw ) = (1 − ρ2xy )σx2 σy2 > 0. (1.115)
From this equation we can see that Λw will not be positive definite if and only if the
correlation coefficient ρxy equals ±1. In either of these cases we can conclude that x
must equal a multiple of y plus a constant, i.e., x = cy + d for some constants c and
d. It is straightforward to check that the sign of c is the same as that of ρxy .
6
See Appendix 1.A for a discussion of positive definite and semidefinite matrices.
24 Probability, Random Vectors, and Vector Spaces Chap. 1
Example 1.2
Continuing Example 1.1, suppose that
3/2 1/2
Λw = .
1/2 3/2
Then the eigenvalues are λ1 = 2 and λ2 = 1, and the corresponding normalized
eigenvectors are
√ √
1/√2 1/ √2
p1 = p2 = .
1/ 2 −1/ 2
Hence, we can conclude that the pair of random variables
√
u = 2 pT w =x +y
√ 1T
v = 2 p2 w = x − y
are uncorrelated and have variances
x +y
var u = var [x + y ] = 2 var √ = 2λ1 = 4
2
x −y
var v = var [x − y ] = 2 var √ = 2λ2 = 2.
2
In this section we define and develop the basic properties of jointly Gaussian or
normal random variables (or, equivalently, Gaussian or normal random vectors).
Gaussian random variables are important for at least two reasons. First, Gaus-
sian random vectors are good models in many physical scenarios. For example,
Sec. 1.6 Gaussian Random Variables 25
px (x)
1/√2πσ2
1/√2πeσ2
One fairly general form of the Central Limit Theorem is stated formally as follows.
Let
x1 , x2 , x3 , . . .
be a sequence of mutually independent zero-mean random variables with distri-
butions
Px1 (x1 ), Px2 (x2 ), Px3 (x3 ), . . .
and variances
σ12 , σ22 , σ32 , . . . ,
respectively. If for any ǫ > 0 there exists a k (depending on ǫ) sufficiently large that
σi < ǫ Sk for i = 1, 2, . . . , k (1.124)
with v
u k
uX
Sk = t σi2 ,
i=1
A couple of points are worth emphasizing. First, the somewhat exotic con-
straint (1.124) essentially ensures that no one term dominates the sum (1.125).7 In
7
To see that this constraint is critical, it suffices to consider the sequence of independent
Bernoulli random variables xi , each of which is ±1/2i with equal probability. Note that this se-
quence does not satisfy (1.124). For this sequence, it is straightforward to verify using (1.36) that
the distribution of the normalized sum (1.125) converges to a uniform rather than Gaussian distri-
bution.
Sec. 1.6 Gaussian Random Variables 27
fact, a simpler special case of this theorem corresponds to the xi being identically-
distributed and having a finite common variance. Second, it is important to em-
phasize that the theorem guarantees convergence in distribution but not in density.
In fact, when the random variables in the sum are discrete, it is impossible to have
convergence in density since arbitrary partial sums will be discrete!
The notion of a Gaussian random vector is a powerful and important one, and
builds on our notion of a Gaussian random variable. Specifically, an N-dimensional
28 Probability, Random Vectors, and Vector Spaces Chap. 1
random vector
x1
x2
x = ..
.
xN
is defined to be a Gaussian random vector, or equivalently {x1 , x2 , . . . , xN } is de-
fined to be a set of jointly Gaussian random variables when for all choices of the
constant vector
a1
a2
a = .. (1.130)
.
aN
the scalar y = aT x is a Gaussian random variable.
Gaussian random vectors have several important properties. In what fol-
lows, suppose x is a Gaussian random vector whose mean is mx and whose co-
variance matrix is Λx .
First, all subsets of {x1 , x2 , . . . , xN } are jointly Gaussian. Deriving this result
simply requires setting some of the ai ’s in (1.130) to zero. As a special case of
this result—corresponding to having only one nonzero component in (1.130)—we
have that all the constituents must be individually Gaussian random variables,
i.e.,
xi ∼ N(mi , λii ) for i = 1, 2, . . . , N
where
λii = [Λx ]ii .
The characteristic function for a Gaussian random vector takes the form
T 1 T
Mx (jv) = exp jv mx − v Λx v . (1.131)
2
To prove (1.131), first note that
h T i
Mx (jv) = E ejv x = E ej(v1 x1 +v2 x2 +···+vN xN )
(1.132)
The density (1.137) can be obtained via direct computation of the inverse
Fourier transform of (1.131); a derivation is as follows. First, observe that
Z +∞
1 T
px (x) = N
Mx (jv)e−jv x dv
(2π) −∞
Z +∞
1 T 1 T
= exp −jv (x − mx ) − v Λx v dv. (1.140)
(2π)N −∞ 2
Then, using the change of variables u = Λ1/2
x v with Λx
1/2
as defined in (1.139), and
noting that the Jacobian of the transformation is, using (1.139b), |du/dv| = |Λx1/2 | =
|Λx |1/2 , we can rewrite (1.140) as
Z +∞
1 T −1/2 1 T du
px (x) = N
exp −ju Λx (x − mx ) − u u ,
(2π) −∞ 2 |Λx |1/2
which when we adopt the convenient notation
x̃ = Λ−1/2
x (x − mx ) (1.141)
8
We have used Λx−1/2 to denote the inverse of this square root matrix. Incidently, it is
straightforward to verify that this matrix is also the positive definite square root matrix of Λ−1
x .
30 Probability, Random Vectors, and Vector Spaces Chap. 1
(1.142)
Finally, recognizing that each of the integrals in (1.142) is unity, and replacing x̃
with its definition (1.141) we obtain, after some simple manipulations, our desired
result (1.137).
Several additional important properties of Gaussian random vectors are worth de-
veloping. First, a pair of jointly Gaussian random vectors x and y are independent
if and only if they are uncorrelated. We established the “only if” part for any pair
of random vectors earlier. To establish the “if” part, let
x
z= (1.143)
y
so that z ∼ N(mz , Λz ), with
mx
mz = (1.144)
my
Λx Λxy
Λz = . (1.145)
Λyx Λy
Then when x and y are uncorrelated, i.e., when Λxy = 0, we have
det Λz = det Λx det Λy (1.146)
and
Λ−1
0
Λ−1 = x
. (1.147)
z
0 Λ−1
y
Using the expression (1.137) for the Gaussian density with these results, one can
easily check that in this case
pz (z) = px,y (x, y) = px (x) py (y). (1.148)
9
Note from (1.150) that Λx|y (y) is a constant matrix, i.e., independent of the value of y.
32 Probability, Random Vectors, and Vector Spaces Chap. 1
where the summation in (1.152) is over all distinct pairings {j1 , j2 }, {j3 , j4 }, . . . ,
{jL−1 , jL } of the set of symbols {i1 , i2 , . . . , iL }. Although we won’t develop it here,
this result may be derived in a relatively straightforward manner using, e.g., a
Taylor series expansion of Mx (jv). As an example application of (1.152) we have
E [x̃i1 x̃i2 x̃i3 x̃i4 ] = λi1 i2 λi3 i4 + λi1 i3 λi2 i4 + λi1 i4 λi2 i3 (1.153)
so that
E [x̃1 x̃2 x̃3 x̃4 ] = λ12 λ34 + λ13 λ24 + λ14 λ23 (1.154)
E x̃12 x̃22 = λ11 λ22 + 2λ212
(1.155)
E x̃14 = 3λ211 .
(1.156)
As a final remark, we point out that the detailed shape of the contours of
equiprobability for the multidimensional Gaussian density can be directly de-
duced from geometry of the covariance matrix as developed in Section 1.5.4. In
particular, from (1.137) we see that the contours of equiprobability are the N-
dimensional ellipsoids defined by
(x − mx )T Λ−1
x (x − mx ) = constant. (1.157)
From this perspective the transformation (1.116) from x to z corresponds to a gen-
eralized coordinate rotation (i.e., length preserving transformation) such that the
components of z represent the principal or major axes of this ellipsoid, i.e.,
(z1 − mz1 )2 (z2 − mz2 )2 (zN − mzN )2
+ +···+ = constant.
λ1 λ2 λN
Note that the λi = var zi describe the proportions of the ellipsoid: they correspond
to (squares of) the relative lengths along the principal axes. Note too that since z
is Gaussian, its components are not only uncorrelated but mutually independent
random variables.
We conclude this section by specializing our results to the case of two-dimensional
Gaussian random vectors, where we let
2
m1 σ1 ρσ1 σ2
mx = Λx = . (1.158)
m2 ρσ1 σ2 σ22
Here
1 1 T −1
px (x) = exp − (x − mx ) Λx (x − mx ) (1.159)
(2π)N/2 |Λx |1/2 2
h i
(x −m ) σ −2(x1 −m1 )(x2 −m2 )ρσ1 σ2 +(x2 −m2 )2 σ12
2 2
exp − 1 1 2 2σ12 σ22 (1−ρ2 )
= (1.160)
2πσ1 σ2 (1 − ρ2 )1/2
Fig. 1.2 depicts the joint density of a pair of Gaussian random variables. In
Fig. 1.3 we have plotted the associated contours of constant values of px (x) which
are the ellipses
(x1 − m1 )2 σ22 − 2(x1 − m1 )(x2 − m2 )ρσ1 σ2 + (x2 − m2 )2 σ12 = constant. (1.161)
Sec. 1.7 Abstract Vector Space, and Spaces of Random Variables 33
z2-m2′
z1-m1′
x2
m2
As indicated in the figure, the components of z define the principal axes of the
ellipses in (1.161), i.e., this equation in the transformed coordinates becomes
(z1 − m′1 )2 (z2 − m′2 )2
+ = constant
λ1 λ2
where λ1 and λ2 are the eigenvalues of Λx and where
′
m1
mz = = Pmx .
m′2
The larger |ρ| is, the more eccentric these ellipses become, degenerating to lines
when |ρ| = 1.
The notion of a vector space is very powerful, and one that we will exploit on
numerous occasions throughout the course. Clearly, we’ve already used certain
34 Probability, Random Vectors, and Vector Spaces Chap. 1
vector space ideas in preceding sections exploiting results from Appendix 1.A. In
particular, we’ve exploited properties of the Euclidean space RN consisting of N-
dimensional vectors. However, while Euclidean space is an important example
of a vector space, there are in fact many other somewhat more abstract vector
spaces that turn out to be at least as important to us in this course. Although more
abstract, many properties carry over from the Euclidean case, and you will often
be able to rely on the geometric picture and intuition you have developed for this
case.
Most generally, a vector space is a collection of elements or objects satisfy-
ing certain properties. This collection of elements may indeed consist of vectors
x as we usually think of them, or they may be other kinds of objects like whole
sequences x[n] or functions x(t), or even random variables x(ω). To avoid a con-
ceptual bias, we’ll just use the generic notation x for one such element.
For our purposes, vector spaces are special classes of metric spaces—i.e., spaces
in which there is some notion of distance between the various elements in the col-
lection.10 A metric space is described by the pair (S, d(·, ·)) where S is the collection
of elements and d(·, ·) is referred to as the metric. It is a measure of distance be-
tween an arbitrary pair of elements in the set; in particular d(x, y) is the distance
between elements x and y in S.
For a metric to be useful, it must satisfy certain key properties that are con-
sistent with our intuition about what distance is. In particular, we must have, for
any elements x, y, and z in S,
d(x, y) ≥ 0 (1.162)
d(x, y) = 0 ⇔ x = y (1.163)
d(x, y) = d(y, x) (1.164)
d(x, y) ≤ d(x, z) + d(z, y) (1.165)
The last of these, i.e., (1.165), is referred to as the triangle inequality.
An obvious (but not unique) example of a metric in RN is the usual Euclidean
distance v
u N
uX
d(x, y) = t (xn − yn )2 . (1.166)
n=1
where xn and yn are the nth elements of x and y, respectively. You can verify that
(1.166) satisfies (1.162)–(1.165).
The metric spaces we’re usually interested in have additional structure.
First, we want to work with spaces that are complete. While the technical def-
inition is beyond the scope of our treatment here, in essence completeness means
10
As a note of caution, our use of the term “vector space” is not universal. Some references
consider the term to be equivalent to the term “linear space.” However, as will become apparent,
we will find it convenient to define vector spaces as linear spaces that are also metric spaces.
Sec. 1.7 Abstract Vector Space, and Spaces of Random Variables 35
the metric space has no “holes.” An example of a metric space that isn’t complete
is S = (0, 1] ⊂ R with d(x, y) = |x−y|. Note that the sequence of elements xn = 1/n
for n = 1, 2, . . . , are all in the space S, but limn→∞ xn = 0 is not. The sequence xn
is an example of what is called a Cauchy sequence, and for a metric space to be
complete, all such Cauchy sequences must converge to an element of S.
A vector space V is a metric space that is linear. In order to talk about lin-
earity, we’ll need to define addition and scalar multiplication operators for objects
in V. For the cases of interest to us, we’ll be using the usual definitions of these
operators. We say V is a vector space if the following two properties hold:11
x, y ∈ V ⇒ x + y ∈ V (1.167)
x ∈ V, α ∈ R ⇒ αx ∈ V. (1.168)
There are lots of important examples of vector spaces. First, there is the usual
Euclidean space RN composed of N-dimensional vectors x. There is also the space
of sequences x[n] with finite energy
∞
X
x2 [n] < ∞
n=−∞
which is usually denoted ℓ2 (Z), and the space of (integrable) functions x(t) with
finite energy
Z +∞
x2 (t) dt < ∞
−∞
which is usually denoted L2 (R). And there is the space of random variables x(ω)
with finite mean-square
Z +∞
2
var x = E x = x2 px (x) dx < ∞,
−∞
11
When the scalar α is restricted to be a real number as (1.168) indicates, the result is referred
to as a real vector space; when it can be a complex number, i.e., α ∈ C, the result is a complex
vector space. Although we will largely focus on the former class in this course to simplify our
development, we remark in advance that we will sometimes need to work with complex vector
spaces. Fortunately, however, there are no significant conceptual differences between the two.
12
Incidently, for every probability space, there is an associated vector space of such random
variables.
36 Probability, Random Vectors, and Vector Spaces Chap. 1
A linear transformation L(·) is a linear mapping from one vector space V to another
vector space U. This means that the powerful principle of superposition is satisfied,
i.e., if xk for k = 1, 2, . . . , K are each elements of V, and if αk for k = 1, 2, . . . , K are
scalars, then
L(α1 x1 + α2 x2 + · · · + αK xK ) = α1 y1 + α2 y2 + · · · + αK yK
where yk = L(xk ).
When the vector spaces are the familiar Euclidean spaces, e.g., V = RN and
U = RM , then L(·) is represented by a matrix, i.e.,
y = L(x) = Ax
where A is an M × N-dimensional matrix. Several properties of matrices are de-
veloped in Appendix 1.A.
1.7.4 Bases
Note that when this set is linearly independent, the αk must be unique.
All bases for a vector space have the same cardinality. This cardinality is
referred to as the dimension of the vector space. As you’ve seen, the Euclidean
space V = RN has dimension N. Hence, these spaces are finite-dimensional.
Other spaces, like the space of finite-energy sequences ℓ2 (Z) and the space of finite-
energy functions L2 (R) are infinite-dimensional. The space of finite mean-square
random variables L2 (Ω) is not only infinite-dimensional, but its dimension is un-
countable (unless the probability space is discrete)! Although infinite-dimensional
spaces are difficult to visualize, much intuition from finite-dimensional Euclidean
space carries over. Furthermore, in many problems involving these vector spaces,
we will often work with finite-dimensional subspaces, for which our geometric
pictures are well-developed.
A normed vector space is a special kind of vector space for which the concept of
length is defined for elements of the space. Let us use kxk to denote the length or
norm of each x ∈ V, so a normed vector space is defined by specifying the pair
(V, k · k). In order for a function k · k to make sense as a norm on V it must satisfy
certain properties. In particular, for x ∈ V and α an arbitrary scalar, it must satisfy:
kxk ≥ 0 (1.170)
kx + yk ≤ kxk + kyk (1.171)
kαxk = |α|kxk (1.172)
kxk = 0 ⇔ x = 0 (1.173)
For normed vector spaces, the following rather natural metric can be defined:
for x and y in V,
d(x, y) = kx − yk. (1.174)
As examples, RN , ℓ2 (Z), L2 (R), and L2 (Ω) are all normed vector spaces. The
corresponding norms are defined by, respectively,
N
X
2
kxk = x2n
n=1
X
2
kx[·]k = x2 [n]
n
Z
kx(·)k2 = x2 (t) dt
kx(·)k2 = E x 2 .
Note however that there are many other norms one can define even for vectors in
RN ; for example,
kxk = max |xn |.
1≤n≤N
is a valid norm, and defines a whole family of normed vector spaces Lp (R) param-
eterized by p. We emphasize that to fully specify a normed vector space we need
both a collection of elements and a norm.
Ultimately, we’re interested in normed vector spaces with even more struc-
ture, as we’ll now develop.
An inner product space is a normed vector space where there is a notion of relative
orientation or “angle” between elements. We use the notation hx, yi to denote the
inner product between two elements x and y in V. An inner product space is
therefore defined by the pair (V, h·, ·i). An inner product defines the operation
of projection of one element onto another. A valid inner product must satisfy the
following properties14 :
14
For simplicity, we’ll restrict our attention to real-valued inner products even though many
important examples are complex-valued.
Sec. 1.7 Abstract Vector Space, and Spaces of Random Variables 39
For each inner product space there is a natural notion of norm. We call this
the induced norm, and it is defined in terms of the inner product as follows:
p
kxk = hx, xi. (1.179)
In turn, from the induced norm we get the associated metric
p
d(x, y) = hx − y, x − yi.
From the inner product and the induced norm, we arrive at a definition of
the angle θ between two elements x, y ∈ V. In particular, we have
hx, yi
cos θ = . (1.180)
kxkkyk
One enormously useful inequality that applies to inner product spaces is the
Cauchy-Schwarz inequality: for any x and y in V,
| hx, yi | ≤ kxkkyk (1.181)
A proof is as follows. For any α, we have, from the properties of a norm and
(1.179),
hx − αy, x − αyi = kx − αyk2 ≥ 0. (1.182)
Exploiting (1.175)–(1.178), we can rewrite the left hand side of (1.182) to get
hx − αy, x − αyi = kxk2 − 2α hx, yi + α2 kyk2 ≥ 0 (1.183)
Then if we let α = hx, yi / hy, yi, (1.183) becomes
2 hx, yi2
kxk − ≥0 (1.184)
kyk2
which can be rewritten in the form (1.181). As a final comment, note from (1.182)
and (1.178) that equality in (1.181) holds if and only if x − αy = 0 for an arbitrary
α.
Inner product spaces that are complete also have a special and arcane name—
they are referred to as Hilbert spaces.
As examples, RN , ℓ2 (Z), L2 (R), and L2 (Ω) are also all complete inner product
(i.e., Hilbert) spaces. The corresponding inner products are defined by, respec-
tively,
N
X
hx, yi = xT y = xn yn
n=1
X
hx[·], y[·]i = x[n]y[n]
n
Z
hx(·), y(·)i = x(t) y(t) dt
With inner product spaces we have enough structure that we can finally talk about
the concept of orthogonality. Specifically, we say that elements x and y in V are
orthogonal, denoted x ⊥ y, when their inner product is zero, i.e.,
x ⊥ y ⇔ hx, yi = 0. (1.186)
where the αk are projections of x onto the basis functions xk , i.e., αk = hx, xk i.
Sec. 1.A Linear Algebra and Euclidean Vector Space 41
In this course vectors will be matrices that are specifically columns, and as such
will be denoted by boldface lowercase characters; for example,
x1
x2
x = .. (1.187)
.
xn
where x1 , x2 , . . . , xn are either real or complex numbers. The set of all such n-
dimensional vectors of real numbers is denoted by Rn . The corresponding set of
all n-dimensional vectors of complex numbers is denoted by Cn . The transpose of
a column vector x is the row vector
xT = x1 x2 · · · xn .
(1.188)
Vector addition and scalar multiplication are defined componentwise, i.e.,
x1 y1 x1 + y1
x2 y2 x2 + y2
.. + .. = .. (1.189)
. . .
xn yn xn + yn
42 Probability, Random Vectors, and Vector Spaces Chap. 1
and
x1 αx1
x2 αx2
α .. = .. (1.190)
. .
xn αxn
where α is a real or complex number.
A set of vectors x1 , x2 , . . . , xr in Rn is linearly independent if15
α1 x1 + α2 x2 + · · · + αr xr = 0 (1.191)
implies that
α1 = α2 = · · · = αr = 0. (1.192)
Otherwise the set of vectors is said to be linearly dependent, and in this case one of
the xi can be written as a linear combination of the others. For example, if α1 6= 0 in
(1.191)
x1 = β2 x2 + · · · + βr xr (1.193)
with
β2 = −α2 /α1 , . . . , βr = −αr /α1 . (1.194)
In Rn there exist sets of at most n linearly independent vectors. Any such set
{x1 , x2 , . . . , xn } forms a basis for Rn . That is, any x ∈ Rn can be written as a linear
combination of x1 , x2 , . . . , xn .
Matrices will in general be denoted by boldface uppercase characters. The
element in the ith row and jth column of A will be denoted by aij or, alternatively,
by [A]ij . If A is m × n, i.e., if A has m rows and n columns, then
a11 a12 · · · a1n
a21 a22 · · · a2n
A = .. .. . (1.195)
.. . .
. . . .
am1 am2 · · · amn
15
The symbol 0 denotes the matrix or vector of appropriate dimension, all of whose compo-
nents are zero.
Sec. 1.A Linear Algebra and Euclidean Vector Space 43
compute |A| by “expanding by minors” using any row or column. For example,
using the ith row,
|A| = ai1 Ai1 + ai2 Ai2 + · · · + ain Ain , (1.214)
or, using the jth column,
|A| = a1j A1j + a2j A2j + · · · + anj Anj , (1.215)
where the cofactors Aij are given by
Aij = (−1)i+j det(Mij ) (1.216)
and where Mij is the (n − 1) × (n − 1) matrix obtained from A by deleting the ith
row and jth column.
As a simple example, we have
a11 a12
= a11 a22 − a12 a21 . (1.217)
a21 a22
As a more complex example we have
2 0 0 3
1 1 0 0
1 1 1 0
5 1 1 9
1 0 0 1 0 0
1+1 1+2
= 2(−1) 1 1 0 + 0(−1) 1 1 0
1 1 9 5 1 9
1 1 0 1 1 0
1+3 1+4
+ 0(−1) 1 1 0 + 3(−1) 1 1 1
5 1 9 5 1 1
1 0 1 1 1 1
= 2 · 1 · (−1)1+1 − 3 · 1 · (−1)1+1 − 3 · 1 · (−1)1+2
1 9 1 1 5 1
= 2 · 9 − 3 · 0 + 3 · (−4) = 6
1. |A| =
6 0
46 Probability, Random Vectors, and Vector Spaces Chap. 1
Consequently, we see that P is orthogonal if and only if its columns are orthonormal,
i.e., if xi ⊥ xj for i 6= j, and if kxi k = 1.
There are also some useful results for block matrices. For example, for a block
diagonal matrix
A = diag(F1 , F2 , . . . , Fr ) ⇔ A−1 = diag(F−1 −1 −1
1 , F2 , . . . , Fr ). (1.230)
Also, we have the formulas
−1
A11 A12
=
A21 A22
(A11 − A12 A−1 −1 −1 −1 −1
22 A 21 ) −(A 11 − A 12 A 22 A 21 ) A 12 A 22
−A−1 −1
22 A21 (A11 − A12 A22 A21 )
−1
A−1 −1 −1 −1
22 + A22 A21 (A11 − A12 A22 A21 ) A12 A22
−1
(1.231)
Sec. 1.A Linear Algebra and Euclidean Vector Space 47
and
A11 A12
det = A11 − A12 A−1
22 A21 |A22 | , (1.232)
A21 A22
which are valid if A22 is nonsingular. These formulas can be verified by exploiting
the identity
(1.234)
Other useful results are obtained by comparing (1.231) and (1.234). For ex-
ample, equating the upper left blocks in these two expressions yields the useful
identity
(A11 − A12 A−1
22 A21 )
−1
= A−1 −1 −1 −1 −1
11 + A11 A12 (A22 − A21 A11 A12 ) A21 A11 . (1.235)
complex conjugate pairs. However, if A is symmetric, the λi are always real. Also
note that
Y n
n n
|A| = (−1) φA (0) = (−1) α0 = λi (1.241)
i=1
so that A is invertible if and only if all of the eigenvalues of A are nonzero. In
addition, one can show that
Xn
tr(A) = −αn−1 = λi . (1.242)
i=1
u = Px (1.249)
v = Py, (1.250)
so that (since x = P−1 u) each component of u, for example, is a weighted sum of
components of x and vice versa. Then
v = Bu (1.251)
with B as given in (1.247). Furthermore,
16
Beware, however—there is no corresponding test for positive semidefiniteness that in-
volves
0 0 examining upper submatrices for nonnegative determinants. Consider, e.g., the matrix
0 −1 which is not positive semidefinite.
Sec. 1.A Linear Algebra and Euclidean Vector Space 51
1.A.6 Subspaces
The dimension of a subspace equals the maximum number of vectors in S that can
form a linearly independent set.
Let K be any subset of Rn . The orthogonal complement of K in Rn is defined as
follows:
K⊥ = {x ∈ Rn | x ⊥ y for all y ∈ K}. (1.271)
Note that K⊥ is a subspace whether or not K is, since if x1 , x2 ∈ K⊥ and y ∈ K,
(x1 + x2 )T y = xT T
1 y + x2 y = 0 (1.272)
(αx1 )T y = αxT1y = 0 (1.273)
so x1 + x2 ∈ K⊥ and αx1 ∈ K⊥ .
Let d be a single nonzero vector in Rn (so {d} is not a subspace), and consider
{d}⊥ . This is a subspace of dimension n − 1. For example, as illustrated in Fig. 1.4,
when n = 2 the set of x such that dT x = 0 is a line through the origin perpendicular
to d. In 3-dimensions this set is a plane through the origin, again perpendicular to
d. Note that the subspace {d}⊥ splits Rn into two half-spaces, one corresponding to
those x for which dT x > 0, the other to dT x < 0.
For additional insights into the concepts and results summarized in this sec-
tion, see, e.g., G. S. Strang, Linear Algebra and its Applications, 3rd ed., Academic
Press, New York, 1988.
Several results from vector calculus, which we briefly summarize here, will prove
useful. First, consider a scalar function of a vector of n real variables
x1
x2
f (x) = f .. = f (x1 , x2 , . . . , xn ). (1.274)
.
xn
17
Here R equals the set of real numbers.
Sec. 1.B Vector Calculus 53
Partial derivatives, integrals, etc., can all be defined in a useful manner. For exam-
ple, it is convenient to define a Jacobian row vector, which consists of first partial
derivatives:
df h
∂f ∂f ∂f
i
(x) = ∇x f (x) = ∂x (x) ∂x2
(x) · · · ∂xn
(x) . (1.275)
dx 1
Note that the Hessian is a symmetric matrix. Furthermore, the Hessian matrix at
x = x0 is positive semidefinite, i.e., d2 f /dx2 (x0 ) ≥ 0 whenever x0 corresponds to
a local minimum of f (·). Similarly, if x = x0 is the location of a local maximum of
f (·), then the Hessian satisfies d2 f /dx2 (x0 ) ≤ 0 which means that −d2 f /dx2(x0 ) is
positive semidefinite.
Using the notation (1.275) and (1.276) we can conveniently express the mul-
tivariable Taylor’s series expansion as
df 1 d2 f
f (x + δx) = f (x) + (x)δx + (δx)T 2 (x)δx + · · · (1.277)
dx 2! dx
where · · · in (1.277) denotes higher order terms.
Finally, we briefly discuss vector-valued functions f(·). Derivatives, inte-
grals, limits, etc., for functions of this type are defined component-wise, e.g., for a
54 Probability, Random Vectors, and Vector Spaces Chap. 1
18
Note that (1.279) is consistent both with (1.275) when m = 1 and with (1.278) when n = 1.
2
55
56 Detection Theory, Decision Theory, and Hypothesis Testing Chap. 2
ing {H0 , H1 , . . . , HM −1}.1 For each of the possible hypotheses, there is a different
model for the observed data, and this is what we will exploit to distinguish among
the hypotheses.
In this chapter, we will restrict our attention to the case in which the observed
data can be represented as a K-dimensional random vector
T
y = y1 y2 · · · yK , (2.1)
with the hypothesis. Our binary communication example fits within this class,
with M = 2 and the a priori probabilities typically being equal. In this case, the
model for the observed data under each hypothesis takes the form of a conditional
probability density, i.e., py|H (y|Hm ) for m = 0, 1, . . . , M −1. As we’ll see, in practice
these conditional probabilities are often specified implicitly rather than explicitly,
and must be inferred from the available information.
In other cases, it is more appropriate to view the valid hypothesis not as a
random variable, but as a deterministic but unknown quantity, which we denote
simply by H. In these situations, a priori probabilities are not associated with the
various hypotheses. The radar detection problem mentioned above is one that
is often viewed this way, since there is typically no natural notion of the a priori
probability of an aircraft being present. For these tests, while the valid hypoth-
esis is nonrandom, the observations still are, of course. In this case, the proba-
bility density model for the observations is parameterized by the valid hypothesis
rather than conditioned on it, so these models are denoted using py (y; Hm), for
m = 0, 1, . . . , M − 1. As in the random hypothesis case, these densities are also
often specified implicitly.
This chapter explores methods applicable to both random and nonrandom
hypothesis tests. However, we begin by focusing on random hypotheses to de-
velop the key ideas, and restrict attention to the binary (M = 2) case.
1
Note that H0 is sometimes referred to as the “null” hypothesis, particularly in asymmetric
problems where it has special significance.
Sec. 2.1 Binary Random Hypothesis Testing: A Bayesian Approach 57
P0 = Pr [H = H0 ]
(2.2)
P1 = Pr [H = H1 ] = 1 − P0 .
These summarize our state of knowledge about the applicable hypothesis before
any observed data is available.
The second is the measurement model, corresponding to the probability den-
sities for y conditioned on each of the hypotheses, i.e.,
H0 : py|H (y|H0)
(2.3)
H1 : py|H (y|H1).
The observation densities in (2.3) are often referred to as likelihood functions. Our
choice of notation suggests that y is continuous-valued; however, y can equally
well be discrete-valued, in which case the corresponding probability mass func-
tions take the form
py|H [y|Hm ] = Pr [y = y | H = Hm ] . (2.4)
For simplicity of exposition, we start by restricting our attention to the continuous
case. Again, it is important to emphasize in many problems this measurement
model information is provided indirectly, as the following example illustrates.
Example 2.1
As a highly simplified scenario, suppose a single bit of information m ∈ {0, 1} is to
be sent over a communication channel by transmitting the scalar sm , where s0 and
s1 are both deterministic, known quantities. Let’s further suppose that the channel
is noisy; specifically, what is received is
y = sm + w ,
Z0 = {y | Ĥ(y) = H0 }
(2.6)
Z1 = {y | Ĥ(y) = H1 }.
and where the expectation in (2.9) is over both y and H, and f (·) is a generic deci-
sion rule.
Often, the context of the specific problem suggests how to choose the costs
Cij . For example, a symmetric cost function of the form Cij = 1 − δ[i − j], i.e.,
C00 = C11 = 0
(2.10)
C01 = C10 = 1
Z1
Z1
Z0
that he or she does, we would typically want to select cost assignments such that
C01 ≫ C10 .2
Having chosen suitable cost assignments, we proceed to our solution by con-
sidering an arbitrary but fixed decision rule f (·). In terms of this generic f (·), the
Bayes risk can be expanded in the form
h i
ϕ(f ) = E C̃(H, f (y))
h i
= E E C̃(H, f (y)) | y = y
Z
= ϕ̃(f (y), y) py(y) dy, (2.11)
with h i
ϕ̃(H, y) = E C̃(H, H) y = y , (2.12)
and where to obtain the second equality in (2.11) we have used iterated expecta-
tion.
From the last equality in (2.11) we obtain a key insight: since py (y) is nonneg-
ative, it is clear that we will minimize ϕ if we minimize ϕ̃(f (y), y) for each particular
value of y. The implication here is that we can determine the optimum decision rule
Ĥ(·) on a point by point basis, i.e., Ĥ(y) for each y.
Let’s consider a particular (observation) point y = y∗ . For this point, if we
choose the assignment
Ĥ(y∗ ) = H0 ,
2
In still other problems, it is difficult to make meaningful cost assignments at all. In this case,
the Neyman-Pearson framework developed later in the chapter is more natural than the Bayesian
framework we develop in this section.
60 Detection Theory, Decision Theory, and Hypothesis Testing Chap. 2
Note that when the two sides of (2.15) are equal, then either assignment is equally
good—both have the same effect on the objective function (2.11).
A minor rearrangement of the terms in (2.15) results in
Ĥ(y)=H1
(C01 − C11 ) Pr [H = H1 | y = y] R (C10 − C00 ) Pr [H = H0 | y = y] . (2.16)
Ĥ(y)=H0
Evidently, the a posteriori probabilities, i.e., the probabilities for each of the two hy-
potheses conditioned on having observed the sample value y of the random vector
y, play an important role in the optimum decision rule (2.16). These probabilities
can be readily computed from our measurement models (2.3) together with the
a priori probabilities (2.2). This follows from a simple application of Bayes’ Rule,
viz.,
py|H (y|Hm ) Pm
Pr [H = Hm | y = y] = . (2.17)
py|H (y|H0) P0 + py|H (y|H1 ) P1
Since for any reasonable choice of cost function the cost of an error is higher
than the cost of being correct, the terms in parentheses in (2.16) are both nonnega-
tive, so we can equivalently write (2.16) in the form3
Ĥ(y)=H1
Pr [H = H1 | y = y] (C10 − C00 )
R . (2.18)
Pr [H = H0 | y = y] Ĥ(y)=H0
(C01 − C11 )
3
Technically, we have to be careful about dividing by zero here. To simplify our exposition,
however, as we discuss in Section 2.1.2, we will generally restrict our attention to the case where
this does not happen.
Sec. 2.1 Binary Random Hypothesis Testing: A Bayesian Approach 61
When we then substitute (2.17) into (2.18) and multiply both sides by P0 /P1 , we
obtain the decision rule in its final form, directly in terms of the measurement
densities:
py|H (y|H1 ) Ĥ(y)=H1 P0 (C10 − C00 )
L(y) , R , η. (2.19)
py|H (y|H0 ) P1 (C01 − C11 )
Ĥ(y)=H0
The left side of (2.19) is a function of the observed data y referred to as the
likelihood ratio—which we denote using L(y)—and is constructed from the mea-
surement model. The right side of (2.19)—which we denote using η—is a precom-
putable threshold which is determined from the a priori probabilities and costs.
The overall decision rule then takes the form of what is referred to as a likelihood
ratio test (LRT).
Several observations lend valuable insights into the optimum decision rule (2.19).
First, note that the likelihood ratio L(·) is a scalar-valued function, i.e., L : RK → R,
regardless of the dimension K of the data. In fact, L(y) is an example of what is
referred to as a sufficient statistic for the problem: it summarizes everything we
need to know about the observation vector in order to make a decision. Phrased
differently, in terms of our ability to make the optimum decision (in the Bayesian
sense in this case), knowledge of L(y) is as good as knowledge of the full data
vector y itself.
We will develop the notion of a sufficient statistic in more detail in subse-
quent chapters; however, at this point it suffices to make two observations with
respect to our detection problem. First, (2.19) tells us an explicit construction for a
scalar sufficient statistic for the Bayesian binary hypothesis testing problem. Sec-
ond, sufficient statistics are not unique. For example, the data y itself is a sufficient
statistic, albeit a trivial one. More importantly, any invertible function of L(y) is
also a sufficient statistic. In fact, for the purposes of implementation or analysis it
is often more convenient to rewrite the likelihood ratio test in the form
Ĥ(y)=H1
ℓ(y) = g(L(y)) R g(η) = γ, (2.20)
Ĥ(y)=H0
It follows immediately from the definition in (2.19) that the likelihood ratio
is a nonnegative quantity. Furthermore, depending on the problem, some values
of y may lead to L(y) being zero or infinite. In particular, the former occurs when
py|H (y|H1 ) = 0 but py|H (y|H0) > 0, which is an indication that values in a neighbor-
hood of y effectively cannot occur under H1 but can under H0 . In this case, there
will be values of y for which we’ll effectively know with certainty that the correct
hypothesis is H0 . When the likelihood ratio is infinite, corresponding a division
by zero scenario, an analogous situation exists, but with the roles of H0 and H1
reversed. These cases where such perfect decisions are possible are referred to as
singular detection scenarios. In some practical problems, these scenarios do in fact
occur. However, in other cases they suggest a potential lack of robustness in the
data modeling, i.e., that some source of inherent uncertainty may be missing from
the model. In any event, to simplify our development for the remainder of the
chapter we will largely restrict our attention to the case where 0 < L(y) < ∞ for
all y.
While the likelihood ratio focuses the observed data into a single scalar for
the purpose of making an optimum decision, the threshold η for the test plays a
complementary role. In particular, from (2.19) we see that η focuses the relevant
features of the cost function and a priori probabilities into a single scalar. Further-
more, this information is combined in a manner that is intuitively satisfying. For
example, as (2.19) also reflects, an increase in P0 means that H0 is more likely, so
that η is increased to appropriately bias the test toward deciding H0 for any partic-
ular observation. Similarly, an increase in C10 means that deciding H1 when H0 is
true is more costly, so η is increased to appropriately bias the test toward deciding
H0 to offset this risk. Finally, note that adding a constant to the cost function (i.e.,
to all Cij ) has, as we would anticipate, no effect on the threshold. Hence, without
loss of generality we may set at least one of the correct decision costs—i.e., C00 or
C11 —to zero.
Finally, it is important to emphasize that the likelihood ratio test (2.19) indi-
rectly determines the decision regions (2.6). In particular, we have
As Fig. 2.1 suggests, while a decision rule expressed in the measurement data
space {y} can be complicated,4 (2.19) tells us that the observations can be trans-
formed into a one-dimensional space defined via L = L(y) where the decision
regions have a particularly simple form: the decision Ĥ(L) = H0 is made when-
ever L lies to the left of some point on the line, and Ĥ(L) = H1 whenever L lies to
the right.
4
Indeed, the respective sets Z0 and Z1 are not even connected in general, even for the case
K = 1.
Sec. 2.1 Binary Random Hypothesis Testing: A Bayesian Approach 63
An important cost assignment for many problems is that given by (2.10), which as
we recall corresponds to a minimum probability-of-error (Pr(e)) criterion. Indeed,
in this case, we have
h i h i
ϕ(Ĥ) = Pr Ĥ(y) = H0 , H = H1 + Pr Ĥ(y) = H1 , H = H0 = Pr(e).
The corresponding decision rule in this case can be obtained by simply specializ-
ing (2.19) to obtain
Ĥ(y)=H1
py|H (y|H1) P0
R . (2.22)
py|H (y|H0) P1
Ĥ(y)=H0
Ĥ(y)=H1
Pr [H = H1 | y = y] R Pr [H = H0 | y = y] . (2.23)
Ĥ(y)=H0
From (2.23) we see that to minimize the probability of a decision error, we should
choose the hypothesis corresponding to the largest a posteriori probability, i.e.,
For this reason, we refer to the test associated with this cost assignment as the
maximum a posteriori (MAP) decision rule.
Still further simplification is possible when the hypotheses are equally likely
(P0 = P1 = 1/2). In this case, (2.22) becomes simply
Ĥ(y)=H1
py|H (y|H1) R py|H (y|H0 ),
Ĥ(y)=H0
and thus we see that our optimum decision rule chooses the hypothesis for which
the corresponding likelihood function is largest, i.e.,
This special case is referred to as the maximum likelihood (ML) decision rule. Max-
imum likelihood detection plays an important role in a large number of applica-
tions, and in particular is widely used in the design of receivers for digital com-
munication systems.
64 Detection Theory, Decision Theory, and Hypothesis Testing Chap. 2
Example 2.2
Continuing with Example 2.1, we can obtain from (2.5) that the likelihood ratio test
for this problem takes the form
1 2 2
√ e−(y−s1 ) /(2σ ) Ĥ(y)=H1
2
L(y) = 2πσ R η. (2.26)
1 −(y−s0 )2 /(2σ2 )
√ e Ĥ(y)=H0
2πσ 2
As (2.26) suggests—and as is generally the case in Gaussian problems—the natural
logarithm of the likelihood ratio is a more convenient sufficient statistic to work with
in this example. In this case, taking logarithms of both sides of (2.26) yields
Ĥ(y)=H1
1
ℓ(y) = 2 (y − s0 )2 − (y − s1 )2 R
ln η. (2.27)
2σ
Ĥ(y)=H0
Expanding the quadratics and cancelling terms in (2.27) we obtain the test in its
simplest form, which for s1 > s0 is given by
Ĥ(y)=H1
s1 + s0 σ 2 ln η
y R + , γ. (2.28)
2 s1 − s0
Ĥ(y)=H0
As we’ll see later in the chapter, this minimum-distance property holds in multidi-
mensional Gaussian problems as well.
Note too that in this problem the decisions regions on the y-axis have a partic-
ularly simple form; for example, for s1 > s0 we obtain
Z0 = {y | y < γ}
(2.29)
Z1 = {y | y > γ}.
In other problems—even Gaussian ones—the decision regions can be more compli-
cated, as our next example illustrates.
Example 2.3
Suppose that a zero-mean Gaussian random variable has one of two possible vari-
ances, σ12 or σ02 , where σ12 > σ02 . Let the costs and prior probabilities be arbitrary.
Then the likelihood ratio test for this problem takes the form
1 2 2
p
2
e−y /(2σ1 ) Ĥ(y)=H1
2πσ1
L(y) = R η.
1 −y 2 /(2σ02 )
p e Ĥ(y)=H0
2πσ02
Sec. 2.1 Binary Random Hypothesis Testing: A Bayesian Approach 65
In this problem, it is a straightforward exercise to show that the test simplifies to one
of the form s
Ĥ(y)=H1
σ02 σ12
σ1
|y| R 2 2 ln η , γ.
σ1 − σ02 σ0
Ĥ(y)=H0
Hence, the decision region Z1 is the union of two disconnected regions in this case,
i.e.,
Z1 = {y | y > γ} ∪ {y | y < −γ}.
where Z0 and Z1 are the decision regions defined via (2.6). Using terminology
that originated in the radar community where H1 refers to the presence of a tar-
get and H0 the absence, the quantities PD and PF are generally referred to as the
“detection” and “false-alarm” probabilities, respectively (and, hence, the choice of
notation). In the statistics community, by contrast, PF is referred to as the size of
the test and PD as the power of the test.
It is worth emphasizing that the characterization in terms of (PD , PF ) is not
unique, however. For example, any invertible linear or affine transformation of
the pair (PD , PF ) is also complete. For instance, the pair of “probabilities of error
of the first and second kind” defined respectively via
h i
PE1 = Pr Ĥ(y) = H1 H = H0 = PF
h i (2.31)
PE2 = Pr Ĥ(y) = H0 H = H1 = 1 − PD , PM
5
Note that the arbitrary decision rule we consider here need not be optimized with respect
to any particular criterion—it might be, but it might also be a heuristically reasonable rule, or even
a bad rule.
66 Detection Theory, Decision Theory, and Hypothesis Testing Chap. 2
η→0
1
η
sing
rea
inc
←
PD
η→ ∞
0 Figure 2.2. Operating characteristic
0 PF 1 associated with a likelihood ratio test.
probability of error of the first kind is the probability of a false alarm, while prob-
ability of error of the second kind is the probability of a miss, which is denoted by
PM .
In general, a good decision rule (detector) is one with a large PD and a small
PF (or equivalently, small PE1 and PE2 ). However, ultimately these are competing
objectives. As an illustration of this behavior, let us examine the performance of a
likelihood ratio test (2.19) when the threshold η is varied. Note that each choice of η
completely specifies a decision rule, with which is associated a particular (PD , PF )
operating point. Hence, each value of η is associated with a single point in the PD –
PF plane. Moreover, as η is varied from 0 to ∞, a curve is traced out in this plane
as illustrated in Fig. 2.2. This curve is referred to as the operating characteristic of
the likelihood ratio test.
As Fig. 2.2 suggests, good PD is generally obtained at the expense of high
PF , and so choosing a threshold η for a particular problem involves making an
acceptable tradeoff. Indeed, as η → 0 we have (PD , PF ) → (1, 1), while as η → ∞
we have (PD , PF ) → (0, 0). From this perspective, the Bayesian test represents a
particular tradeoff, and corresponds to a single point on this curve. To obtain this
tradeoff, we effectively selected as our objective function a linear combination of
PD and PF . More specifically, we performed the optimization (2.8) using
ϕ(f ) = αPF − βPD + γ,
where the choice of α and β is, in turn, determined by the cost assignment (Cij ’s)
and the a priori probabilities (Pm ’s). In particular, rewriting (2.9) in the form
X h i
ϕ(f ) = Cij Pr Ĥ(y) = Hi H = Hj Pj
i,j
Sec. 2.1 Binary Random Hypothesis Testing: A Bayesian Approach 67
we obtain
α = (C10 − C00 )P0 β = (C01 − C11 )P1 γ = (C00 P0 + C01 P1 ).
Example 2.4
Let us consider the following special case of our simple scalar Gaussian detection
problem from Example 2.1:
H0 : y ∼ N (0, σ 2 )
(2.32)
H1 : y ∼ N (m, σ 2 ), m ≥ 0,
so that
Z ∞
PD = py |H (y|H1 ) dy (2.33a)
γ
Z ∞
PF = py |H (y|H0 ) dy. (2.33b)
γ
6
For any reasonable performance criterion, a curve above and to the left of another is always
preferable.
68 Detection Theory, Decision Theory, and Hypothesis Testing Chap. 2
0.9
0.7
0.6
0.5
0.4
0.3
0.2
Figure 2.3. Operating characteristic of
0.1 the likelihood ratio test for the scalar
0
Gaussian detection problem. The suc-
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 cessively higher curves correspond to
False−Alarm Probability PF
d = m/σ = 0, 1, . . . , 5.
Example 2.5
As motivation, let us reconsider our simple radar or communications scenario (2.32)
of Example 2.4. For this problem we showed that the decision rule that minimizes
the Bayes risk for a cost assignment {Cij } and a set of a priori probabilities {Pm }
reduces to
Ĥ(y)=H1
m σ2
(C10 − C00 )P0
y R + ln . (2.36)
2 m (C01 − C11 )P1
Ĥ(y)=H0
For scenarios such as that in Example 2.5, let us assess the impact on per-
formance of using a test optimized for the a priori probabilities {1 − p, p} when
the corresponding true a priori probabilities are {P0 , P1 }. The Bayes risk for this
“mismatched” system is
ϕ(p, P1 ) = C00 (1−P1 )+C01 P1 +(C10 −C00 )(1−P1 )PF (p)−(C01 −C11 )P1 PD (p), (2.37)
where we have explicitly included the dependence of PF and PD on p to empha-
size that these conditional probabilities are determined from a likelihood ratio test
whose threshold is computed using the incorrect prior p.
The system performance in this situation has a convenient geometrical inter-
pretation, as we now develop. With the notation (2.37), ϕ(P1 , P1 ) denotes the Bayes
risk when the likelihood ratio test corresponding to the correct priors is used. By
the optimality properties of this latter test, we know that
ϕ(p, P1 ) ≥ ϕ(P1 , P1 ). (2.38)
with, of course, equality if p = P1 . Moreover, for a fixed p, the mismatch Bayes
risk ϕ(p, P1 ) is a linear function of the true prior P1 . Hence, we can conclude, as
depicted in Fig. 2.4, that when plotted as functions of P1 for a particular p, the
mismatch risk ϕ(p, P1) is a line tangent to ϕ(P1 , P1 ) at P1 = p. Furthermore, since
70 Detection Theory, Decision Theory, and Hypothesis Testing Chap. 2
ϕ(p, P1)
ϕ ϕ(P1, P1 )
C11
Figure 2.4. Bayes risk objective func-
tions for min-max tests. The solid
C00 curve indicates the objective function
as a function of the true prior probabil-
ity for a correctly matched detector—
i.e., the detector has correct knowl-
0 p 1 edge of the prior probabilities. The
dashed curve indicates the objective
P1 function with mismatched detector; for
this curve the assumed prior is p.
this must be true for all choices of p, we can further conclude that the optimum
risk ϕ(P1 , P1 ) must be a concave function of P1 as Fig. 2.4 also reflects.
From the geometric picture of Fig. 2.4, it is clear that the performance ϕ(p, P1 )
obtained using a fixed assumed prior p when the correct prior is P1 , depends on
the value of P1 . Moreover, because the mismatch risk is a linear function, we see
that the poorest performance, corresponding to worst case mismatch, is obtained
either when P1 = 0 or P1 = 1, depending on the sign of the slope of the mismatch
risk function. For example, for the p shown in Fig. 2.4, this worst case performance
takes place when the true prior is P1 = 1.
Given this behavior, a conservative design strategy is to use a likelihood ratio
test based on an assumed prior p chosen so that the worst-case performance is as
good as possible. Mathematically, this corresponds to choosing the assumed prior
according to a “min-max” criterion: the prior P̂1 obtained in this manner is given
by
n o
P̂1 = arg min max ϕ(p, P1) . (2.39)
p P1
This minimizes the sensitivity of ϕ(p, P1) to variations in P1 , and hence leads to a
decision rule that is robust with respect to uncertainty in the prior P1 .
The solution to this min-max problem follows readily from the geometrical
picture of Fig. 2.4. Our solution depends on on the details of the shape of ϕ(P1 , P1 )
over the range 0 ≤ P1 ≤ 1. There are three cases, which we consider separately.
An example of this case is depicted in Fig. 2.5(a). In this case, the slope of any
line tangent to ϕ(P1 , P1 ) is nonpositive, so for any p the maximum of ϕ(p, P1 ) lies
Sec. 2.2 Min-Max Hypothesis Testing 71
An example of this case is depicted in Fig. 2.5(b). In this case, the slope of any line
tangent to ϕ(P1 , P1 ) is nonnegative, so for any p the maximum of ϕ(p, P1 ) lies at
P1 = 1. Hence, as the solution to (2.39) we obtain P̂1 = 1. This corresponds to
using a likelihood ratio test with threshold
1 − P̂1 (C10 − C00 )
η̂ = = 0.
P̂1 (C01 − C11 )
This detector makes the decision Ĥ = H1 regardless of the observed data, and
therefore corresponds to the operating point (PD , PF ) = (1, 1).
An example of this case is depicted in our original Fig. 2.4. In this case, ϕ(P1 , P1 )
has a point of zero slope (and hence a maximum) at an interior point (0 < P1 < 1),
so that P̂1 is the value of p for which the slope of ϕ(p, P1 ) is zero.
The corresponding point on the operating characteristic (and in turn, implic-
itly, the threshold η̂) can be determined geometrically. In particular, substituting
(2.37) with our zero-slope condition, it follows that P̂1 satisfies
d
ϕ(p, P1 ) = (C01 − C00 ) − (C10 − C00 )PF (P̂1 ) − (C01 − C11 )PD (P̂1 ) = 0. (2.40)
dP1 p=P̂1
C00 ϕ(P1, P1 )
C11
0 P1 1
(a)
ϕ(P1, P1 ) C11
C00
0 P1 1
(b) Figure 2.5. Examples of possible
ϕ(P1 , P1 ) curves without points of zero
slope.
Sec. 2.3 Neyman-Pearson Binary Hypothesis Testing 73
As a final comment, examining the endpoints of the curve in Fig. 2.4 we see
that for some cost assignments we can guarantee that this Case 3 will apply. For
example, this happens when C00 = C11 = 0. In this case, we see that the particular
operating point on the likelihood ratio test operating characteristic is defined as
the point for which the cost of a miss is the same as the cost of a false alarm, i.e.,
C01 (1 − PD ) = C10 PF . (2.42)
As C01 and C10 are varied relative to one another, the likelihood ratio test threshold
η̂ is varied accordingly.
Both our basic Bayesian and min-max hypothesis testing formulations require that
we choose suitable cost assignments Cij . As we saw, these cost assignments di-
rectly influence the (PD , PF ) operating point of the optimum decision rule. How-
ever, in many applications there is no obvious set of cost assignments.
In this kind of situation, an optimization criterion that is frequently more
natural is to choose the decision rule so as to maximize PD subject to a constraint
on the maximum allowable PF , i.e.,
max PD such that PF ≤ α.
Ĥ(·)
Now recall from our earlier discussion that specifying Z0 fully determines
Ĥ(·), so we can view our problem as one of determining the optimum Z0 . From
this perspective it is clear we want to choose Z0 so that it contains precisely those
values of y for which the term in brackets inside the integral in (2.43) is negative,
since this choice makes ϕ(Ĥ) smallest. This statement can be expressed in the form
Ĥ(y)=H1
py|H (y|H1) − λpy|H (y|H0 ) R 0,
Ĥ(y)=H0
A large number of detection and decision problems take the form of binary hy-
pothesis tests in which the data are jointly Gaussian under each hypothesis. In
this section, we explore some of the particular properties of likelihood ratio tests
for these problems. The general scenario we consider takes the form
H0 : y ∼ N(m0 , Λ0 )
(2.45)
H1 : y ∼ N(m1 , Λ1 ).
Sec. 2.4 Gaussian Hypothesis Testing 75
Note that y could be a random vector obtained from some array of sensors, or it
could be a collection of samples obtained from a random discrete-time signal y [n],
e.g., T
y = y [0] y [1] · · · y [N − 1] . (2.46)
For the hypotheses (2.45), and provided Λ0 and Λ1 are nonsingular, the like-
lihood ratio test takes the form
1 T −1
1
(2π)N/2 |Λ1 | 1/2 exp − 2
(y − m1 ) Λ 1 (y − m1 ) Ĥ(y)=H1
L(y) = R η. (2.47)
1
exp − 12 (y − m0 )T Λ−1
(2π)N/2 |Λ0 |1/2 0 (y − m0 ) Ĥ(y)=H 0
Defining T
∆m = ∆m[0] ∆m[1] · · · ∆m[K − 1] = m1 − m0 , (2.52)
we see from simplifying (2.48) that a sufficient statistic for making the optimum
Bayesian decision at the detector (regardless of the cost and a priori probability as-
signments) is
K−1
X
ℓ(y) = yT ∆m = y[n] ∆m[n], (2.53)
n=0
76 Detection Theory, Decision Theory, and Hypothesis Testing Chap. 2
Note that the correlation computation (2.53) that defines ℓ(y) is equivalent to
a “convolution and sample” operation. Specifically, it is a straightforward exercise
to verify that ℓ can be expressed in the form
or
Ĥ(y)=H1
1
ℓ′ (y) , (m1 − m0 )T Λ−1 y R 2 ln η + mT −1 T −1 ′
1 Λ m1 − m0 Λ m0 , η . (2.59)
2
Ĥ(y)=H0
7
We leave it as an exercise for you to verify that we have defined a valid inner product.
Sec. 2.4 Gaussian Hypothesis Testing 77
where our distance metric is defined via our norm (2.61). In addition, via (2.59),
when km0 k = km1 k, solving for the minimum distance is equivalent to solving for
the maximum projection; in particular, specializing (2.59) we obtain
Ĥ(y) = Hm̂ where m̂ = arg max hy, mm i ,
m∈{0,1}
where our inner product is that defined in (2.62). In our vector space {y}, the
corresponding decisions regions are separated by a hyperplane equidistant from
m0 and m1 and perpendicular to the line connecting them.
The performance of these likelihood ratio tests can be readily calculated since
′
ℓ (y) is a linear function of the Gaussian random vector y under each hypothesis.
As a result, ℓ′ (y) is Gaussian under each hypothesis. Specifically,
H0 : ℓ′ ∼ N(m′0 , σℓ2′ )
(2.63)
H1 : ℓ′ ∼ N(m′1 , σℓ2′ ),
where
PD = Pr [ℓ′ ≥ η ′ | H = H1 ]
′
ℓ − m′1 η ′ − m′1 η ′ − m′1
= Pr ≥ H = H1 = Q (2.66a)
σℓ′ σℓ′ σℓ′
PF = Pr [ℓ′ ≥ η ′ | H = H0 ]
′
ℓ − m′0 η ′ − m′0 η ′ − m′0
= Pr ≥ H = H0 = Q . (2.66b)
σℓ′ σℓ′ σℓ′
While computing actual values of PD and PF requires that we evaluate the Q (·)
function numerically, we can obtain bounds on these quantities that are explic-
itly computable. For example, using (1.129) and (1.128) with (2.66) we obtain, for
78 Detection Theory, Decision Theory, and Hypothesis Testing Chap. 2
As in Example 2.6, the sufficient statistic for this problem can be imple-
mented either directly as a correlator, or as a matched filter. In the latter case, the
filter h[n] captures the relevant information about the signal and the noise, taking
the form
0
n≥K
h[n] = ∆m′ [K − 1 − n] 0 ≤ n ≤ K − 1 (2.68)
0 n ≤ −1,
where
T
∆m′ = ∆m′ [0] ∆m′ [1] · · · ∆m′ [K − 1] = Λ−1 (m1 − m0 ).
(2.69)
Example 2.7
Suppose that we are trying to detect a random signal x[n] at an antenna based on
observations of the samples n = 0, 1, . . . K − 1. Depending on whether the signal is
absent or present, the noisy observations y [n] take the form
H0 : y [n] = w [n]
(2.70)
H1 : y [n] = x[n] + w [n],
where under both hypotheses w [n] is an IID sequence of N (0, σ 2 ) random variables
that is independent of the sequence x[n].
If x[n] was a known signal x[n], then our optimum detector would be that de-
termined in Example 2.6 with ∆m[n] = x[n]; i.e., our sufficient statistic is obtained
by correlating our received signal y[n] with x[n]. However, in this example we as-
sume we know only the statistics of x[n]—specifically, we know that
T
x = x[0] x[1] · · · x[K − 1] ∼ N (0, Λx )
where Λx is the signal covariance matrix.
In this case, we have that the likelihood ratio test (2.48) simplifies to the test
!
Ĥ(y)=H1
T 2 |σ 2 I + Λx |1/2
ℓ(y) , y x̂(y) R 2σ ln η ,γ (2.71)
σK
Ĥ(y)=H0
where −1 T
x̂(y) = Λx σ 2 I + Λx
y = x̂[0] x̂[1] · · · x̂[K − 1] . (2.72)
Sec. 2.5 Tests with Discrete-Valued Observations 79
Thus far we have focussed in this chapter on the case of observation vectors y
that are specifically continuous-valued. However, there are many important deci-
sion and detection problems that involve inherently discrete-valued observations.
Conveniently, the preceding theory carries over to the discrete case in a largely
straightforward manner, with the likelihood ratio test continuing to play a central
role, as we now develop. However, there are at least some important differences,
which we will emphasize.
With the notation (2.4), the optimum Bayesian decision rule takes a form analo-
gous to that for the case of continuous-valued data; specifically,
Ĥ(y)=H1
py|H [y|H1 ] P0 (C10 − C00 )
L(y) , R , η. (2.73)
py|H [y|H0 ] Ĥ(y)=H0
P1 (C01 − C11 )
The derivation of this result closely mimics that for the case of continuous-valued
data, the verification of which we leave as an exercise for the reader.
80 Detection Theory, Decision Theory, and Hypothesis Testing Chap. 2
A couple of special issues arise in the case of likelihood ratio tests of the
form (2.73). In particular, when the observations are discrete-valued, the like-
lihood function L(y) is also discrete-valued; for future convenience, we denote
these values by 0 ≤ η0 < η1 < η2 < · · · . This property has some significant
consequences. First, it means that for many sets of costs and a priori probability
assignments, the resulting η formed in (2.73) will often not coincide with one of
the possible values of L = L(y), i.e., we will often have η 6= ηi for all i. In such
cases, the case of equality in the likelihood ratio test will not arise, and the min-
imum Bayes risk is achieved by a likelihood ratio test that corresponds, as usual,
to a unique (PD , PF ) point on the associated operating characteristic.
However, for some choices of the costs and a priori probabilities, the resulting
threshold η will satisfy η = ηi for some particular i. This means that, unlike with
continuous-valued data, in this case equality in the likelihood ratio test will occur
with nonzero probability. Nevertheless, it can be readily verified that the Bayes
risk is the same no matter how a decision is made in this event. Hence, when
equality occurs, the decision can still be made arbitrarily, but these choices will
correspond to different points on the operating characteristic.
2.5.2 The Operating Characteristic of the Likelihood Ratio Test, and Neyman-
Pearson Tests
Let us more generally examine the form of the operating characteristic associated
with the likelihood ratio test in the case of discrete-valued observations. In partic-
ular, let us begin by sweeping the threshold η in (2.73) from 0 to ∞ and examining
the (PD , PF ) values that are obtained. As our development in the last section re-
vealed, in order to ensure that each threshold η maps to a unique (PD , PF ), we
need to choose an arbitrary but fixed convention for handling the case of equality
in (2.73). For this purpose let us associate equality with the decision Ĥ(y) = H1 ,
expressing the likelihood ratio test in the form
Ĥ(y)=H1
py|H [y|H1 ] ≥
L(y) = η. (2.74)
py|H [y|H0 ] <
Ĥ(y)=H0
Example 2.8
Suppose in an optical communication system bits are signaled by turning on and
off a laser. At the detector, the measured photon arrival rate is used to determine
whether the laser is on or off. When the laser is off (0-bit), the photons arrive accord-
ing to a Poisson process with average arrival rate m0 ; when the laser is on (1-bit),
the rate is m1 with m1 > m0 . Suppose during a bit period we count the number
Sec. 2.5 Tests with Discrete-Valued Observations 81
of photons y that arrive and use this observed data to make a decision. Then the
likelihood functions for this decision problem are
y −mi
mi e y = 0, 1, . . .
py |H [y|Hi ] = Pr [y = y | H = Hi ] = y! .
0 otherwise
y Ĥ(y)=H1
m1 −(m1 −m0 ) ≥
L(y) = e η,
m0 <
Ĥ(y)=H0
The discrete nature of the this hypothesis testing problem means that the op-
erating characteristic associated with the likelihood ratio test (2.75) is a discrete col-
lection of points rather than a continuous curve of the type we encountered in an
earlier example involving Gaussian data. Indeed, while the left-hand side of (2.75)
is integer-valued, the right side is not in general. As a result, we have that PD and
PF are given in terms of γ by the expressions
X my e−m1
1
PD = Pr [y ≥ γ | H = H1 ] = Pr [y ≥ ⌈γ⌉ | H = H1 ] = (2.76a)
y!
y≥⌈γ⌉
X my e−m0
0
PF = Pr [y ≥ γ | H = H0 ] = Pr [y ≥ ⌈γ⌉ | H = H0 ] = . (2.76b)
y!
y≥⌈γ⌉
The resulting operating characteristic is depicted in Fig. 2.6. In Figure 2.6, only
the isolated (PD , PF ) points indicated by circles are achievable by simple likelihood
ratio tests. For example, the uppermost point is achieved for all γ ≤ 0, the next
highest point for all 0 < γ ≤ 1, the next for 1 < γ ≤ 2, and so on. This behavior is
representative of such discrete decision problems.
As we discussed in Section 2.5.1, for Bayesian problems the specific cost and
a priori probability assignments determine a threshold η in (2.74), which in turn
corresponds to one of the isolated points comprising the operating characteristic.
For Neyman-Pearson problems with discrete observations, it is also straight-
forward to show that a likelihood ratio test of the form (2.74) is the optimum de-
terministic decision rule. Moreover, as we might expect, for this rule the threshold
η is chosen so as to achieve the largest PD subject to our contraint on the maximum
allowable PF (i.e., α). When α corresponds to at least one of the discrete points of
the operating characteristic, then the corresponding PD indicates the achievable
detection probability. More typically, however, α will lie strictly between the PF
82 Detection Theory, Decision Theory, and Hypothesis Testing Chap. 2
0.9
0.7
0.6
0.5
0.4
values of points on the operating characteristic. In this case, the appropriate oper-
ating point corresponds to the (PD , PF ) whose PF is the largest PF that is smaller
than α. This operating point then uniquely specifies the decision rule.
While this decision rule is a straightforward extension of the result for conti-
nuous-valued observations, it is possible to develop a more sophisticated decision
rule that typically yields at least somewhat better performance as measured by
the Neyman-Pearson criterion. To see this, we first note that the straightforward
likelihood ratio test defined above was obtained as the optimum decision rule
among all possible deterministic tests. Specifically, the likelihood ratio test was the
best decision rule of the form
(
H0 y ∈ Z0
Ĥ(y) =
H1 y ∈ Z1
where Z0 and Z1 are mutually exclusive, collectively exhaustive sets in the obser-
vation space {y}. For these rules, each observation y maps to a unique decision
Ĥ(y).
When the regular likelihood ratio test cannot meet the false-alarm constraint
with equality, better performance can be achieved by employing a decision rule
chosen from outside the class of deterministic rules. In particular, if we allow for
some randomness in the decision process, so that a particular observation does not
always produce the same decision, we can obtain improved detection performance
while meeting our false-alarm constraint. This can be accomplished as follows.
Consider the sequence of thresholds η0 , η1 , . . . that correspond to values that the
likelihood ratio can take on, and denote the corresponding operating points by
(PD (ηi ), PF (ηi )) for i = 0, 1, . . . . Determine ı̂ such that ηı̂ is the threshold value that
results in the likelihood ratio test with the smallest false-alarm probability that is
greater than α. Then as illustrated in Fig. 2.7, PF (ηı̂ ) and PF (ηı̂+1 ) “bracket” α:
using ηı̂+1 results in test with the largest false-alarm probability that is less than α.
Sec. 2.5 Tests with Discrete-Valued Observations 83
PD
ηî
η î+1
Now consider the following randomized decision rule. We flip a biased coin for
which the probability of “heads” is p and that of “tails” is 1−p. If “heads” turns
up, we use the likelihood ratio test with threshold ηı̂ ; if “tails” turns up, we use
the likelihood ratio test with threshold ηı̂+1 . Then the resulting test achieves
and corresponds to a point on the line segment connecting (PD (ηı̂ ), PF (ηı̂ )) and
(PD (ηı̂+1 ), PF (ηı̂+1 )). This line segment is indicated with the dashed line segment
in Fig. 2.7. In particular, as p is varied from 0 to 1, the operating point moves from
(PD (ηı̂+1 ), PF (ηı̂+1 )) to (PD (ηı̂ ), PF (ηı̂ )). Hence, by choosing p appropriately, we
can achieve the false-alarm probability constraint with equality, i.e., there exists a
p̂ such that
PF = p̂PF (ηı̂ ) + (1 − p̂)PF (ηı̂+1 ) = α.
Our results in this section have suggested that at least in some problems
involving discrete-valued data and the Neyman-Pearson criterion a randomized
test can lead to better performance than a deterministic test. This observation, in
turn, raises several natural and important questions. For example, we considered
a very particular class of randomized tests—specifically, a simple random choice
between the outcomes of two (deterministic) likelihood ratio tests. Would some
other type of randomized test be able to perform better still? And could some
more general form of randomized test be able to improve performance in the case
of continuous-valued data with Bayesian or Neyman-Pearson criteria? To answer
these questions, in the next section we develop decision rules optimized over of a
broad class of randomized tests.
For a randomized test, the decision rule is a random function of the data, which
we denote using Ĥ(·). Hence, even for a deterministic argument y, the decision
Ĥ(y) is a random quantity. However, Ĥ(y) has the property that conditioned on
knowledge of y, the function is independent of the hypothesis H. Such a test is
fully described by the probabilities
h i h i
Q0 (y) = Pr Ĥ(y) = H0 y = y = Pr Ĥ(y) = H0 y = y, H = Hi
h i h i (2.79)
Q1 (y) = Pr Ĥ(y) = H1 y = y = Pr Ĥ(y) = H1 y = y, H = Hi ,
With this notation, we see that deterministic rules are a special case, corre-
sponding to
(
1 y ∈ Z1
Q1 (y) = . (2.80)
0 y ∈ Z0
Moreover, it also follows immediately that tests formed by a random choice among
two likelihood ratio tests, such as were considered in Section 2.5.2 are also special
cases. For example, the test described via (2.78) corresponds to
1 L(y) ≥ ηı̂+1
Q1 (y) = p̂ L(y) = ηı̂ (2.81)
0 L(y) ≤ ηı̂−1 .
More generally, for a randomized test that corresponds to the random choice be-
tween two likelihood ratio tests with respective thresholds η1 and η2 such that
η2 > η1 , and where the probability of selecting the first test is p, it is straightfor-
ward to verify that
1 L(y) ≥ η2
Q1 (y) = p η1 ≤ L(y) < η2 (2.82)
0 L(y) < η1 .
We begin by establishing a more general version of our Bayesian result for the
case of randomized tests. We consider the case of continuous-valued data; the
derivation in the discrete case is analogous. We begin by writing our Bayes risk in
the form Z
ϕ(Q0 ) = ϕ̃(y) py (y) dy,
where h i
ϕ̃(y) = E C(H, Ĥ(y)) y = y .
86 Detection Theory, Decision Theory, and Hypothesis Testing Chap. 2
Again we see it suffices to minimize ϕ̃(y) for each y. Applying, in turn, (2.79) and
Bayes’ Rule, we can write ϕ̃(y) in the form
X h i
ϕ̃(y) = Cij Pr H = Hj , Ĥ(y) = Hi y = y
i,j
X h i
= Cij Pr [H = Hj | y = y] Pr Ĥ = Hi y = y
i,j
X Pj py|H (y|Hj )
= Cij Qi (y) , (2.83)
i,j
py (y)
We next establish a more general version of our Neyman-Pearson result for the
case of randomized tests. We begin with the case of continuous-valued data. As in
the deterministic case, we follow a Lagrange multiplier approach, expressing our
objective function as
ϕ(Q0 ) = 1 − PD + λ(PF − α′ )
h i h i
= λ(1 − α′ ) + Pr Ĥ(y) = H0 H = H1 − λ Pr Ĥ(y) = H0 H = H0
(2.86)
for some α′ ≤ α. This time, though, (2.86) expands as
Z h i
′
ϕ(Q0 ) = λ(1 − α ) + Pr Ĥ(y) = H0 y = y py|H (y|H1 ) − λpy|H (y|H0 ) dy
Sec. 2.6 Randomized Tests 87
Hence, we can conclude that our optimum rule is almost a simple (deterministic)
likelihood ratio test. In particular, at least except when L(y) = λ we have that
(
H1 L(y) > λ
Ĥ(y) = (2.89)
H0 L(y) < λ.
It remains only to determine what the nature of the decision (i.e., Q0 (y))
when L(y) = λ and the choice of λ. These quantities are determined by meet-
ing the constraint PF = α′ . Specifically,
Z
PF = [1 − Q0 (y)] py|H (y|H0) dy
Z Z
= py|H (y|H0) dy + [1 − Q0 (y)] py|H (y|H0 ) dy
{y|L(y)>λ} {y|L(y)=λ}
Z
= Pr [L(y) > λ | H = H0 ] + [1 − Q0 (y)] py|H (y|H0) dy, (2.90)
{y|L(y)=λ}
where the second equality follows from substituting for Q0 (y) using (2.88). Ob-
serve that PF = 1 if λ < 0, so it suffices to restrict our attention to λ ≥ 0. Further-
more, note that the first term in (2.90) is a nonincreasing function of λ.
Now since y is a continuous-valued random variable, then the second term in
(2.90) is zero and the remaining term, which is a continuous of λ, can be chosen so
that it equals α′ . In this case, it does not matter how we choose Q0 (y) when L(y) =
λ, so the optimum randomized rule degenerates to a deterministic likelihood ratio
test again. It remains only to show that for optimum PD we want α′ = α. Given
our preceding results, it suffices to exploit the fact that for likelihood ratio tests PD
is a monotonically nondecreasing function of PF .
88 Detection Theory, Decision Theory, and Hypothesis Testing Chap. 2
Discrete Case
When the data is discrete, some important differences arise, which we now ex-
plore. Proceeding as we did at the outset of Section 2.6.2, we obtain
X
ϕ(Q0 ) = λ(1 − α′ ) +
Q0 (y) py|H [y|H1 ] − λpy|H [y|H0 ]
y
X
= λ(1 − α′ ) +
Q0 (y) L(y) − λ py|H [y|H0 ]. (2.91)
y
from which we analogously conclude that for L(y) 6= λ we have (2.88). Hence, it
remains only to determine λ and the decision when L(y) = λ from the false alarm
constraint.
An important difference from the continuous case is that when y is discrete-
valued, so that L(y) is also discrete-valued, the first term in (2.90) is not a continu-
ous function of λ, but is piecewise constant.
As before, let us denote the values that L(y) takes on by η0 , η1 , η2 , . . . , where
0 < η0 < η1 < η2 < · · · . And let us choose ı̂ so that λ = ηı̂ is the smallest threshold
such that
Pr [L(y) > λ | H = H0 ] = Pr [L(y) ≥ ηı̂+1 | H = H0 ] ≤ α′ .
Then we can obtain PF = α′ by choosing
Q1 (y) = 1 − Q0 (y)
appropriately for all y such that L(y) = ηı̂ . In particular, 0 < Q1 (·) < 1 in this
range must be chosen so that
X
α′ − Pr [L(y) ≥ ηı̂+1 | H = H0 ] = Q1 (y) py|H [y|H0 ]
{y|L(y)=ηı̂ }
More generally, when q > 0, we have that the decision rule involves a random
choice between the decision rule (2.94) and
Ĥ(y)=H1
≥
L(y) ηı̂ . (2.95)
<
Ĥ(y)=H0
Sec. 2.7 Properties of the Likelihood Ratio Test Operating Characteristic 89
In particular, the test (2.95) is chosen with probability q, and the test (2.94) with
probability 1 − q. In terms of the operating characteristic of the likelihood ratio
test, this is a point on the line segment connecting the points corresponding to the
two deterministic tests (2.94) and (2.95). This, of course, is precisely the form of
the heuristically designed randomized test we explored at the end of Section 2.5.2,
which we now see is optimal. As a result, when we include randomized tests,
it makes sense to redefine the operating characteristic for the discrete case as the
isolate likelihood ratio test operating points together with the line segments that
connect these discrete points.
It remains only to verify that for optimum PD we want α′ = α in the discrete
case as well. However, for likelihood ratio tests involving discrete data, the PD
values form a nondecreasing sequence as a function of PF . Then since the perfor-
mance of randomized tests corresponds to points on the line segments connecting
the performance points associated with deterministic tests, PD is a nondecreasing
function of PF for our more general class of randomized tests as well.
In summary, optimum Neyman-Pearson decision rules always take the form
of either a deterministic rule in the form of a likelihood ratio test or a randomized
rule in the form of a simple randomization between two likelihood ratio tests. In
both cases, the likelihood ratio test and its associated operating characteristic play
a central role. Accordingly, we explore their properties further.
By exploiting the special role that the likelihood ratio test plays in both determin-
istic and randomized optimum decisions rules, we can develop a number of key
properties of the PD –PF operating characteristic associated with the likelihood ra-
tio test. For future reference, recall that for continuous-valued data the test takes
the form
py|H (y|H1 ) Ĥ(y)=H1
L(y) = R η, (2.96)
py|H (y|H0 )
Ĥ(y)=H0
We emphasize at the outset that the detailed shape of the operating charac-
teristic is determined by the measurement model for the data—for example, by
py|H (y|H0 ) and py|H (y|H1) in the continuous-case—since it is this information that
is used to construct the likelihood ratio L(y). However, all operating character-
istics share some important characteristics in common, and it is these that we ex-
plore in this section. As a simple example, which was mentioned earlier, we have
90 Detection Theory, Decision Theory, and Hypothesis Testing Chap. 2
that the (PD , PF ) points (0, 0) and (1, 1) always lie on the operating characteristic,
and correspond to η → ∞ and η → 0, respectively.
It is also straightforward to verify that PD ≥ PF , i.e., that the operating char-
acteristic always lies above the diagonal in the PD –PF plane. This can be verified
using a randomization argument. In particular, suppose our decision rule ignores
the data y and bases its decision solely on the outcome of a biased coin flip, where
the probability of “heads” is p. If the coin comes up “heads” we make the decision
Ĥ(y) = H1 , while if it comes up “tails” we make the decision Ĥ(y) = H0 . Then for
this rule we have
h i h i
PD = Pr Ĥ(y) = H1 H = H1 = Pr Ĥ(y) = H1 = p
h i h i
PF = Pr Ĥ(y) = H1 H = H0 = Pr Ĥ(y) = H1 = p,
8
Indeed, we wouldn’t expect to obtain better performance by ignoring the data y!
Sec. 2.7 Properties of the Likelihood Ratio Test Operating Characteristic 91
show that at those points where it is defined, the slope of the operating character-
istic is numerically equal to the corresponding threshold η, i.e.,
dPD
= η. (2.99)
dPF
Note that since η ≥ 0, this is another way of verifying that the operating charac-
teristic is nondecreasing. A proof is as follows. With
Z1 (η) = {y | L(y) > η} (2.100)
we have Z
PD (η) = py|H (y|H1 ) dy,
Z1 (η)
which after applying the definition of the likelihood function (2.96) yields
Z
PD (η) = L(y) py|H (y|H0) dy. (2.101)
Z1 (η)
Next, note that with u(·) denoting the unit step function, i.e.,
(
1 x>0
u(x) = ,
0 otherwise
we have (
1 y ∈ Z1 (η)
u(L(y) − η) = . (2.102)
0 otherwise
we know
dPF
= −pL|H (η|H0). (2.105)
dη
Hence, dividing (2.104) by (2.105) we obtain (2.99).
92 Detection Theory, Decision Theory, and Hypothesis Testing Chap. 2
More generally, if Q0 (·) describes a randomized decision rule, then the correspond-
ing decision rule, which we denote using Q0 (·), is defined via
Q0 (y) = Q1 (y) = 1 − Q0 (y).
Using this property, it follows that those points lying on the curve correspond-
ing to the operating characteristic reflected across the PD = 1/2 and PF = 1/2
lines are achievable by likelihood ratio tests whose decisions are reversed. In turn,
all points between the diagonal and this “reflected operating characteristic” are
Sec. 2.8 M -ary Hypothesis Testing 93
PD
0
0 PF 1 Figure 2.8. Achievable region in the
PD –PF plane.
Thus far we have focussed on the case of binary hypothesis testing in this chap-
ter. From this investigation, we have developed important insights that apply to
decision problems involving multiple hypotheses more generally. However, some
special considerations and issues arise in the more general M-ary hypothesis test-
ing problem. We explore a few of these issues in this section, but emphasize that
our treatment is an especially introductory one.
94 Detection Theory, Decision Theory, and Hypothesis Testing Chap. 2
For this problem, we explore a procedure for choosing one of the M hy-
potheses based on the observed data y so as to minimize the associated expected
cost. This corresponds to designing an M-valued decision rule Ĥ(·) : RK →
{H0 , H1 , . . . , HM −1 }. By analogy to the binary case, we begin by noting that if
y = y and Ĥ(y) = Hm̂ , then the expected cost is
M
X −1
ϕ̃(Hm̂ , y) = Cm̂m Pr [H = Hm | y = y] . (2.106)
m=0
If the Pm are all equal, it follows immediately that (2.112) can be simplified to the
maximum likelihood (ML) rule
Ĥ(y) = Hm̂ where m̂ = arg max py|H (y|Hm ). (2.113)
m∈{0,1,...,M −1}
Example 2.9
Suppose that y is a K-dimensional Gaussian vector under each of the M hypotheses,
with
py|H (y|Hm ) = N (y; mm , I). (2.114)
In this case, applying (2.112)—and recognizing that we may work with the loga-
rithm of the quantity we are maximizing—yields
K 1 T
m̂ = arg max − log(2π) − (y − mm ) (y − mm ) + log Pm . (2.115)
m∈{0,1,...,M −1} 2 2
This can be simplified to
1
m̂ = arg max ℓm (y) − mT mm + log Pm (2.116)
m∈{0,1,...,M −1} 2 m
where
ℓm (y) = hmm , yi , i = 0, 1, . . . , M − 1, (2.117)
with the inner product and associated norm defined by
hx, yi = xT y
96 Detection Theory, Decision Theory, and Hypothesis Testing Chap. 2
and
p p
kyk = hy, yi = yT y.
The Bayesian decision rule in the M-ary case is a natural generalization of the
corresponding binary rule. To see this, we substitute Bayes’ rule (2.111) into the
general rule (2.107) and multiply through by the denominator in (2.111) to obtain
the following rule
M
X −1
Ĥ(y) = Hm̂ where m̂ = arg min Cmj Pj py|H (y|Hj ). (2.119)
m∈{0,1,...,M −1} j=0
py|H (y|Hj )
Lj (y) = , j = 1, 2, . . . , M − 1, (2.120)
py|H (y|H0)
Let us illustrate the resulting structure of the rule for the case of three hy-
potheses. In this case, minimization inherent in the rule (2.121) requires access
to a subset of the results from three comparisons, each of which eliminates one
Sec. 2.8 M -ary Hypothesis Testing 97
L2 (y)
ˆ
H(y)=H1
Ĥ(y) = H2 or H0
(i.e., Ĥ(y)6=H1 )
P1 (C11 − C21 )L1 (y) R P0 (C20 − C10 ) + P2 (C22 − C12 )L2 (y) (2.122b)
Ĥ(y) = H1 or H0
(i.e., Ĥ(y)6=H2 )
Ĥ(y) = H0 or H1
(i.e., Ĥ(y)6=H2 )
P1 (C21 − C01 )L1 (y) R P0 (C00 − C20 ) + P2 (C02 − C22 )L2 (y) (2.122c)
Ĥ(y) = H2 or H1
(i.e., Ĥ(y)6=H0 )
and note that if the last term in (2.123) were not present, this would be exactly
the binary hypothesis test for deciding between H0 and H1 . Assuming that C01 >
C11 (i.e., that it is always more costly to make a mistake than to be correct) and
that C12 < C02 , we see that the last term is negative, biasing the comparison in
favor of H1 . Furthermore, this bias increases as L2 (y) increases, i.e., when the data
indicates H2 to be more and more likely.
Let us next explore some aspects of the performance of the optimum decision rule
for M-ary Bayesian hypothesis testing problems. Generalizing our approach from
the binary case, we begin with an expression for the expected cost obtained by
enumerating all the possible scenarios (i.e., deciding Ĥ(y) = Hi when H = Hj is
correct for all possible values of i and j):
M
X −1 M
X −1 h i
E [C] = Cij Pr Ĥ(y) = Hi H = Hj Pj . (2.124)
i=0 j=0
Sec. 2.8 M -ary Hypothesis Testing 99
so that
pℓ|H (l|Hj ) = N (l; Mmj , MMT ) (2.130)
In this case
h i
Pr Ĥ(y) = Hi H = Hj = Pr [ℓi ≥ ℓm , m = 0, 1, . . . , M − 1 | H = Hj ] (2.131)
which is an M -dimensional integral of the distribution in (2.130) over the set in
which the ith coordinate is at least as large as all of the others.
In order to simply our development, let us restrict our attention to the costs asso-
ciated with the minimum probability-of-error criterion. Specifically, suppose that
the Cij are given by (2.108), and in addition that the hypotheses are equally likely.
In this case (2.124) becomes
M −1
1 X h i
E [C] = Pr(e) = Pr Ĥ(y) 6= Hj H = Hj (2.132)
M j=0
Furthermore, from (2.110) we see that the event
Ej = {Ĥ(y) 6= Hj }
is a union, for k 6= j, of the events
Ekj = {Pr [H = Hk | y] > Pr [H = Hj | y]}. (2.133)
Therefore
h i
Pr Ĥ(y) 6= Hj H = Hj = Pr [Ej | H = Hj ]
" #
[
= Pr Ekj H = Hj
k6=j
X
= Pr [Ekj | H = Hj ]
k6=j
X
− Pr [Ekj ∩ Eij | H = Hj ]
k6=j
i6=j
i6=k
X
+ Pr [Ekj ∩ Eij ∩ Enj | H = Hj ]
k6=j, i6=j, n6=j
k6=i, k6=n, i6=n
−··· (2.134)
Sec. 2.8 M -ary Hypothesis Testing 101
where to obtain the last equality in (2.134) we have used the natural generalization
of the equality (1.5).9
Eq. (2.134) leads us to a natural approximation strategy. Specifically, it fol-
lows that " #
[ X
Pr Ekj H = Hj ≤ Pr [Ekj | H = Hj ] (2.135)
k6=j k6=j
since the sum on the right-hand side adds in more than once the probabilities of
intersection of the Ekj .
The bound (2.135) is referred to as the union bound and is in fact the simplest
(and loosest) of a sequence of possible bounds. Specifically, while (2.135) tells us
that the first term on the right of (2.134) is an upper bound to the desired proba-
bility, the first two terms together are a lower bound10 , the first three terms again
form another (somewhat tighter) upper bound, etc. It is the first of these bounds,
however, that leads to the simplest computations and is of primary interest, yield-
ing
M −1
1 XX
Pr(e) ≤ Pr [Ekj | H = Hj ] . (2.136)
M j=0 k6=j
The union bound is widely used in bit error rate calculations for communications
applications.
We conclude this section by illustrating the application of this bound in the
context of an example.
Example 2.11
Again, we return to Example 2.9 and its continuation 2.10. In this case, via (2.118)
or, equivalently, (2.131)), we have
Ekj = {ℓk > ℓj } = {ℓk − ℓj > 0}
so that
M −1
1 XX
Pr(e) ≤ Pr [ℓk − ℓj > 0 | H = Hj ] (2.137)
M
j=0 k6=j
9
For example, we have
Pr(A ∪ B ∪ C) = Pr(A) + Pr(B) + Pr(C) − Pr(A ∩ B) − Pr(A ∩ C) − Pr(B ∩ C) + Pr(A ∩ B ∩ C).
10
To see this, note that we’ve subtracted out probabilities of intersections of pairs of the Ekj
but in the process have now missed probabilities of intersections of three of the Ekj together.
102 Detection Theory, Decision Theory, and Hypothesis Testing Chap. 2
As a specific numerical example, suppose that there are three hypotheses and
1 0 0
m1 = 0 , m2 = 1 , m3 = 0
0 0 1
In this case, for any k 6= j we have
pℓk −ℓj |H (l|Hj ) = N (l; −1, 2) (2.139)
so that
∞
1 1
Z
2
Pr(e) ≤ 2 √ e−(l+1) /4 dl = 2Q √ (2.140)
0 4π 2
or, equivalently,
Ĥ(y) = Hm̂ if for all m we have (cm̂ − cm )T π(y) ≤ 0 (2.144)
As before, this rule takes the form of a set of comparisons: for each i and k
we have
(ck − ci )T π(y) < 0 =⇒ Ĥ(y) 6= Hi (2.145a)
(ck − ci )T π(y) > 0 =⇒ Ĥ(y) 6= Hk (2.145b)
Sec. 2.8 M -ary Hypothesis Testing 103
π2
(0,1)
π1 +π2 = 1
π1
(1,0)
(a)
π3
(0,0,1)
π1 +π2 +π3 = 1
(0,1,0)
π2
(1,0,0)
π1
(b)
Figure 2.10. Geometry optimum M -
ary Bayesian decision rules in π-space.
104 Detection Theory, Decision Theory, and Hypothesis Testing Chap. 2
(0,0,1)
ˆ
H(y)=H3
ˆ ˆ
H(y)=H
H(y)=H1 2
(1,0,0) (0,1,0)
Figure 2.11. Optimum decision re-
gions with triangle of Fig. 2.10(b).
From the vector space geometry developed in Section 1.7, we see that each of the
equations (ck − ci )T π(y) = 0 defines a subspace that separates the space into two
half-spaces corresponding to (2.145a) and (2.145b). By incorporating all of these
comparisons, the set of subspaces (ck − ci )T π(y) = 0 for all choices of k and i
partition the space into the optimum decision regions. The set of possible proba-
bility vectors, which is a subset of the space, is therefore partitioned into decision
regions. For example, for M = 3, there are three planes through the origin in
Fig. 2.10(b) that partition the space. In Fig. 2.11 we illustrate what this partitioning
looks like restricted to the triangle in Fig. 2.10(b) of possible probability vectors.
Throughout this chapter, recall that we have assumed that the hypotheses Hm were
inherently outcomes of a random variable H. Indeed, we described the probabil-
ity density for the data y under each hypothesis Hm as conditional densities of
the form py|H (y|Hm ). And, in addition, with each Hm we associated an a priori
probability Pm for the hypothesis.
However, it is important to reemphasize that in many problems it may not be
appropriate to view the hypotheses as random—the notion of a priori probabilities
may be rather unnatural. Rather, as we discussed at the outset of the chapter, the
true hypothesis H may be a completely deterministic but unknown quantity. In
these situations, it often makes more sense to view the density for the observed
data not as being conditioned on the unknown hypothesis but rather as being
parameterized by the unknown hypothesis. For such tests it is then appropriate to
Sec. 2.9 Random and Nonrandom Hypotheses, and Selecting Tests 105
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.
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 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
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)
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.
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).
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
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).
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 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.
N
X
2 T
C(a, â) = ka − âk = (a − â) (a − â) = (ai − âi )2 (3.26)
i=1
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 −∞
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
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)
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)
To prove the converse, let us rewrite (3.42) using (3.43) and (3.44) as
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
Λe = E (x̂(y) − x − b) (x̂(y) − x − b)T
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).
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
( x-x )
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
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
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
= λ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)
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)
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)
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
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
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
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)
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.
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
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 + 2 (y − my ) = mx + ρxy (y − my ). (3.124)
σy σ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
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
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
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
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 )
..
.
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
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
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!
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)
−∞
= 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
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
λ 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)
Jy (x)
where the nonnegative quantity Jy (x) is referred to as the Fisher information in y
about x, which is defined by
" 2 #
∂
Jy (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 Jy (x) can’t be calculated.
Second, the notion of referring to (3.168) as an information measure comes
from the fact that Jy (x) is both nonnegative and additive, i.e., whenever
T
y = y1 y2 · · · yM
consists of mutually independent components we have
N
X
Jy (x) = Jyi (x).
i=1
Sec. 3.3 Nonrandom Parameter Estimation 145
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
Jy (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) = Jy (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
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.
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)
Jy (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)
Jy (x)
148 Estimation Theory Chap. 3
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→∞
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
Jy (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
Jy (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/Jy (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
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)
Jy (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 Jy (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
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 σw2
λML (x) = E = (3.219)
h2 h2
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/Jy (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.
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
∂
Jy (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.
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
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
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
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
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 + J−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
∂ 1 ∂
ln py (y; x) = py (y; x), (3.241)
∂x py (y; x) ∂x
158 Estimation Theory Chap. 3
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 J−1
T T −1
y (x)c c Λx̂ (x)c − c J y (x)c ≥ 0. (3.249)
Sec. 3.3 Nonrandom Parameter Estimation 159
However, since J−1 y (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) − J−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) J−1
y (x) f(y). (3.250)
However, since (3.250) holds for any choice of c we must have
e(y) = k(x) J−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) = J−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) J−1 = J−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
But since Jy (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
ϕ(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
Jy (x) = − ϕ(x) = HT Λ−1w 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,
11
This result is referred to as the Gauss-Markov theorem.
162 Estimation Theory Chap. 3
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
Jy (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)
Jy (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
Analog Communication
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.
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
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.
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 Θ′ ,
y [n] = yn (t∗ )
d
= A cos Ω0 cos φ n + Θ + wn (t∗ )
c
= A cos(ω0 n + Θ) + w [n], (3.285)
where Θ = Θ′ + Ω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 φ.
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
[Jy (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
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
Jy (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
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).
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
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.
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
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
σ 2 = 10 log10 (E |WN (ejω )|2 ) corre-
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
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
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
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
10
Periodogram (dB)
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
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.
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)
Jy (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 Jy (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
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)
where
ϕ(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
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
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
d2 ∂2
mr (a) =E r(y; a)
da2 a=x ∂a2
a=x
2
∂
=E ln py (y; a)
∂a2
a=x
= −Jy (x). (3.335)
d2
mr (a) = −Jy (x). (3.336)
da2 a=x
180 Estimation Theory Chap. 3
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
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.
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.
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)
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)
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
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
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)
As we’ll see, this function and its first and second derivatives, respectively
N −1
′ 2 X j(2ω0 n+2Θ)
ξ (ω0 ) = ne (3.352)
N n=0
and
N −1
′′ 4 X 2 j(2ω0 n+2Θ)
ξ (ω0 ) = n e , (3.353)
N n=0
∂2
[Jy (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 )} , (3.354)
2σ 2 2
Sec. 3.B Fisher Information Calculations for Sinusoid Estimation 185
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 ) ∼ o(N), ξ ′′ (ω0 ) ∼ o(N 2 ). (3.364)
where
ϕ1 (ω0 ) = ϕ(α̂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)
where −1
ϕ2 (ω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
corresponding to (3.302).