The Bayesian Philosophy

Reference:
Kay S M. Fundamentals of statistical signal processing[M]. Prentice Hall PTR, 1993. (Ch. 10 - 10.6)
Slides of ET4386, TUD

We now depart from the classical approach to statistical estimation in which the parameter $\theta$ of interest is assumed to be a deterministic but unknown constant. Instead, we assume that $\theta$ is a random variable whose particular realization we must estimate.

Motivation:

• The Bayesian approach can incorporate the available prior knowledge into our estimator.
• Bayesian estimation is useful in situations where an MVU estimator cannot be found, as for example, when the variance of an unbiased estimator may not be uniformly less than that of all other unbiased estimators. In this instance, it may be true that for most values of the parameter an estimator can be found whose mean square error may be less than that of all other estimators.

Prior Knowledge and Estimation

It is a fundamental rule of estimation theory that the use of prior knowledge will lead to a more accurate estimator.

Reconsider the DC in WGN problem: $x[n]=A+w[n],n=0,\cdots N-1,w[n]\sim \mathcal{N}(0,{\sigma }^2)$. It can be shown that the MVU estimator of $A$ is the sample mean $\bar{x}$. However, this assumed that $A$ could take on any value in the interval $-\infty <A<\infty$. Due to physical constraints it may be more reasonable to assume that $A$ can take on only values in the finite interval $-A_0\le A\le A_0$. Obviously, we would expect to improve our estimation if we used the truncated sample mean estimator
$$\check{A}=\{\begin{array}{cccc}& -A_0& & \bar{x}<-A_0\\ & \bar{x}& & -A_0\le \bar{x}\le A_0\\ & A_0& & \bar{x}>A_0\end{array}$$
which would be consistent with the known constraints. Such an estimator would have the PDF
p A ˇ ( ξ ; A ) = Pr ⁡ { x ˉ ≤ − A 0 } δ ( ξ + A 0 ) + p A ^ ( ξ ; A ) [ u ( ξ + A 0 ) − u ( ξ − A 0 ) ] + Pr ⁡ { x ˉ ≥ A 0 } δ ( ξ − A 0 ) p_{\check A}(\xi;A)=\Pr\{\bar x \le -A_0 \}\delta(\xi+A_0)+p_{\hat A}(\xi;A)[u(\xi+A_0)-u(\xi-A_0)]+\Pr\{\bar x \ge A_0 \}\delta(\xi-A_0) pAˇ​(ξ;A)=Pr{xˉ≤−A0​}δ(ξ+A0​)+pA^​(ξ;A)[u(ξ+A0​)−u(ξ−A0​)]+Pr{xˉ≥A0​}δ(ξ−A0​)
where $u(x)$ is the unit step function.

It is seen that $\check{A}$ is a biased estimator. However, if we compare the MSE of the two estimators, we note that for any $A$ in the interval $-A_0\le A\le A_0$
mse ⁡ ( A ^ ) = ∫ − ∞ ∞ ( ξ − A ) 2 p A ^ ( ξ ; A ) d ξ = ∫ − ∞ − A 0 ( ξ − A ) 2 p A ^ ( ξ ; A ) d ξ + ∫ − A 0 A 0 ( ξ − A ) 2 p A ^ ( ξ ; A ) d ξ + ∫ A 0 ∞ ( ξ − A ) 2 p A ^ ( ξ ; A ) d ξ > ∫ − ∞ − A 0 ( − A 0 − A ) 2 p A ^ ( ξ ; A ) d ξ + ∫ − A 0 A 0 ( ξ − A ) 2 p A ^ ( ξ ; A ) d ξ + ∫ A 0 ∞ ( A 0 − A ) 2 p A ^ ( ξ ; A ) d ξ = ( − A 0 − A ) 2 Pr ⁡ { x ˉ ≤ − A 0 } + ∫ − A 0 A 0 ( ξ − A ) 2 p A ^ ( ξ ; A ) d ξ + ( A 0 − A ) 2 Pr ⁡ { x ˉ ≥ A 0 } = mse ⁡ ( A ~ ) \begin{aligned} \operatorname{mse}(\hat{A})=& \int_{-\infty}^{\infty}(\xi-A)^{2} p_{\hat{A}}(\xi ; A) d \xi \\ =& \int_{-\infty}^{-A_{0}}(\xi-A)^{2} p_{\hat{A}}(\xi ; A) d \xi+\int_{-A_{0}}^{A_{0}}(\xi-A)^{2} p_{\hat{A}}(\xi ; A) d \xi +\int_{A_{0}}^{\infty}(\xi-A)^{2} p_{\hat{A}}(\xi ; A) d \xi \\>&\int_{-\infty}^{-A_{0}}(-A_0-A)^{2} p_{\hat{A}}(\xi ; A) d \xi+\int_{-A_{0}}^{A_{0}}(\xi-A)^{2} p_{\hat{A}}(\xi ; A) d \xi +\int_{A_{0}}^{\infty}(A_0-A)^{2} p_{\hat{A}}(\xi ; A) d \xi \\ =& \left(-A_{0}-A\right)^{2} \Pr\{\bar x \le -A_0 \}+\int_{-A_{0}}^{A_{0}}(\xi-A)^{2} p_{\hat{A}}(\xi ; A) d \xi +\left(A_{0}-A\right)^{2} \Pr\{\bar x \ge A_0 \} \\ =& \operatorname{mse}(\tilde{A}) \end{aligned} mse(A^)==>==​∫−∞∞​(ξ−A)2pA^​(ξ;A)dξ∫−∞−A0​​(ξ−A)2pA^​(ξ;A)dξ+∫−A0​A0​​(ξ−A)2pA^​(ξ;A)dξ+∫A0​∞​(ξ−A)2pA^​(ξ;A)dξ∫−∞−A0​​(−A0​−A)2pA^​(ξ;A)dξ+∫−A0​A0​​(ξ−A)2pA^​(ξ;A)dξ+∫A0​∞​(A0​−A)2pA^​(ξ;A)dξ(−A0​−A)2Pr{xˉ≤−A0​}+∫−A0​A0​​(ξ−A)2pA^​(ξ;A)dξ+(A0​−A)2Pr{xˉ≥A0​}mse(A~)​
Hence, $\check{A}$, the truncated sample mean estimator, is better than the sample mean estimator in terms of MSE. Although $\hat{A}$ is still the MVU estimator, we have been able to reduce the mean square error by allowing the estimator to be biased.

---

In as much as we have been able to produce a better estimator, the question arises as to whether an optimal estimator exists for this problem. (In the classical case the MSE criterion of optimality usually led to unrealizable estimators. We shall see that this is not a problem in the Bayesian approach.)

With knowledge only of the interval and no inclination as to whether $A$ should be nearer any particular value, it makes sense to assign a $\mathcal{U}[-A_0,A_0]$ PDF to the random variable $A$. The overall data model is different from the classical approach

Now we can incorporate our knowledge of how $A$ was chosen. Define the Bayesian MSE as
$$\mathrm{B}\mathrm{m}\mathrm{s}\mathrm{e}(\hat{A})=E[(A-\hat{A}{)}^2]\text{\hspace{1em}(BP.1)}$$
To appreciate the difference compare the classical MSE
$$\mathrm{m}\mathrm{s}\mathrm{e}(\hat{A})=\int (\hat{A}-A{)}^{2}p(\mathbf{x};A)d\mathbf{x}\text{\hspace{1em}(BP.2)}$$
to the Bayesian MSE
$$\mathrm{B}\mathrm{m}\mathrm{s}\mathrm{e}(\hat{A})=\int (A-\hat{A}{)}^2\left[\int p(\mathbf{x},A)d\mathbf{x}\right]dA=\int \int (A-\hat{A}{)}^{2}p(\mathbf{x},A)d\mathbf{x}dA\text{\hspace{1em}(BP.3)}$$
Note that whereas the classical MSE will depend on $A$, and hence estimators that attempt to minimize the MSE will usually depend on $A$, the Bayesian MSE will not. In effect, we have integrated the parameter dependence away.

Now we can derive the estimator that minimizes the Bayesian MSE. First, we use Bayes’ theorem to write
$$p(\mathbf{x},A)=p(A|\mathbf{x})p(\mathbf{x})$$
so that
$$\mathrm{Bmse}(\hat{A})=\int \left[\int (A-\hat{A}{)}^{2}p(A|\mathbf{x})dA\right]p(\mathbf{x})d\mathbf{x}\text{\hspace{1em}(BP.4)}$$
Now since $p(\mathbf{x})\ge 0$ for all $\mathbf{x},$ if the integral in brackets can be minimized for each $\mathbf{x},$ then the Bayesian MSE will be minimized. Hence, fixing $\mathbf{x}$ so that $\hat{A}$ is a scalar variable, we have
$$\begin{array}{rl}\frac{\partial }{\partial \hat{A}}\int (A-\hat{A}{)}^{2}p(A|\mathbf{x})dA& =\int \frac{\partial }{\partial \hat{A}}(A-\hat{A}{)}^{2}p(A|\mathbf{x})dA\\ & =\int -2(A-\hat{A})p(A|\mathbf{x})dA\\ & =-2\int Ap(A|\mathbf{x})dA+2\hat{A}\int p(A|\mathbf{x})dA\end{array}$$
which when set equal to zero results in
$$\hat{A}=\int Ap(A|\mathbf{x})dA=E(A|\mathbf{x})\text{\hspace{1em}(BP.5)}$$
since the conditional PDF must integrate to $1$. It is seen that the optimal estimator in terms of minimizing the Bayesian MSE is the mean of the posterior PDF $p(A|\mathbf{x})$. We will henceforth term the estimator that minimizes the Bayesian MSE the minimum mean square error (MMSE) estimator.

---

To determine the MMSE, we first require the posterior PDF. We can use Bayes’ rule to determine it as
$$p(A|\mathbf{x})=\frac{p(\mathbf{x}|A)p(A)}{p(\mathbf{x})}=\frac{p(\mathbf{x}|A)p(A)}{\int p(\mathbf{x}|A)p(A)dA}\text{\hspace{1em}(BP.6)}$$
Note that the denominator is just a normalizing factor, independent of $A$, needed to ensure that $p(A|\mathbf{x})$ integrates to $1$. Then the MMSE
$$\hat{A}=E(A|\mathbf{x})=\frac{\int Ap(\mathbf{x}|A)p(A)dA}{\int p(\mathbf{x}|A)p(A)dA}\text{\hspace{1em}(BP.7)}$$
If we continue our example, we recall that the prior PDF $p(A)$ is $\mathcal{U}[-A_0,A_0]$. To specify the conditional PDF $p(\mathbf{x}|A)$ we need to further assume that the choice of $A$ via $p(A)$ does not affect the PDF of the noise samples or that $w[n]$ is independent of $A$. Then for $n=0,1,\cdots ,N-1$
$$\begin{array}{rl}{p}_x(x[n]|A)& =p_w(x[n]-A|A)\\ & =p_w(x[n]-A)\\ & =\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \left[-\frac{1}{2{\sigma }^2}(x[n]-A{)}^2\right]\end{array}$$
and therefore
$$p(\mathbf{x}|A)=\frac{1}{(2\pi {\sigma }^2{)}^{N/2}}\exp \left[-\frac{1}{2{\sigma }^2}\sum _{n=0}^{N-1}(x[n]-A{)}^2\right]\text{\hspace{1em}(BP.8)}$$
It is apparent that the PDF is identical in form to the usual classical PDF $p(\mathbf{x};A)$. In the Bayesian case, however, the PDF is a conditional PDF, hence the “$|$” separator, while in the classical case, it represents an unconditional PDF, albeit parameterized by $A$, hence the separator "$;$"

Using $(BP.7)$ and $(BP.8)$, the posterior PDF becomes
$$p(A|\mathbf{x})=\{\begin{array}{rlrl}& \frac{\frac{1}{2A_0{\left(2\pi {\sigma }^2\right)}^{N/2}}\exp \left[-\frac{1}{2{\sigma }^2}{\sum }_{n=0}^{N-1}(x[n]-A{)}^2\right]}{{\int }_{-A_0}^{A_0}\frac{1}{2A_0{\left(2\pi {\sigma }^2\right)}^{N/2}}\exp \left[-\frac{1}{2{\sigma }^2}{\sum }_{n=0}^{N-1}(x[n]-A{)}^2\right]dA}& & |A|\le A_0\\ & 0& & |A|>A_0\end{array}$$
But
$$\begin{array}{rl}\sum _{n=0}^{N-1}(x[n]-A{)}^2& =\sum _{n=0}^{N-1}{x}^2[n]-2NA\bar{x}+N{A}^2=N(A-\bar{x}{)}^2+\sum _{n=0}^{N-1}{x}^2[n]-N{\bar{x}}^2\end{array}$$
Cancel out the terms that are irrelevant to $A$ and write the form as truncated Gaussian distribution
$$p(A|\mathbf{x})=\{\begin{array}{rlrl}& \frac{\frac{1}{\sqrt{2\pi {\sigma }^2/N}}\exp \left[-\frac{1}{2{\sigma }^2/N}(A-\bar{x}{)}^2\right]}{{\int }_{-A_0}^{A_0}\frac{1}{\sqrt{2\pi {\sigma }^2/N}}\exp \left[-\frac{1}{2{\sigma }^2/N}(A-\bar{x}{)}^2\right]dA}& & |A|\le A_0\\ & 0,& & |A|>A_0\end{array}\text{\hspace{1em}(BP.9)}$$
The final MMSE estimate is then given by
$$\hat{A}={\int }_{-\infty }^{\infty }Ap(A|\mathbf{x})dA=\frac{{\int }_{-A_0}^{A_0}A\frac{1}{\sqrt{2\pi {\sigma }^2/N}}\exp \left[-\frac{1}{2{\sigma }^2/N}(A-\bar{x}{)}^2\right]dA}{{\int }_{-A_0}^{A_0}\frac{1}{\sqrt{2\pi {\sigma }^2/N}}\exp \left[-\frac{1}{2{\sigma }^2/N}(A-\bar{x}{)}^2\right]dA}\text{\hspace{1em}(BP.10)}$$

• In absence of data (if $x$ gives no information about $A$),
  $$\hat{A}=E(A|\mathbf{x})=E(A)=0$$

• If $A_0\gg \sqrt{{\sigma }^2/N}$ (no truncation), e.g., as $N$ increases, $\hat{A}\to \bar{x}$
  $$\hat{A}=\underset{A_0\to \infty }{\lim }\frac{{\int }_{-A_0}^{A_0}A\frac{1}{\sqrt{2\pi {\sigma }^2/N}}\exp \left[-\frac{1}{2{\sigma }^2/N}(A-\bar{x}{)}^2\right]dA}{{\int }_{-A_0}^{A_0}\frac{1}{\sqrt{2\pi {\sigma }^2/N}}\exp \left[-\frac{1}{2{\sigma }^2/N}(A-\bar{x}{)}^2\right]dA}=\frac{\bar{x}}{1}=\bar{x}$$

The effect of the data is to position the posterior mean between $A=0$ and $A=\bar{x}$ in a compromise between the prior knowledge and that contributed by the data. To further appreciate this weighting consider what happens as $N$ becomes large so that the data knowledge becomes more important. As shown in Figure 10.4, as $N$ increases, we have from $(BP.9)$ that the posterior PDF becomes more concentrated about $\bar{x}$ (since ${\sigma }^2/N$ decrease). Hence, it becomes nearly Gaussian, and its mean becomes just $\bar{x}$. The MMSE estimator relies less and less on the prior knowledge and more on the data. It is said that the data “swamps out” the prior knowledge.

Choosing a Prior PDF

As shown in the previous section, once a prior PDF has been chosen, the MMSE estimator follows directly from $(BP.10)$. The only practical stumbling block that remains, however, is whether or not $E(\theta |\mathbf{x})$ can be determined in closed form. This motivates us to assume the Gaussian prior PDF instead of the uniform prior PDF
$$p(A)=\frac{1}{\sqrt{2\pi {\sigma }_A^2}}\exp \left[-\frac{1}{2{\sigma }_A^2}(A-{\mu }_A{)}^2\right]$$
With ${\mu }_A=0$ and $3{\sigma }_A=A_0$ the Gaussian prior PDF could be thought of as incorporating the knowledge that $|A|\le A_0$. From $(BP.6)$ and $(BP.8)$, we can formulate
$$\begin{array}{rl}p(A|\mathbf{x})=& \frac{p(\mathbf{x}|A)p(A)}{\int p(\mathbf{x}|A)p(A)dA}\\ =& \frac{\frac{1}{{\left(2\pi {\sigma }^2\right)}^{N/2}\sqrt{2\pi {\sigma }_A^2}}\exp \left[-\frac{1}{2{\sigma }^2}{\sum }_{n=0}^{N-1}{x}^2[n]\right]\exp \left[-\frac{1}{2{\sigma }^2}\left(N{A}^2-2NA\bar{x}\right)\right]\exp \left[-\frac{1}{2{\sigma }_A^2}{\left(A-{\mu }_A\right)}^2\right]}{{\int }_{-\infty }^{\infty }\frac{1}{{\left(2\pi {\sigma }^2\right)}^{N/2}\sqrt{2\pi {\sigma }_A^2}}\exp \left[-\frac{1}{2{\sigma }^2}{\sum }_{n=0}^{N-1}{x}^2[n]\right]\exp \left[-\frac{1}{2{\sigma }^2}\left(N{A}^2-2NA\bar{x}\right)\right]\exp \left[-\frac{1}{2{\sigma }_A^2}{\left(A-{\mu }_A\right)}^2\right]dA}\\ =& \frac{\exp \left[-\frac{1}{2}\left(\frac{1}{{\sigma }^2}(N{A}^2-2NA\bar{x})+\frac{1}{{\sigma }_A^2}{\left(A-{\mu }_A\right)}^2\right)\right]}{{\int }_{-\infty }^{\infty }\exp \left[-\frac{1}{2}\left(\frac{1}{{\sigma }^2}(N{A}^2-2NA\bar{x})+\frac{1}{{\sigma }_A^2}{\left(A-{\mu }_A\right)}^2\right)\right]dA}\\ =& \frac{\exp \left[-\frac{1}{2}Q(A)\right]}{{\int }_{-\infty }^{\infty }\exp \left[-\frac{1}{2}Q(A)\right]dA}\end{array}$$
Continuing, we have for $Q(A)$
$$\begin{array}{rl}Q(A)& =\frac{N}{{\sigma }^2}{A}^2-\frac{2NA\bar{x}}{{\sigma }^2}+\frac{A^2}{{\sigma }_A^2}-\frac{2{\mu }_{A}A}{{\sigma }_A^2}+\frac{{\mu }_A^2}{{\sigma }_A^2}\\ & =\left(\frac{N}{{\sigma }^2}+\frac{1}{{\sigma }_A^2}\right)A^2-2\left(\frac{N}{{\sigma }^2}\bar{x}+\frac{{\mu }_A}{{\sigma }_A^2}\right)A+\frac{{\mu }_A^2}{{\sigma }_A^2}\end{array}$$
Let
$$\begin{array}{rl}{\sigma }_{A\mid x}^2& =\frac{1}{\frac{N}{{\sigma }^2}+\frac{1}{{\sigma }_A^2}}\\ {\mu }_{A\mid x}& =\left(\frac{N}{{\sigma }^2}\bar{x}+\frac{{\mu }_A}{{\sigma }_A^2}\right){\sigma }_{A\mid x}^2\end{array}\text{\hspace{1em}(BP.11)}$$
Then, by completing the square we have
$$Q(A)=\frac{1}{{\sigma }_{A\mid x}^2}{\left(A-{\mu }_{A\mid x}\right)}^2-\frac{{\mu }_{A\mid x}^2}{{\sigma }_{A\mid x}^2}+\frac{{\mu }_A^2}{{\sigma }_A^2}$$
so that
$$\begin{array}{rl}p(A|\mathbf{x})& =\frac{\exp \left[-\frac{1}{2{\sigma }_{A\mid x}^2}{\left(A-{\mu }_{A\mid x}\right)}^2\right]\exp \left[-\frac{1}{2}\left(\frac{{\mu }_A^2}{{\sigma }_A^2}-\frac{{\mu }_{A\mid x}^2}{{\sigma }_{A\mid x}^2}\right)\right]}{{\int }_{-\infty }^{\infty }\exp \left[-\frac{1}{2{\sigma }_{A\mid x}^2}{\left(A-{\mu }_{A\mid x}\right)}^2\right]\exp \left[-\frac{1}{2}\left(\frac{{\mu }_A^2}{{\sigma }_A^2}-\frac{{\mu }_{A\mid x}^2}{{\sigma }_{A\mid x}^2}\right)\right]dA}\\ & =\frac{1}{\sqrt{2\pi {\sigma }_{A\mid x}^2}}\exp \left[-\frac{1}{2{\sigma }_{A\mid x}^2}{\left(A-{\mu }_{A\mid x}\right)}^2\right]\end{array}$$
where the last step follows from the requirement that $p(A|\mathbf{x})$ integrate to $1.$ The posterior PDF is also Gaussian, $p(A|\mathbf{x})\sim \mathcal{N}({\mu }_{A|x},{\sigma }_{A|x}^2)$.

In this form the MMSE estimator is readily found as
$$\begin{array}{rl}\hat{A}& =E(A|\mathbf{x})={\mu }_{A|x}=\frac{\frac{N}{{\sigma }^2}\bar{x}+\frac{{\mu }_A}{{\sigma }_A^2}}{\frac{N}{{\sigma }^2}+\frac{1}{{\sigma }_A^2}}\\ & =\frac{{\sigma }_A^2}{{\sigma }_A^2+\frac{{\sigma }^2}{N}}\bar{x}+\frac{\frac{{\sigma }^2}{N}}{{\sigma }_A^2+\frac{{\sigma }^2}{N}}{\mu }_A\\ & =\alpha \bar{x}+(1-\alpha ){\mu }_A\end{array}\text{\hspace{1em}(BP.12)}$$
where
$$\alpha =\frac{{\sigma }_A^2}{{\sigma }_A^2+\frac{{\sigma }^2}{N}}.$$
Note that $\alpha$ is a weighting factor since $0<\alpha <1$. It is interesting to examine the interplay between the prior knowledge and the data.

• If there is little data so that ${\sigma }_A^2\ll {\sigma }^2/N$, $\alpha$ is small and $\hat{A}\approx {\mu }_A$
• If more data are observed so that ${\sigma }_A^2\gg {\sigma }^2/N$, $\alpha \approx 1$ and $\hat{A}\approx \bar{x}$.

Alternatively, we may view this process by examining the posterior PDF. As $N$ increases,
$$\mathrm{var}(A|\mathbf{x})={\sigma }_{A\mid x}^2=\frac{1}{\frac{N}{{\sigma }^2}+\frac{1}{{\sigma }_A^2}}$$
will decrease. Also, the posterior mean or $\hat{A}$ $(BP.12)$ will approach $\bar{x}$.

Observe that if there is no prior knowledge, which can be modeled by letting ${\sigma }_A^2\to \infty$, then $\hat{A}\to \bar{x}$ for any data record length. The classical estimator is obtained.

Finally, we are going to prove that by using prior knowledge we could improve the estimation accuracy. Recall that
$$\mathrm{Bmse}(\hat{A})=\int \left[\int (A-\hat{A}{)}^{2}p(A|\mathbf{x})dA\right]p(\mathbf{x})d\mathbf{x}\text{\hspace{1em}(BP.4)}$$
Since $\hat{A}=E(A|\mathbf{x})$, we have
$$\begin{array}{rl}\mathrm{Bmse}(\hat{A})& =\int \left[\int (A-E(A|\mathbf{x}){)}^{2}p(A|\mathbf{x})dA\right]p(\mathbf{x})d\mathbf{x}\\ & =\int \mathrm{var}(A|\mathbf{x})p(\mathbf{x})d\mathbf{x}\\ & =\int \frac{1}{\frac{N}{{\sigma }^2}+\frac{1}{{\sigma }_A^2}}p(\mathbf{x})d\mathbf{x}\\ & =\frac{1}{\frac{N}{{\sigma }^2}+\frac{1}{{\sigma }_A^2}}\end{array}\text{\hspace{1em}(BP.13)}$$
This can be rewritten as
$$\mathrm{Bmse}(\hat{A})=\frac{{\sigma }^2}{N}(\frac{{\sigma }_A^2}{{\sigma }_A^2+\frac{{\sigma }^2}{N}})<\frac{{\sigma }^2}{N}\text{\hspace{1em}(BP.14)}$$
where ${\sigma }^2/N$ is the minimum MSE obtained when no prior knowledge is available (let ${\sigma }_A^2\to \infty$). Clearly, any prior knowledge when modeled in the Bayesian sense will improve our Bayesian estimator.

Properties of the Gaussian PDF

We now generalize the results of the previous section by examining the properties of the Gaussian PDF. The bivariate Gaussian PDF is first investigated to illustrate the important properties. Then, the corresponding results for the general multivariate Gaussian PDF are described.

Consider a jointly Gaussian random vector $[xy{]}^T$ whose PDF is
$$p(x,y)=\frac{1}{2\pi {\det }^{\frac{1}{2}}(\mathbf{C})}\exp \left[-\frac{1}{2}{\left[\begin{array}{c}x-E(x)\\ y-E(y)\end{array}\right]}^T{\mathbf{C}}^{-1}\left[\begin{array}{l}x-E(x)\\ y-E(y)\end{array}\right]\right]\text{\hspace{1em}(BF.15)}$$
This is also termed as bivariate Gaussian PDF. The mean vector and covariance matrix are
$$\begin{array}{rl}E\left(\left[\begin{array}{c}x\\ y\end{array}\right]\right)& =\left[\begin{array}{c}E(x)\\ E(y)\end{array}\right]\\ \mathbf{C}& =\left[\begin{array}{cc}\mathrm{var}(x)& \mathrm{cov}(x,y)\\ \mathrm{cov}(x,y)& \mathrm{var}(y)\end{array}\right]\end{array}$$
Note that the marginal PDFs $p(x)$ and $p(y)$ are also Gaussian, as can be verified by the integrations
$$\begin{array}{r}p(x)={\int }_{-\infty }^{\infty }p(x,y)dy=\frac{1}{\sqrt{2\pi \mathrm{var}(x)}}\exp \left[-\frac{1}{2\mathrm{var}(x)}(x-E(x){)}^2\right]\\ p(y)={\int }_{-\infty }^{\infty }p(x,y)dx=\frac{1}{\sqrt{2\pi \mathrm{var}(y)}}\exp \left[-\frac{1}{2\mathrm{var}(y)}(y-E(y){)}^2\right]\end{array}\text{\hspace{1em}(BF.16)}$$
The contours along which the PDF $p(x,y)$ is constant are those values of $x$ and $y$ for which
$${\left[\begin{array}{l}x-E(x)\\ y-E(y)\end{array}\right]}^T{\mathbf{C}}^{-1}\left[\begin{array}{l}x-E(x)\\ y-E(y)\end{array}\right]$$
is a constant. They are shown in Figure 10.6 as elliptical contours.

Once $x,$ say $x_0,$ is observed, the conditional PDF of $y$ becomes
$$p\left(y|x_0\right)=\frac{p\left(x_0,y\right)}{p\left(x_0\right)}=\frac{p\left(x_0,y\right)}{{\int }_{-\infty }^{\infty }p\left(x_0,y\right)dy}\text{\hspace{1em}(BF.17)}$$
Since from $(BF.15)$ the exponential argument is quadratic in $y$, $p(x_0,y)$ has the Gaussian form in $y$, and thus the conditional PDF $p(y|x_0)$ must also be Gaussian. To summarize,

Theorem 1 (Conditional PDF of Bivariate Gaussian) If $x$ and $y$ are distributed according to a bivariate Gaussian $PDF$ with mean vector $[E(x)E(y){]}^T$ and covariance matrix
$$\mathbf{C}=\left[\begin{array}{cc}\mathrm{var}(x)& \mathrm{cov}(x,y)\\ \mathrm{cov}(y,x)& \mathrm{var}(y)\end{array}\right]$$
so that
$$p(x,y)=\frac{1}{2\pi {\det }^{\frac{1}{2}}(\mathbf{C})}\exp \left[-\frac{1}{2}{\left[\begin{array}{l}x-E(x)\\ y-E(y)\end{array}\right]}^T{\mathbf{C}}^{-1}\left[\begin{array}{l}x-E(x)\\ y-E(y)\end{array}\right]\right]$$
then the conditional PDF $p(y|x)$ is also Gaussian and
$$E(y|x)=E(y)+\frac{\mathrm{cov}(x,y)}{\mathrm{var}(x)}(x-E(x))\text{\hspace{1em}(BF.18)}$$

$$\mathrm{var}(y|x)=\mathrm{var}(y)-\frac{{\mathrm{cov}}^2(x,y)}{\mathrm{var}(x)}\text{\hspace{1em}(BF.19)}$$

Assuming that $x$ and $y$ are not independent and hence $\mathrm{cov}(x,y)\ne 0,$ the posterior PDF becomes more concentrated since there is less uncertainty about $y.$ To verify this, note from $(BF.19)$ that
$$\begin{array}{rl}\mathrm{var}(y|x)& =\mathrm{var}(y)\left[1-\frac{{\mathrm{cov}}^2(x,y)}{\mathrm{var}(x)\mathrm{var}(y)}\right]\\ & =\mathrm{var}(y)\left(1-{\rho }^2\right)\end{array}\text{\hspace{1em}(BF.20)}$$
where
$$\rho =\frac{\mathrm{cov}(x,y)}{\sqrt{\mathrm{var}(x)\mathrm{var}(y)}}\text{\hspace{1em}(BF.21)}$$
is the correlation coefficient satisfying $|\rho |\le 1$. From our previous discussions, we also realize that $E(y|x)$ is the MMSE estimator of $y$ after observing $x$, so that from $(BF.18)$
$$\hat{y}=E(y)+\frac{\mathrm{cov}(x,y)}{\mathrm{var}(x)}(x-E(x))\text{\hspace{1em}(BF.22)}$$
In normalized form (a random variable with zero mean and unity variance) this becomes
$$\frac{\hat{y}-E(y)}{\sqrt{\mathrm{var}(y)}}=\frac{\mathrm{cov}(x,y)}{\sqrt{\mathrm{var}(x)\mathrm{var}(y)}}\frac{x-E(x)}{\sqrt{\mathrm{var}(x)}}$$
or
$${\hat{y}}_n=\rho x_n\text{\hspace{1em}(BF.23)}$$
If the random variables are already normalized $(E(x)=E(y)=0,\mathrm{var}(x)=\mathrm{var}(y)=1)$, the constant PDF contours appear as in Figure 10.7.

The locations of the peaks of $p(x,y)$, when considered as a function of $y$ for each $x$, is the dashed line $y=\rho x$, and it is readily shown that $\hat{y}=E(y|x)=\rho x$.

The MMSE estimator therefore exploits the correlation between the random variables to estimate the realization of one based on the realization of the other. The minimum MSE is, from $(BF.13)$ and $(BF.20)$,
$$\begin{array}{rl}\mathrm{Bmse}(\hat{y})=& \int \mathrm{var}(y|x)p(x)dx\\ =& \mathrm{var}(y|x)\\ =& \mathrm{var}(y)(1-{\rho }^2)\end{array}\text{\hspace{1em}(BF.24)}$$
Hence, the quality of our estimator also depends on the correlation coefficient, which is a measure of the statistical dependence between $x$ and $y$.

---

To generalize these results consider a jointly Gaussian vector ${\left[{\mathbf{x}}^T{\mathbf{y}}^T\right]}^T,$ where $\mathbf{x}$ is $k\times 1$ and $\mathbf{y}$ is $l\times 1.$ In other words, ${\left[{\mathbf{x}}^T{\mathbf{y}}^T\right]}^T$ is distributed according to a multivariate Gaussian PDF. Then, the conditional PDF of $\mathbf{y}$ for a given $\mathbf{x}$ is also Gaussian. as summarized in the following theorem (see Appendix 10A for proof).

Theorem 10.2 (Conditional PDF of Multivariate Gaussian) If $\mathbf{x}$ and $\mathbf{y}$ are jointly Gaussian, where $\mathbf{x}$ is $k\times 1$ and $\mathbf{y}$ is $l\times 1,$ with mean vector ${\left[E(\mathbf{x}{)}^{T}E(\mathbf{y}{)}^T\right]}^T$ and
partitioned covariance matrix
$$\mathbf{C}=\left[\begin{array}{ll}{\mathbf{C}}_{xx}& {\mathbf{C}}_{xy}\\ {\mathbf{C}}_{yx}& {\mathbf{C}}_{yy}\end{array}\right]=\left[\begin{array}{ll}k\times k& k\times l\\ l\times k& l\times l\end{array}\right]\text{\hspace{1em}(BF.25)}$$
so that
$$p(\mathbf{x},\mathbf{y})=\frac{1}{(2\pi {)}^{\frac{k+1}{2}}{\det }^{\frac{1}{2}}(\mathbf{C})}\exp \left[-\frac{1}{2}{\left(\left[\begin{array}{l}\mathbf{x}-E(\mathbf{x})\\ \mathbf{y}-E(\mathbf{y})\end{array}\right]\right)}^T{\mathbf{C}}^{-1}\left(\left[\begin{array}{l}\mathbf{x}-E(\mathbf{x})\\ \mathbf{y}-E(\mathbf{y})\end{array}\right]\right)\right]$$
Then the conditional PDF $p(\mathbf{y}|\mathbf{x})$ is also Gaussian and
$$E(\mathbf{y}|\mathbf{x})=E(\mathbf{y})+{\mathbf{C}}_{yx}{\mathbf{C}}_{xx}^{-1}(\mathbf{x}-E(\mathbf{x}))\text{\hspace{1em}(BF.26)}$$

$${\mathbf{C}}_{y|x}={\mathbf{C}}_{yy}-{\mathbf{C}}_{yx}{\mathbf{C}}_{xx}^{-1}{\mathbf{C}}_{xy}\text{\hspace{1em}(BF.27)}$$

Bayesian Linear Model

Let the data be modeled as
$$\mathbf{x}=\mathbf{H}\mathbf{\theta }+\mathbf{w}\text{\hspace{1em}(BF.28)}$$
where $\mathbf{x}$ is an $N\times 1$ data vector, $\mathbf{H}$ is a known $N\times p$ matrix, $\mathbf{\theta }$ is a $p\times 1$ random vector with prior PDF $\mathcal{N}({\mathbf{\mu }}_{\theta },{\mathbf{C}}_{\theta })$, and $\mathbf{w}$ is a $N\times 1$ noise vector with PDF $\mathcal{N}(0,{\mathbf{C}}_w)$ and independent of $\mathbf{\theta }$. This data model is termed the Bayesian general linear model. It differs from the classical general linear model in that $\mathbf{\theta }$ is modeled as a random variable with a Gaussian prior PDF.

From Theorem 10.2 we know that if $\mathbf{x}$ and $\mathbf{\theta }$ are jointly Gaussian, then the posterior PDF is also Gaussian. Hence, it only remains to verify that this is indeed the case. Let $\mathbf{z}=[{\mathbf{x}}^T{\mathbf{\theta }}^T{]}^T$, so that from $(BF.28)$ we have
$$\mathbf{z}=\left[\begin{array}{c}\mathbf{x}\\ \mathbf{w}\end{array}\right]=\left[\begin{array}{c}\mathbf{H}\mathbf{\theta }+\mathbf{w}\\ \mathbf{w}\end{array}\right]=\left[\begin{array}{cc}\mathbf{H}& \mathbf{I}\\ \mathbf{I}& 0\end{array}\right]\left[\begin{array}{c}\mathbf{\theta }\\ \mathbf{w}\end{array}\right]$$
Since $\mathbf{\theta }$ and $\mathbf{w}$ are independent of each other and each one is Gaussian, they are jointly Gaussian. Furthermore, because $\mathbf{z}$ is a linear transformation of a Gaussian vector, it too is Gaussian. Hence, Theorem 10.2 applies directly, and we need only determine the mean and covariance of the posterior PDF.

We can obtain the means and covariances
$$\begin{array}{rl}E(\mathbf{x})=& E(\mathbf{H}\mathbf{\theta }+\mathbf{w})=\mathbf{H}E(\mathbf{\theta })=\mathbf{H}{\mathbf{\mu }}_{\theta }\\ E(\mathbf{y})=& E(\mathbf{\theta })={\mathbf{\mu }}_{\theta }\\ {\mathbf{C}}_{xx}=& E[(\mathbf{x}-E(\mathbf{x}))(\mathbf{x}-E(\mathbf{x}){)}^T]\\ =& E[(\mathbf{H}(\mathbf{\theta }-{\mathbf{\mu }}_{\theta })+\mathbf{w})(\mathbf{H}(\mathbf{\theta }-{\mathbf{\mu }}_{\theta })+\mathbf{w}{)}^T]\\ =& \mathbf{H}{\mathbf{C}}_{\theta }{\mathbf{H}}^T+{\mathbf{C}}_w\\ {\mathbf{C}}_{yx}=& E[(\mathbf{y}-E(\mathbf{y}))(\mathbf{x}-E(\mathbf{x}){)}^T]\\ =& E[(\mathbf{\theta }-{\mathbf{\mu }}_{\theta })(\mathbf{H}(\mathbf{\theta }-{\mathbf{\mu }}_{\theta })+\mathbf{w}{)}^T]\\ =& {\mathbf{C}}_{\theta }{\mathbf{H}}^T\end{array}$$
We can now summarize our results for the Bayesian general linear model.

Theorem 10.3 (Posterior PDF for the Bayesian General Linear Model) If the observed data $\mathbf{x}$ can be modeled as
$$\mathbf{x}=\mathbf{H}\mathbf{\theta }+\mathbf{w}$$
where $\mathbf{x}$ is an $N\times 1$ data vector, $\mathbf{H}$ is a known $N\times p$ matrix, $\mathbf{\theta }$ is a $p\times 1$ random vector with prior PDF $\mathcal{N}\left({\mu }_{\theta },{\mathbf{C}}_{\theta }\right),$ and $\mathbf{w}$ is an $N\times 1$ noise vector with PDF $\mathcal{N}\left(0,{\mathbf{C}}_w\right)$ and independent of $\mathbf{\theta },$ then the posterior PDF $p(\mathbf{\theta }|\mathbf{x})$ is Gaussian with mean
$$E(\mathbf{\theta }|\mathbf{x})={\mathbf{\mu }}_{\theta }+{\mathbf{C}}_{\theta }{\mathbf{H}}^T{\left(\mathbf{H}{\mathbf{C}}_{\theta }{\mathbf{H}}^T+{\mathbf{C}}_w\right)}^{-1}\left(\mathbf{x}-\mathbf{H}{\mathbf{\mu }}_{\theta }\right)\text{\hspace{1em}(BF.29)}$$
and covariance
$${\mathbf{C}}_{\theta \mid x}={\mathbf{C}}_{\theta }-{\mathbf{C}}_{\theta }{\mathbf{H}}^T{\left(\mathbf{H}{\mathbf{C}}_{\theta }{\mathbf{H}}^T+{\mathbf{C}}_w\right)}^{-1}\mathbf{H}{\mathbf{C}}_{\theta }\text{\hspace{1em}(BF.30)}$$
In contrast to the classical general linear model, $\mathbf{H}$ need not be full rank to ensure the invertibility of ${\mathbf{H}\mathbf{C}}_{\theta }{\mathbf{H}}^T+{\mathbf{C}}_w.$

Alternative formulation using Matrix inversion lemma:
$$E(\mathbf{\theta }|\mathbf{x})={\mathbf{\mu }}_{\theta }+{\left({\mathbf{C}}_{\theta }^{-1}+{\mathbf{H}}^T{\mathbf{C}}_w^{-1}\mathbf{H}\right)}^{-1}{\mathbf{H}}^T{\mathbf{C}}_w^{-1}\left(\mathbf{x}-\mathbf{H}{\mathbf{\mu }}_{\theta }\right)\text{\hspace{1em}(BF.31)}$$

$${\mathbf{C}}_{\theta \mid x}={\left({\mathbf{C}}_{\theta }^{-1}+{\mathbf{H}}^T{\mathbf{C}}_w^{-1}\mathbf{H}\right)}^{-1}\text{\hspace{1em}(BF.32)}$$

We illustrate the use of these formulas by applying them to Example DC Level in WGN:

Let us assume now that the prior distribution of $A$ is Gaussian: $A\sim \mathcal{N}\left({\mu }_A,{\sigma }_A^2\right)$, and $w[n]$ is white Gaussian noise, i.e., for $n=0,\dots ,N-1$, $w[n]\sim \mathcal{N}\left(0,{\sigma }^2\right)$
$$\mathbf{x}=1A+\mathbf{w}$$
then, $\mathbf{x}$ and $A$ are jointly Gaussian $(k=N$ and $l=1),$ with mean
$$\begin{array}{r}E(A|\mathbf{x})={\mu }_A+{\sigma }_A^2{1}^T(1{\sigma }_A^2{1}^T+{\sigma }^{2}I{)}^{-1}(\mathbf{x}-1{\mu }_A)\end{array}$$
Using Woodbury’s identity
$${\left(\mathbf{I}+\frac{{\sigma }_A^2}{{\sigma }^2}11^T\right)}^{-1}=\mathbf{I}-\frac{\frac{{\sigma }_A^2}{{\sigma }^2}11^T}{1+N\frac{{\sigma }_A^2}{{\sigma }^2}}$$
so that
$$\begin{array}{rl}\hat{A}=E(A|\mathbf{x})& ={\mu }_A+\frac{{\sigma }_A^2}{{\sigma }^2}{1}^T\left(\mathbf{I}-\frac{11^T}{N+\frac{{\sigma }^2}{{\sigma }_A^2}}\right)\left(\mathbf{x}-1{\mu }_A\right)\\ & ={\mu }_A+\frac{{\sigma }_A^2}{{\sigma }^2}\left(1^T-\frac{N}{N+\frac{{\sigma }^2}{{\sigma }_A^2}}{1}^T\right)\left(\mathbf{x}-1{\mu }_A\right)\\ & ={\mu }_A+\frac{{\sigma }_A^2}{{\sigma }^2}\left(1-\frac{N}{N+\frac{{\sigma }^2}{{\sigma }_A^2}}\right)\left(N\bar{x}-N{\mu }_A\right)\\ & ={\mu }_A+\frac{{\sigma }_A^2}{{\sigma }_A^2+\frac{{\sigma }^2}{N}}\left(\bar{x}-{\mu }_A\right)\\ & =\frac{{\sigma }_A^2}{{\sigma }_A^2+\frac{{\sigma }^2}{N}}\bar{x}+\frac{\frac{{\sigma }^2}{N}}{{\sigma }_A^2+\frac{{\sigma }^2}{N}}{\mu }_A\\ & =\alpha \bar{x}+(1-\alpha ){\mu }_A\end{array}$$
which is exactly the same as $(BP.12)$. And
$$\begin{array}{rl}\mathrm{Bmse}(\hat{A})& =\int \left[\int (A-E(A|\mathbf{x}){)}^{2}p(A|\mathbf{x})dA\right]p(\mathbf{x})d\mathbf{x}\\ & =\int \mathrm{var}(A|\mathbf{x})p(\mathbf{x})d\mathbf{x}=\mathrm{var}(A|\mathbf{x})\\ & ={\sigma }_A^2-\frac{{\sigma }_A^2}{{\sigma }^2}{1}^T\left(\mathbf{I}-\frac{11^T}{N+\frac{{\sigma }^2}{{\sigma }_A^2}}\right)1{\sigma }_A^2\\ & =\frac{1}{\frac{N}{{\sigma }^2}+\frac{1}{{\sigma }_A^2}}\end{array}$$

which is just $(BP.13)$.