4.1.1 MMSE Estimation

The solution based on the quadratic cost function (4.6) is called the minimum mean square error estimator, also called the minimum variance estimator for reasons that will become clear in a moment. Substitution of Equations (4.6) and (4.9) in (4.10) gives

(4.12)

Differentiating the function between braces with respect to (see Appendix B.4) and equating this result to zero yields a system of M linear equations, the solution of which is

(4.13)

The conditional risk of this solution is the sum of the variances of the estimated parameters:

(4.14)

and hence the name ‘minimum variance estimator’.

4.1.2 MAP Estimation

If the uniform cost function is chosen, the conditional risk (4.9) becomes

(4.15)

The estimate that now minimizes the risk is called the maximum a posteriori (MAP) estimate:

(4.16)

This solution equals the mode (maximum) of the posterior probability. It can also be written entirely in terms of the prior probability densities and the conditional probabilities:

(4.17)

This expression is similar to the one of MAP classification (see Equation (3.12)).

4.1.3 The Gaussian Case with Linear Sensors

Suppose that the parameter vector x has a normal distribution with expectation vector E[x] = and covariance matrix C_x. In addition, suppose that the measurement vector can be expressed as a linear combination of the parameter vector corrupted by additive Gaussian noise:

(4.18)

where v is an N-dimensional random vector with zero expectation and covariance matrix C_v, x and v are uncorrelated, and H is an N × M matrix.

The assumption that both x and v are normal implies that the conditional probability density of z is also normal. The conditional expectation vector of z equals

(4.19)

The conditional covariance matrix of z is C_z|x = C_v.

Under these assumptions the posterior distribution of x is normal as well. Application of Equation (4.2) yields the following expressions for the MMSE estimate and the corresponding covariance matrix:

(4.20)

The proof is left as an exercise for the reader (see Exercise 4). Note that , being the posterior expectation, is the MMSE estimate .

The posterior expectation, E[x|z], consists of two terms. The first term is linear in z. It represents the information coming from the measurement. The second term is linear in , representing the prior knowledge. To show that this interpretation is correct it is instructive to see what happens at the extreme ends: either no information from the measurement or no prior knowledge:

• The measurements are useless if the matrix H is virtually zero or if the noise is too large, that is C^{− 1}_v is too small. In both cases, the second term in Equation (4.20) dominates. In the limit, the estimate becomes with covariance matrix C_x, that is the estimate is purely based on prior knowledge.

• On the other hand, if the prior knowledge is weak, that is if the variances of the parameters are very large, the inverse covariance matrix C^{− 1}_x tends to zero. In the limit, the estimate becomes

  (4.21)

  In this solution, the prior knowledge, that is , is completely ruled out.

Note that the mode of a normal distribution coincides with the expectation. Therefore, in the linear-Gaussian case, MAP estimation and MMSE estimation coincide: .

4.1.4 Maximum Likelihood Estimation

In many practical situations the prior knowledge needed in MAP estimation is not available. In these cases, an estimator that does not depend on prior knowledge is desirable. One attempt in that direction is the method referred to as maximum likelihood estimation (ML estimation). The method is based on the observation that in MAP estimation, Equation (4.17), the peak of the first factor p(z|x) is often in an area of x in which the second factor p(x) is almost constant. This holds true especially if little prior knowledge is available. In these cases, the prior density p(x) does not affect the position of the maximum very much. Discarding the factor and maximizing the function p(z|x) solely gives the ML estimate:

(4.22)

Regarded as a function of x the conditional probability density is called the likelihood function; hence the name ‘maximum likelihood estimation’.

Another motivation for the ML estimator is when we change our viewpoint with respect to the nature of the parameter vector x. In the Bayesian approach x is a random vector, statistically defined by means of probability densities. In contrast, we may also regard x as a non-random vector whose value is simply unknown. This is the so-called Fisher approach. In this view, there are no probability densities associated with x. The only density in which x appears is p(z|x), but here x must be regarded as a parameter of the density of z. From all estimators discussed so far, the only estimator that can handle this deterministic point of view on x is the ML estimator.

Example 4.4 Maximum likelihood estimation of the backscattering coefficient The maximum likelihood estimator for the backscattering coefficient (see the previous examples) is found by maximizing Equation (4.5):

(4.23)

The estimator is depicted in Figure 4.6 together with the MAP estimator. The figure confirms the statement above that in areas of flat prior probability density the MAP estimator and the ML estimator coincide. However, the figure also reveals that the ML estimator can produce an estimate of the backscattering coefficient that is larger than one – a physical impossibility. This is the price that we have to pay for not using prior information about the physical process.

If the measurement vector z is linear in x and corrupted by additive Gaussian noise, as given in Equation (4.18), the likelihood of x is given in Equation (4.21). Thus, in that case:

(4.24)

A further simplification is obtained if we assume that the noise is white, that is C_v ∼ I:

(4.25)

The operation (H^TH)^–1H^T is the pseudo inverse of H. Of course, its validity depends on the existence of the inverse of H^TH. Usually, such is the case if the number of measurements exceeds the number of parameters, that is N > M, which is if the system is overdetermined.

4.1.5 Unbiased Linear MMSE Estimation

The estimators discussed in the previous sections exploit full statistical knowledge of the problem. Designing such an estimator is often difficult. The first problem that arises is adequate modelling of the probability densities. Such modelling requires detailed knowledge of the physical process and the sensors. Once we have the probability densities, the second problem is how to minimize the conditional risk. Analytic solutions are often hard to obtain. Numerical solutions are often burdensome.

If we constrain the expression for the estimator to some mathematical form, the problem of designing the estimator boils down to finding the suitable parameters of that form. An example is the unbiased linear MMSE estimator, with the following form:²

(4.26)

The matrix K and the vector a must be optimized during the design phase in order to match the behaviour of the estimator to the problem at hand. The estimator has the same optimization criterion as the MMSE estimator, that is a quadratic cost function. The constraint results in an estimator that is not as good as the (unconstrained) MMSE estimator, but it requires only knowledge of moments up to the order two, that is expectation vectors and covariance matrices.

The starting point is the overall risk expressed in Equation (4.11). Together with the quadratic cost function we have

(4.27)

The optimal unbiased linear MMSE estimator is found by minimizing R with respect to K and a. Hence we differentiate R with respect to a and equate the result to zero:

yielding

(4.28)

with and the expectations of x and z.

Substitution of a back into Equation (4.27), differentiation with respect to K and equating the result to zero (see also Appendix B.4):

yields

(4.29)

where C_z is the covariance matrix of z and the cross-covariance matrix between x and z.

Example 4.5 Unbiased linear MMSE estimation of the backscattering coefficient In the scalar case, the linear MMSE estimator takes the form:

(4.30)

where cov[x, z] is the covariance of x and z. In the backscattering problem, the required moments are difficult to obtain analytically. However, they are easily estimated from the population of the 500 realizations shown in Figure 4.6 using techniques from Chapter 5. The resulting estimator is shown in Figure 4.6. Matlab^® code to plot the ML and unbiased linear MMSE estimates of the backscattering coefficient on a data set is given in Listing 4.2.

Listing 4.2 Matlab^® code for an unbiased linear MMSE estimation.

 prload scatter;           % Load dataset (zset,xset)
 z = 0.005:0.005:1.5;      % Interesting range of z
 x_ml = z;                 % Maximum likelihood
 mu_x = mean(xset); mu_z = mean(zset);
 K = ((xset-mu_x)'*(zset-mu_z))*inv((zset-mu_z)'*(zset-mu_z));
 a = mu_x K*mu_z;
 x_ulmse = K*z + a;        % Unbiased linear MMSE
 figure; clf; plot(zset,xset,'.'); hold on;
 plot(z,x_ml,'k-'); plot(z,x_ulmse,'k- -');

Linear Sensors

The linear MMSE estimator takes a particular form if the sensory system is linear and the sensor noise is additive:

(4.31)

This case is of special interest because of its crucial role in the Kalman filter (to be discussed in Chapter 5). Suppose that the noise has zero mean with covariance matrix C_v. In addition, suppose that x and v are uncorrelated, that is C_xv = 0. Under these assumptions the moments of z are as follows:

(4.32)

The proof is left as an exercise for the reader.

Substitution of Equations (4.32), (4.28) and (4.29) in (4.26) gives rise to the following estimator:

(4.33)

This version of the unbiased linear MMSE estimator is the so-called Kalman form of the estimator.

Examination of Equation (4.20) reveals that the MMSE estimator in the Gaussian case with linear sensors is also expressed as a linear combination of and z. Thus, in this special case (that is, Gaussian densities + linear sensors) is a linear estimator. Since and are based on the same optimization criterion, the two solutions must be identical here: . We conclude that Equation (4.20) is an alternative form of (4.33). The forms are mathematically equivalent (see Exercise 5).

The interpretation of is as follows. The term represents the prior knowledge. The term is the prior knowledge that we have about the measurements. Therefore, the factor is the informative part of the measurements (called the innovation). The so-called Kalman gain matrix K transforms the innovation into a correction term that represents the knowledge that we have gained from the measurements.

4.2.1 Bias and Covariance

The error e is composed of two parts. One part is the one that does not change value if we repeat the experiment over and over again. It is the expectation of the error, called the bias. The other part is the random part and is due to sensor noise and other random phenomena in the sensory system. Hence, we have

(4.35)

If x is a scalar, the variance of an estimator is the variance of e. Thus the variance quantifies the magnitude of the random part. If x is a vector, each element of e has its own variance. These variances are captured in the covariance matrix of e, which provides an economical and also a more complete way to represent the magnitude of the random error.

The application of the expectation and variance operators to e needs some discussion. Two cases must be distinguished. If x is regarded as a non-random, unknown parameter, then x is not associated with any probability density. The only randomness that enters the equations is due to the measurements z with density p(z|x). However, if x is regarded as random, it does have a probability density. We then have two sources of randomness, x and z.

We start with the first case, which applies to, for instance, the maximum likelihood estimator. Here, the bias b(x) is given by

(4.35)

The integral extends over the full space of z. In general, the bias depends on x. The bias of an estimator can be small or even zero in one area of x, whereas in another area the bias of that same estimator might be large.

In the second case, both x and z are random. Therefore, we define an overall bias b by taking the expectation operator now with respect to both x and z:

(4.36)

The integrals extend over the full space of x and z.

The overall bias must be considered as an average taken over the full range of x. To see this, rewrite p(x, z) = p(z|x)p(x) to yield

(4.37)

where b(x) is given in Equation (4.35).

If the overall bias of an estimator is zero, then the estimator is said to be unbiased. Suppose that in two different areas of x the biases of an estimator have opposite sign; then these two opposite biases may cancel out. We conclude that, even if an estimator is unbiased (i.e. its overall bias is zero), then this does not imply that the bias for a specific value of x is zero. Estimators that are unbiased for every x are called absolutely unbiased.

The variance of the error, which serves to quantify the random fluctuations, follows the same line of reasoning as that of the bias. First, we determine the covariance matrix of the error with non-random x:

(4.38)

As before, the integral extends over the full space of z. The variances of the elements of e are at the diagonal of C_e(x).

The magnitude of the full error (bias + random part) is quantified by the so-called mean square error (the second-order moment matrix of the error):

(4.39)

It is straightforward to prove that

(4.40)

This expression underlines the fact that the error is composed of a bias and a random part.

The overall mean square error M_e is found by averaging M_e(x) over all possible values of x:

(4.41)

Finally, the overall covariance matrix of the estimation error is found as

(4.42)

The diagonal elements of this matrix are the overall variances of the estimation errors.

The MMSE estimator and the unbiased linear MMSE estimator are always unbiased. To see this, rewrite Equation (4.36) as follows:

(4.43)

The inner integral is identical to zero and thus b must be zero. The proof of the unbiasedness of the unbiased linear MMSE estimator is left as an exercise.

Other properties related to the quality of an estimator are stability and robustness. In this context, stability refers to the property of being insensitive to small random errors in the measurement vector. Robustness is the property of being insensitive to large errors in a few elements of the measurements (outliers) (see Section 4.3.2). Often, the enlargement of prior knowledge increases both the stability and the robustness.

Example 4.6 Bias and variance in the backscattering application Figure 4.7 shows the bias and variance of the various estimators discussed in the previous examples. To enable a fair comparison between bias and variance in comparable units, the square root of the latter, that is the standard deviation, has been plotted. Numerical evaluation of Equations (4.37), (4.41) and (4.42) yields:³

From this, and from Figure 4.7, we observe that:

• The overall bias of the ML estimator appears to be zero. Therefore, in this example, the ML estimator is unbiased (together with the two MMSE estimators, which are intrinsically unbiased). The MAP estimator is biased.
• Figure 4.7 shows that for some ranges of x the bias of the MMSE estimator is larger than its standard deviation. Nevertheless, the MMSE estimator outperforms all other estimators with respect to overall bias and variance. Hence, although a small bias is a desirable property, sometimes the overall performance of an estimator can be improved by allowing a larger bias.
• The ML estimator appears to be linear here. As such, it is comparable with the unbiased linear MMSE estimator. Of these two linear estimators, the unbiased linear MMSE estimator outperforms the ML estimator. The reason is that – unlike the ML estimator – the ulMMSE estimator exploits prior knowledge about the parameter. In addition, the ulMMSE estimator is more applicable to the evaluation criterion.
• Of the two non-linear estimators, the MMSE estimator outperforms the MAP estimator. The obvious reason is that the cost function of the MMSE estimator matches the evaluation criterion.
• Of the two MMSE estimators, the non-linear MMSE estimator outperforms the linear one. Both estimators have the same optimization criterion, but the constraint of the ulMMSE degrades its performance.

4.2.2 The Error Covariance of the Unbiased Linear MMSE Estimator

We now return to the case of having linear sensors, z = Hx + v, as discussed in Section 4.1.5. The unbiased linear MMSE estimator appeared to be (see Equation (4.33))

where C_v and C_x are the covariance matrices of v and x; is the (prior) expectation vector of x. As stated before, the is unbiased.

Due to the unbiasedness of , the mean of the estimation error – x is zero. The error covariance matrix, C_e, of e expresses the uncertainty that remains after having processed the measurements. Therefore, C_e is identical to the covariance matrix associated with the posterior probability density. It is given by (Equation (4.20))

(4.44)

The inverse of a covariance matrix is called an information matrix. For instance, C^{− 1}_e is a measure of the information provided by the estimate . If the norm of C^{− 1}_e is large, then the norm of C_e must be small, implying that the uncertainty in is small as well. Equation (4.44) shows that C^{− 1}_e is made up of two terms. The term C^{− 1}_x represents the prior information provided by . The matrix C^{− 1}_v represents the information that is given by z about the vector H_x. Therefore, the matrix H^TC^{− 1}_vH represents the information about x provided by z. The two sources of information add up, so the information about x provided by is C^{− 1}_e = C_x^{− 1} + H^TC^{− 1}_vH.

Using the matrix inversion lemma (4.10) the expression for the error covariance matrix can be given in an alternative form:

(4.45)

The matrix S is called the innovation matrix because it is the covariance matrix of the innovations . The factor is a correction term for the prior expectation vector . Equation (4.45) shows that the prior covariance matrix C_x is reduced by the covariance matrix KSK^T of the correction term.

4.3.1 Least Squares Fitting

The most common error norm is the squared Euclidean norm, also called the sum of squared differences (SSD), or simply the LS norm (least squared error norm):

(4.48)

The least squares fit, or least squares estimate (LS), is the parameter vector that minimizes this norm:

(4.49)

If v is random with a normal distribution, zero mean and covariance matrix C_v = σ²_vI, the LS estimate is identical to the ML estimate: . To see this, it suffices to note that in the Gaussian case the likelihood takes the form

(4.50)

The minimization of Equation (4.48) is equivalent to the minimization of Equation (4.50). If the measurement function is linear, that is z = Hx + v, and H is an N × M matrix having a rank M with M < N, then according to Equation (4.25)

(4.51)

Example 4.7 Repeated measurements Suppose that a scalar parameter x is N times repeatedly measured using a calibrated measurement device: z_n = x + v_n. These repeated measurements can be represented by a vector z = [z₁, …, z_N]^T. The corresponding measurement matrix is H = [1, …, 1]^T. Since (H^TH)^{− 1} = 1/N, the resulting least squares fit is

In other words, the best fit is found by averaging the measurements.

Non-linear sensors

If h(.) is non-linear, an analytic solution of Equation (4.49) is often difficult. One is compelled to use a numerical solution. For that, several algorithms exist, such as ‘Gauss–Newton’, ‘Newton–Raphson’, ‘steepest descent’ and many others. Many of these algorithms are implemented within Matlab^®'s optimization toolbox. The ‘Gauss–Newton’ method will be explained shortly.

Assuming that some initial estimate x_ref is available we expand Equation (4.46) in a Taylor series and neglect all terms of order higher than two:

(4.52)

where H_ref is the Jacobian matrix of h(.) evaluated at x_ref (see Appendix B.4). With such a linearization, Equation (4.51) applies. Therefore, the following approximate value of the LS estimate is obtained:

(4.53)

A refinement of the estimate could be achieved by repeating the procedure with the approximate value as reference. This suggests an iterative approach. Starting with some initial guess , the procedure becomes as follows:

(4.54)

In each iteration, the variable i is incremented. The iterative process stops if the difference between x(i + 1) and x(i) is smaller than some predefined threshold. The success of the method depends on whether the first initial guess is already close enough to the global minimum. If not, the process will either diverge or get stuck in a local minimum.

Example 4.8 Estimation of the diameter of a blood vessel In vascular X-ray imaging, one of the interesting parameters is the diameter of blood vessels. This parameter provides information about a possible constriction. As such, it is an important aspect in cardiologic diagnosis.

Figure 4.8(a) is a (simulated) X-ray image of a blood vessel of the coronary circulation. The image quality depends on many factors. Most important are the low-pass filtered noise (called quantum mottle) and the image blurring due to the image intensifier.

Figure 4.8(b) shows the one-dimensional, vertical cross-section of the image at a location, as indicated by the two black arrows in Figure 4.8(a). Suppose that our task is to estimate the diameter of the imaged blood vessel from the given cross-section. Hence, we define a measurement vector z whose elements consist of the pixel grey values along the cross-section.

The parameter of interest is the diameter D. However, other parameters might be unknown as well, for example the position and orientation of the blood vessel, the attenuation coefficient or the intensity of the X-ray source. This example will be confined to the case where the only unknown parameters are the diameter D and the position y₀ of the image blood vessel in the cross-section. Thus, the parameter vector is two-dimensional: x^T = [D y₀]. Since all other parameters are assumed to be known, a radiometric model can be worked out to a measurement function h(x), which quantitatively predicts the cross-section and thus also the measurement vector z for any value of the parameter vector x.

With this measurement function it is straightforward to calculate the LS norm ‖ϵ‖²₂ for a couple of values of x. Figure 4.8(c) is a graphical representation of that. It appears that the minimum of ‖ϵ‖²₂ is obtained if . The true diameter is D = 0.40 mm. The obtained fitted cross-section is also shown in Figure 4.8(b).

Note that the LS norm in Figure 4.8(c) is a smooth function of x. Hence, the convergence region of a numerical optimizer will be large.

4.3.2 Fitting Using a Robust Error Norm

Suppose that the measurement vector in an LS estimator has a few elements with large measurement errors, the so-called outliers. The influence of an outlier is much larger than one of the others because the LS estimator weights the errors quadratically. Consequently, the robustness of LS estimation is poor.

Much can be improved if the influence is bounded in one way or another. This is exactly the general idea of applying a robust error norm. Instead of using the sum of squared differences, the objective function of (4.48) becomes

(4.55)

where ρ(.) measures the size of each individual residual . This measure should be selected such that above a given level of ϵ_n its influence is ruled out. In addition, one would like to have ρ(.) being smooth so that numerical optimization of ‖ϵ‖_robust is not too difficult. A suitable choice (among others) is the so-called Geman–McClure error norm:

(4.56)

A graphical representation of this function and its derivative is shown in Figure 4.9. The parameter σ has a soft threshold value. For values of smaller than about σ, the function follows the LS norm. For values larger than σ, the function becomes saturated. Consequently, for small values of the derivative of ρ(.) is nearly a constant. However, for large values of , that is for outliers, it becomes nearly zero. Therefore, in a Gauss–Newton style of optimization, the Jacobian matrix is virtually zero for outliers. Only residuals that are about as large as σ or smaller than that play a role.

Example 4.9 Robust estimation of the diameter of a blood vessel If in Example 4.8 the diameter is estimated near the bifurcation (as indicated in Figure 4.8(a) by the white arrows) a large modelling error occurs because of the branching vessel (see the cross-section in Figure 4.10(a)). These modelling errors are large compared to the noise and they should be considered as outliers. However, Figure 4.10(b) shows that the error landscape ‖ϵ‖²₂ has its minimum at . The true value is D = 0.40 mm. Furthermore, the minimum is less pronounced than the one in Figure 4.8(c), and therefore also less stable.

Note also that in Figure 4.10(a) the position found by the LS estimator is in the middle between the two true positions of the two vessels.

Figure 4.11 shows the improvements that were obtained by applying a robust error norm. The threshold σ is selected just above the noise level. For this setting, the error landscape clearly shows two pronounced minima corresponding to the two blood vessels. The global minimum is reached at . The estimated position now clearly corresponds to one of the two blood vessels, as shown in Figure 4.11(a).

4.3.3 Regression

Regression is the act of deriving an empirical function from a set of experimental data. Regression analysis considers the situation involving pairs of measurements (t, z). The variable t is regarded as a measurement without any appreciable error; t is called the independent variable. We assume that some empirical function f(.) is chosen that (hopefully) can predict z from the independent variable t. Furthermore, a parameter vector x can be used to control the behaviour of f(.). Hence, the model is

(4.57)

where f(.,.) is the regression curve and ϵ represents the residual, that is the part of z that cannot be predicted by f (.,.). Such a residual can originate from sensor noise (or other sources of randomness), which makes the prediction uncertain, but it can also be caused by an inadequate choice of the regression curve.

The goal of regression is to determine an estimate of the parameter vector x based on N observations (t_n, z_n), n = 0,…, N – 1 such that the residuals ϵ_n are as small as possible. We can stack the observations z_n in a vector z. Using Equation (4.57), the problem of finding can be transformed to the standard form of Equation (4.47):

(4.58)

where is the vector that embodies the residuals ϵ_n.

Since the model is in the standard form, x can be estimated with a least squares approach as in Section 4.3.1. Alternatively, we use a robust error norm as defined in Section 4.3.2. The minimization of such a norm is called robust regression analysis.

In the simplest case, the regression curve f(t, x) is linear in x. With that, the model takes the form , and thus the solution of Equation (4.51) applies. As an example, we consider polynomial regression for which the regression curve is a polynomial of order M – 1:

(4.59)

If, for instance, M = 3, then the regression curve is a parabola described by three parameters. These parameters can be found by least squares estimation using the following model:

(4.60)

Example 4.10 The calibration curve of a level sensor In this example, the goal is to determine a calibration curve of a level sensor to be used in a water tank. For that purpose, a second measurement system is available with a much higher precision than the ‘sensor under test’. The measurement results of the second system serve as a reference. Figure 4.12 shows the observed errors of the sensor versus the reference values. Here, 46 pairs of measurements are shown. A zero-order (fit with a constant), a first-order (linear fit) and a tenth-order polynomial are fitted to the data. As can be seen, the constant fit appears to be inadequate for describing the data (the model is too simple). The first-order polynomial describes the data reasonably well and is also suitable for extrapolation. The tenth-order polynomial follows the measurement points better and also the noise. It cannot be used for extrapolation because it is an example of overfitting the data. This occurs whenever the model has too many degrees of freedom compared to the number of data samples.

Listing 4.3 illustrates how to fit and evaluate polynomials using Matlab^®’s polyfit() and polyval() routines.

Listing 4.3 Matlab^® code for polynomial regression.

 prload levelsensor;                      % Load dataset (t,z)
 figure; clf; plot(t,z,`k.'); hold on;    % Plot it
 y = 0:0.2:30; M = [1210];plotstring = {`k- -',`k-',`k:'};
 for m = 1:3
   p = polyfit(t,z,M(m)-1);               % Fit polynomial
   z_hat = polyval(p,y);                  % Calculate plot points
   plot(y,z_hat,plotstring{m});           % and plot them
 end;
 axis([0 30 0.3 0.2]);

4.5 Selected Bibliography

1. Kay, S.M., Fundamentals of Statistical Signal Processing – Estimation Theory, Prentice Hall, New Jersey, 1994.
2. Papoulis, A., Probability, Random Variables and Stochastic Processes, Third Edition, McGraw-Hill, New York, 1991.
3. Sorenson, H.W., Parameter Estimation: Principles and Problems, Marcel Dekker, New York, 1980.
4. Tarantola, A., Inverse Problem Theory: Methods for Data Fitting and Model Parameter Estimation, Elsevier, New York, 1987.
5. Van Trees, H.L., Detection, Estimation and Modulation Theory, Vol. 1, John Wiley & Sons, Inc., New York, 1968.