Chapter 6 Exercises and Solutions

Thirteen exercises, worked in full. Every distribution, posterior and covariance below is checked numerically — against a grid, a simulation, or both — and the gaps are reported.

Four of them (6.6, 6.7, 6.11, 6.13) are proofs rather than computations. Each one is short, each is used earlier in the chapter, and each is worth doing rather than skipping.

How to use this page

6.1 and 6.2 are drills on §6.2 and §6.4 — do them by hand. 6.3 and 6.4 are the two flavours of Bayes: conjugate updating and a discrete two-hypothesis problem. 6.5 and 6.12 are the payoff: between them they derive the Kalman filter and the Gaussian-linear-model posterior using nothing but §6.5’s closure rules, and they are the reason the chapter spent so long on marginals and conditionals. 6.6 to 6.11 and 6.13 fill in the results the chapter used without proof.

A convention used throughout

Where a covariance shrinks I report the trace before and after, plus the eigenvalues of the difference. Conditioning must produce a positive semidefinite difference (§6.5’s Equation 6.67), so a negative eigenvalue there would mean an error. Every such check below passes with the smallest eigenvalue at exactly $0$ — which is itself informative: one direction learned nothing.

6.1 — Marginals and conditionals of a joint table

Given $p(x,y)$ over $X$ with five states and $Y$ with three:

	$x_1$	$x_2$	$x_3$	$x_4$	$x_5$	$p(y)$
$y_1$	$0.01$	$0.02$	$0.03$	$0.10$	$0.10$	$0.26$
$y_2$	$0.05$	$0.10$	$0.05$	$0.07$	$0.20$	$0.47$
$y_3$	$0.10$	$0.05$	$0.03$	$0.05$	$0.04$	$0.27$
$p(x)$	$0.16$	$0.17$	$0.11$	$0.22$	$0.34$	$1.00$

a. The marginals, by the sum rule (Equation 6.20) — add along the axis you are removing:

$$p(x) = \begin{bmatrix}0.16 & 0.17 & 0.11 & 0.22 & 0.34\end{bmatrix}, \qquad p(y) = \begin{bmatrix}0.26 & 0.47 & 0.27\end{bmatrix}$$

Both sum to $1.0000000000$, and the whole table sums to $1$ — so it is a valid joint.

b. The conditionals, by Equations 6.13 and 6.14 — divide by the relevant total.

Conditioning on $Y = y_1$ fixes a row and divides by that row’s total, $p(y_1) = 0.26$:

$$p(x \mid Y=y_1) = \frac{1}{0.26}\begin{bmatrix}0.01 & 0.02 & 0.03 & 0.10 & 0.10\end{bmatrix} = \begin{bmatrix}0.038462 & 0.076923 & 0.115385 & 0.384615 & 0.384615\end{bmatrix}$$

Conditioning on $X = x_3$ fixes a column and divides by $p(x_3) = 0.11$:

$$p(y \mid X=x_3) = \frac{1}{0.11}\begin{bmatrix}0.03 \\ 0.05 \\ 0.03\end{bmatrix} = \begin{bmatrix}0.272727 \\ 0.454545 \\ 0.272727\end{bmatrix}$$

Each sums to $1.0000000000$ — the check that catches a transposed conditional.

Note how different the two are in character: given $y_1$, the mass is concentrated on $x_4$ and $x_5$ ($0.385$ each); given $x_3$, $Y$ is nearly uniform. Conditioning on one variable tells you a lot here; conditioning on the other tells you little.

6.2 — The bimodal mixture

$$0.4\,\mathcal{N}\!\left(\begin{bmatrix}10\\2\end{bmatrix},\begin{bmatrix}1&0\\0&1\end{bmatrix}\right) + 0.6\,\mathcal{N}\!\left(\begin{bmatrix}0\\0\end{bmatrix},\begin{bmatrix}8.4&2.0\\2.0&1.7\end{bmatrix}\right)$$

a. The marginals. Two facts compose. A Gaussian’s marginal is Gaussian with the corresponding block (Equation 6.68), and marginalising is linear, so it passes through the mixture weights. So each marginal is a one-dimensional mixture of the marginals of the components:

$$p(x_1) = 0.4\,\mathcal{N}(10,\ 1) + 0.6\,\mathcal{N}(0,\ 8.4)$$ $$p(x_2) = 0.4\,\mathcal{N}(2,\ 1) + 0.6\,\mathcal{N}(0,\ 1.7)$$

Note the off-diagonal $2.0$ disappears from both — marginalising discards exactly the cross-covariance information.

b. Mean, mode and median of each marginal. The means are exact by linearity (Equation 6.34); the modes and medians come from the density on a fine grid.

	mean (exact)	mean (measured)	mode	median	modes found
$x_1$	$4.0$	$4.000008$	$9.998390$	$2.803905$	2
$x_2$	$0.8$	$0.800000$	$1.330770$	$0.868710$	1

Three different answers per marginal, and the two marginals do not even agree on how many humps there are — $p(x_1)$ is bimodal, $p(x_2)$ is unimodal. That is the book’s point in Figure 6.4, and it is why the mode of $p(x_2)$, $1.3308$, is at neither component mean ($2$ and $0$): the two humps have merged into one whose peak sits between them.

c. The two-dimensional mean and mode.

$$\mathbb{E}[\mathbf{x}] = 0.4\begin{bmatrix}10\\2\end{bmatrix} + 0.6\begin{bmatrix}0\\0\end{bmatrix} = \begin{bmatrix}4.0\\0.8\end{bmatrix}$$

exactly. The mode, by grid search, is at $(9.9958,\ 1.9942)$ — essentially the first component’s centre, $(10, 2)$, despite that component carrying only $0.4$ of the weight. Measured densities: $0.063735$ at $(10,2)$ against $0.029783$ at $(0,0)$. The first component wins because it is far tighter (covariance $\mathbf{I}$ against a much broader block), so its mass is concentrated into a taller peak.

And the mean is not near either mode. $(4.0, 0.8)$ sits in the low-density gap between the two clusters — a point the distribution almost never produces.

6.3 — A conjugate prior for the compiler

A program compiles or does not, modelled as $p(x\mid\mu) = \mu^x(1-\mu)^{1-x}$.

Choose the prior. Table 6.2: the conjugate prior for a Bernoulli likelihood is the Beta. Take $\mu \sim \operatorname{Beta}(\alpha,\beta)$.

Compute the posterior. For $N$ independent runs with $s = \sum_n x_n$ successes, the likelihood is $\mu^{s}(1-\mu)^{N-s}$, so

$$p(\mu \mid x_1,\ldots,x_N) \propto \mu^{s}(1-\mu)^{N-s}\cdot\mu^{\alpha-1}(1-\mu)^{\beta-1} = \mu^{\alpha+s-1}(1-\mu)^{\beta+N-s-1}$$ $$\boxed{p(\mu \mid x_1,\ldots,x_N) = \operatorname{Beta}(\alpha+s,\ \beta+N-s)}$$

The successes go into $\alpha$ and the failures into $\beta$ — pseudo-counts again.

Checked: with $\operatorname{Beta}(2,2)$ and a simulated $34$ successes in $50$ runs, the posterior is $\operatorname{Beta}(36,18)$ with mean $0.666667$. Against a $200\,001$-point grid computation: worst density gap $8.8\times10^{-14}$.

Note the whole run history collapses to two numbers. You never need to keep the individual $x_n$.

6.4 — Which bag did the mango come from?

Bag 1: four mangos, two apples. Bag 2: four mangos, four apples. The coin shows heads with probability $0.6$; heads means bag 1.

Set it up as Bayes. The prior is the coin: $P(\text{bag 1}) = 0.6$, $P(\text{bag 2}) = 0.4$. The likelihoods are the bag compositions:

$$P(\text{mango}\mid\text{bag 1}) = \tfrac46 = \tfrac23, \qquad P(\text{mango}\mid\text{bag 2}) = \tfrac48 = \tfrac12$$

The evidence (Equation 6.27):

$$P(\text{mango}) = 0.6\cdot\tfrac23 + 0.4\cdot\tfrac12 = 0.4 + 0.2 = 0.6$$

The posterior:

$$P(\text{bag 2}\mid\text{mango}) = \frac{0.4 \cdot \tfrac12}{0.6} = \frac{0.2}{0.6} = \boxed{\tfrac13 \approx 0.3333}$$

and $P(\text{bag 1}\mid\text{mango}) = \tfrac23$. Simulated over four million trials: $0.333517$.

Read the direction of the update. The prior on bag 2 was $0.4$; the posterior is $0.333$. Seeing a mango made bag 2 less likely, because bag 1 is the mango-richer bag ($\tfrac23$ against $\tfrac12$). Both the prior and the likelihood point the same way here, so they reinforce.

6.5 — The Kalman filter, from §6.5 alone

$$\mathbf{x}_{t+1} = \mathbf{A}\mathbf{x}_t + \mathbf{w},\quad \mathbf{w}\sim\mathcal{N}(\mathbf{0},\mathbf{Q}), \qquad \mathbf{y}_t = \mathbf{C}\mathbf{x}_t + \mathbf{v},\quad \mathbf{v}\sim\mathcal{N}(\mathbf{0},\mathbf{R})$$

with $p(\mathbf{x}_0) = \mathcal{N}(\boldsymbol{\mu}_0,\boldsymbol{\Sigma}_0)$.

a. What is the form of $p(\mathbf{x}_0,\ldots,\mathbf{x}_T)$?

Gaussian. The justification is two closure properties and nothing else. $\mathbf{A}\mathbf{x}_t$ is an affine map of a Gaussian, so Gaussian (Equation 6.88). Adding independent Gaussian noise keeps it Gaussian (Equation 6.78). So each $\mathbf{x}_{t+1}$ is Gaussian given $\mathbf{x}_t$, and the whole trajectory is jointly Gaussian — it is built by repeatedly applying operations that never leave the family.

b. Assume $p(\mathbf{x}_t\mid\mathbf{y}_1,\ldots,\mathbf{y}_t) = \mathcal{N}(\boldsymbol{\mu}_t,\boldsymbol{\Sigma}_t)$.

1. The prediction. Push the current belief through the dynamics:

$$p(\mathbf{x}_{t+1}\mid\mathbf{y}_{1:t}) = \mathcal{N}\bigl(\mathbf{A}\boldsymbol{\mu}_t,\ \mathbf{A}\boldsymbol{\Sigma}_t\mathbf{A}^\top + \mathbf{Q}\bigr)$$

Equation 6.88 for the map, Equation 6.78 for the noise. Note the covariance grows — prediction adds uncertainty.

2. The joint with the next observation. $\mathbf{y}_{t+1} = \mathbf{C}\mathbf{x}_{t+1} + \mathbf{v}$ is another affine map, so stack them:

$$p(\mathbf{x}_{t+1},\mathbf{y}_{t+1}\mid\mathbf{y}_{1:t}) = \mathcal{N}\!\left(\begin{bmatrix}\mathbf{A}\boldsymbol{\mu}_t\\\mathbf{C}\mathbf{A}\boldsymbol{\mu}_t\end{bmatrix},\ \begin{bmatrix}\mathbf{P} & \mathbf{P}\mathbf{C}^\top\\ \mathbf{C}\mathbf{P} & \mathbf{C}\mathbf{P}\mathbf{C}^\top+\mathbf{R}\end{bmatrix}\right)$$

writing $\mathbf{P} = \mathbf{A}\boldsymbol{\Sigma}_t\mathbf{A}^\top + \mathbf{Q}$. The cross-covariance is $\operatorname{Cov}[\mathbf{x},\mathbf{y}] = \mathbf{P}\mathbf{C}^\top$ by Equation 6.52.

3. Observe $\mathbf{y}_{t+1} = \hat{\mathbf{y}}$ and condition. Equations 6.66 and 6.67, verbatim:

$$\boldsymbol{\mu}_{t+1} = \mathbf{A}\boldsymbol{\mu}_t + \underbrace{\mathbf{P}\mathbf{C}^\top(\mathbf{C}\mathbf{P}\mathbf{C}^\top+\mathbf{R})^{-1}}_{\text{the Kalman gain }\mathbf{K}}\bigl(\hat{\mathbf{y}} - \mathbf{C}\mathbf{A}\boldsymbol{\mu}_t\bigr)$$ $$\boldsymbol{\Sigma}_{t+1} = \mathbf{P} - \mathbf{K}\,\mathbf{C}\mathbf{P}$$

That is the Kalman filter. No new mathematics appeared: predict with 6.78 and 6.88, update with 6.66 and 6.67.

Worked with $\mathbf{A} = \begin{bmatrix}0.9&0.2\\-0.1&0.95\end{bmatrix}$, $\mathbf{C} = \begin{bmatrix}1&0.5\end{bmatrix}$ and an observation $\hat{y} = 0.8$:

	value
predicted mean	$(0.800,\ -0.575)$
predicted trace	$1.127750$
Kalman gain	$(0.629390,\ 0.321205)$
posterior mean	$(0.980950,\ -0.482653)$
posterior trace	$0.533358$
eigenvalues of $\mathbf{P}-\boldsymbol{\Sigma}_{t+1}$	$[0,\ 0.594392]$

The trace more than halved, and the difference is positive semidefinite as Equation 6.67 requires. The zero eigenvalue is meaningful: with a scalar observation, only one direction in the two-dimensional state can be informed. The other learned nothing.

6.6 — Prove Equation 6.44

Expand the definition and use linearity:

$$\mathbb{V}[x] = \mathbb{E}\bigl[(x-\mu)^2\bigr] = \mathbb{E}\bigl[x^2 - 2\mu x + \mu^2\bigr] = \mathbb{E}[x^2] - 2\mu\,\mathbb{E}[x] + \mu^2$$

Now $\mathbb{E}[x] = \mu$, so the last two terms are $-2\mu^2 + \mu^2 = -\mu^2$:

$$\boxed{\mathbb{V}[x] = \mathbb{E}[x^2] - \bigl(\mathbb{E}[x]\bigr)^2}$$

The only step is linearity of expectation (Equation 6.34), and the only subtlety is remembering that $\mu$ is a constant so it comes out of the expectation. Measured on eight numbers: $4.000000000000$ both ways, gap $0.0\text{e}{+}00$.

6.7 — Prove Equation 6.45

Expand the square inside the double sum:

$$\frac{1}{N^2}\sum_{i,j}(x_i-x_j)^2 = \frac{1}{N^2}\sum_{i,j}\bigl(x_i^2 - 2x_ix_j + x_j^2\bigr)$$

Take the three pieces separately. $\sum_{i,j}x_i^2 = N\sum_i x_i^2$ (the inner sum over $j$ contributes $N$ copies), and likewise $\sum_{i,j}x_j^2 = N\sum_j x_j^2$. The cross term factorises: $\sum_{i,j}x_ix_j = \bigl(\sum_i x_i\bigr)^2$. So

$$= \frac{1}{N^2}\left(2N\sum_i x_i^2 - 2\Bigl(\sum_i x_i\Bigr)^{\!2}\right) = 2\left[\frac{1}{N}\sum_i x_i^2 - \left(\frac{1}{N}\sum_i x_i\right)^{\!2}\right]$$

which is twice the raw-score variance of Equation 6.44. Measured: $8.000000000000$ against $2\times4.000000000000$, gap $0.0\text{e}{+}00$ — and it took $N^2 = 64$ terms to compute what $N = 8$ terms give.

The geometric reading the book offers: the pairwise distances between points and the distances from their centre carry the same information.

6.8 — The Bernoulli in natural parameter form

Take logs and regroup:

$$\mu^x(1-\mu)^{1-x} = \exp\bigl[x\log\mu + (1-x)\log(1-\mu)\bigr] = \exp\!\left[x\log\frac{\mu}{1-\mu} + \log(1-\mu)\right]$$

Matching against $h(x)\exp\bigl(\langle\theta,\phi(x)\rangle - A(\theta)\bigr)$:

$$h(x) = 1, \qquad \phi(x) = x, \qquad \theta = \log\frac{\mu}{1-\mu}, \qquad A(\theta) = -\log(1-\mu) = \log(1+e^\theta)$$

Verified at $\mu\in\{0.25, 0.6\}$ and $x\in\{0,1\}$: worst gap $5.6\times10^{-17}$.

6.9 — Binomial, Beta, and their product

Binomial, at fixed $N$:

$$\binom{N}{m}\mu^m(1-\mu)^{N-m} = \binom{N}{m}\exp\!\left[m\log\frac{\mu}{1-\mu} + N\log(1-\mu)\right]$$

so $h(m) = \binom{N}{m}$, $\phi(m) = m$, $\theta = \log\frac{\mu}{1-\mu}$, and $A(\theta) = N\log(1+e^\theta)$ — the Bernoulli’s, scaled by $N$. Verified across all $13$ outcomes at $N=12$, $\mu=0.35$: worst gap $1.7\times10^{-16}$.

Beta. Here the random variable is $\mu$:

$$\frac{1}{B(\alpha,\beta)}\mu^{\alpha-1}(1-\mu)^{\beta-1} = \exp\bigl[(\alpha-1)\log\mu + (\beta-1)\log(1-\mu) - \log B(\alpha,\beta)\bigr]$$

so $h(\mu) = 1$, $\phi(\mu) = \begin{bmatrix}\log\mu\\\log(1-\mu)\end{bmatrix}$, $\boldsymbol{\theta} = \begin{bmatrix}\alpha-1\\\beta-1\end{bmatrix}$, and $A(\boldsymbol{\theta}) = \log B(\alpha,\beta)$. Verified over a grid at $\alpha=3$, $\beta=4$: worst gap $1.3\times10^{-15}$.

Their product is exponential-family too. Both are exponentials of linear functions of $(\log\mu,\ \log(1-\mu))$, so multiplying adds the exponents — which means adding the natural parameters:

$$\boldsymbol{\theta}_{\text{product}} = \begin{bmatrix}\alpha-1\\\beta-1\end{bmatrix} + \begin{bmatrix}m\\N-m\end{bmatrix}$$

still of the same form. That is conjugacy, restated. The posterior is $\operatorname{Beta}(\alpha+m,\ \beta+N-m)$ because natural parameters add — the same fact Example 6.11 derived by staring at the exponents.

6.10 — The product of two Gaussians, two ways

a. By completing the square. The product’s exponent is $-\tfrac12\bigl[(\mathbf{x}-\mathbf{a})^\top\mathbf{A}^{-1}(\mathbf{x}-\mathbf{a}) + (\mathbf{x}-\mathbf{b})^\top\mathbf{B}^{-1}(\mathbf{x}-\mathbf{b})\bigr]$. Collect powers of $\mathbf{x}$: the quadratic term is $\mathbf{x}^\top(\mathbf{A}^{-1}+\mathbf{B}^{-1})\mathbf{x}$, so the new precision is $\mathbf{A}^{-1}+\mathbf{B}^{-1}$ and hence

$$\mathbf{C} = (\mathbf{A}^{-1}+\mathbf{B}^{-1})^{-1}$$

The linear term is $-2\mathbf{x}^\top(\mathbf{A}^{-1}\mathbf{a}+\mathbf{B}^{-1}\mathbf{b})$, and matching it against $-2\mathbf{x}^\top\mathbf{C}^{-1}\mathbf{c}$ gives

$$\mathbf{c} = \mathbf{C}(\mathbf{A}^{-1}\mathbf{a}+\mathbf{B}^{-1}\mathbf{b})$$

Everything left over is free of $\mathbf{x}$, and collecting it gives Equation 6.76’s constant.

b. In exponential-family form. From Example 6.13, a Gaussian’s natural parameters are $\bigl(\boldsymbol{\Sigma}^{-1}\boldsymbol{\mu},\ -\tfrac12\boldsymbol{\Sigma}^{-1}\bigr)$. Multiplying two densities adds natural parameters (as in 6.9), so immediately

$$\boldsymbol{\Sigma}_{\text{new}}^{-1} = \mathbf{A}^{-1}+\mathbf{B}^{-1}, \qquad \boldsymbol{\Sigma}_{\text{new}}^{-1}\boldsymbol{\mu}_{\text{new}} = \mathbf{A}^{-1}\mathbf{a}+\mathbf{B}^{-1}\mathbf{b}$$

which is the same answer with no algebra at all. That is the point of asking for both routes.

Verified at $\mathbf{a}=(1,-0.5)$, $\mathbf{b}=(-0.4,1.1)$: the log scaling constant is $-3.7048850193$ from Equation 6.76, and integrating the product over a grid gives $2.4603046081\times10^{-2}$ against the predicted $2.4603046081\times10^{-2}$ — a relative gap of $2.8\times10^{-15}$.

6.11 — Iterated expectations

Show $\mathbb{E}_X[x] = \mathbb{E}_Y\bigl[\mathbb{E}_X[x\mid y]\bigr]$.

$$\mathbb{E}_Y\bigl[\mathbb{E}_X[x\mid y]\bigr] = \sum_y p(y)\sum_x x\,p(x\mid y) = \sum_y\sum_x x\,p(x\mid y)p(y)$$

The product rule (Equation 6.22) turns $p(x\mid y)p(y)$ into $p(x,y)$:

$$= \sum_y\sum_x x\,p(x,y) = \sum_x x\sum_y p(x,y) = \sum_x x\,p(x) = \mathbb{E}_X[x]$$

using the sum rule to collapse the inner sum. Two rules, three lines.

Verified on a $3\times3$ joint: $\mathbb{E}[x] = 0.530000000000$ both ways, gap $1.1\times10^{-16}$. The conditional means were $(0.2353,\ 0.6744,\ 0.5000)$ — very different from each other and from $0.53$. The identity is about their weighted average, not about any one of them.

6.12 — Manipulation of Gaussian random variables

$\mathbf{x}\sim\mathcal{N}(\boldsymbol{\mu}_x,\boldsymbol{\Sigma}_x)$ with $\mathbf{x}\in\mathbb{R}^D$, and $\mathbf{y} = \mathbf{A}\mathbf{x}+\mathbf{b}+\mathbf{w}$ with $\mathbf{w}\sim\mathcal{N}(\mathbf{0},\mathbf{Q})$ independent.

a. The likelihood. Conditioning on $\mathbf{x}$ freezes $\mathbf{A}\mathbf{x}+\mathbf{b}$, leaving only the noise:

$$p(\mathbf{y}\mid\mathbf{x}) = \mathcal{N}\bigl(\mathbf{y} \mid \mathbf{A}\mathbf{x}+\mathbf{b},\ \mathbf{Q}\bigr)$$

b. The marginal $p(\mathbf{y})$. Rather than doing the integral, use the fact that $\mathbf{y}$ is a sum of two independent Gaussians — $\mathbf{A}\mathbf{x}+\mathbf{b}$ and $\mathbf{w}$. By Equations 6.50 and 6.51 the first has mean $\mathbf{A}\boldsymbol{\mu}_x+\mathbf{b}$ and covariance $\mathbf{A}\boldsymbol{\Sigma}_x\mathbf{A}^\top$; by Equation 6.78 the sum adds both:

$$\boldsymbol{\mu}_y = \mathbf{A}\boldsymbol{\mu}_x+\mathbf{b}, \qquad \boldsymbol{\Sigma}_y = \mathbf{A}\boldsymbol{\Sigma}_x\mathbf{A}^\top + \mathbf{Q}$$

c. Then $\mathbf{z} = \mathbf{C}\mathbf{y}+\mathbf{v}$. Identically: $p(\mathbf{z}\mid\mathbf{y}) = \mathcal{N}(\mathbf{C}\mathbf{y},\mathbf{R})$, and

$$\boldsymbol{\mu}_z = \mathbf{C}\boldsymbol{\mu}_y, \qquad \boldsymbol{\Sigma}_z = \mathbf{C}\boldsymbol{\Sigma}_y\mathbf{C}^\top + \mathbf{R}$$

d. The posterior $p(\mathbf{x}\mid\hat{\mathbf{y}})$. Following the book’s hint, build the joint. The cross-covariance is Equation 6.52:

$$\operatorname{Cov}[\mathbf{x},\mathbf{y}] = \boldsymbol{\Sigma}_x\mathbf{A}^\top$$

(the noise $\mathbf{w}$ is independent of $\mathbf{x}$, so it contributes nothing). Now apply Gaussian conditioning, Equations 6.66 and 6.67:

$$\boldsymbol{\mu}_{x\mid\hat{\mathbf{y}}} = \boldsymbol{\mu}_x + \boldsymbol{\Sigma}_x\mathbf{A}^\top\bigl(\mathbf{A}\boldsymbol{\Sigma}_x\mathbf{A}^\top+\mathbf{Q}\bigr)^{-1}\bigl(\hat{\mathbf{y}}-\mathbf{A}\boldsymbol{\mu}_x-\mathbf{b}\bigr)$$ $$\boldsymbol{\Sigma}_{x\mid\hat{\mathbf{y}}} = \boldsymbol{\Sigma}_x - \boldsymbol{\Sigma}_x\mathbf{A}^\top\bigl(\mathbf{A}\boldsymbol{\Sigma}_x\mathbf{A}^\top+\mathbf{Q}\bigr)^{-1}\mathbf{A}\boldsymbol{\Sigma}_x$$

This is Bayesian linear regression, and Chapter 9 is largely this formula.

Verified with $D=3$, $E=2$, $F=2$ against two million simulated draws:

quantity	worst gap against simulation
$\boldsymbol{\mu}_y$	$0.0060$
$\boldsymbol{\Sigma}_y$	$0.0143$
$\boldsymbol{\mu}_z$	$0.0027$
$\boldsymbol{\Sigma}_z$	$0.0056$
$\operatorname{Cov}[\mathbf{x},\mathbf{y}]$	$0.0180$

and the posterior trace fell $17.391171 \to 7.110717$, with $\boldsymbol{\Sigma}_x-\boldsymbol{\Sigma}_{x\mid\hat{\mathbf{y}}}$ having eigenvalues $[0,\ 3.5938,\ 6.6866]$ — positive semidefinite, with one direction unlearned because $\mathbf{y}$ has only two components against $\mathbf{x}$‘s three.

6.13 — The probability integral transform

Given continuous $x$ with cdf $F_x$, show $y = F_x(x)$ is uniform.

Take $u\in[0,1]$ and compute the cdf of $y$:

$$P(y\leqslant u) = P\bigl(F_x(x)\leqslant u\bigr) = P\bigl(x\leqslant F_x^{-1}(u)\bigr) = F_x\bigl(F_x^{-1}(u)\bigr) = u$$

The middle step needs $F_x$ to be strictly monotonic so the inverse exists and applying it preserves the inequality. A random variable whose cdf equals $u$ on $[0,1]$ is uniform on $[0,1]$, so $y\sim\mathcal{U}[0,1]$. $\blacksquare$

Verified on three distributions:

$X$ drawn from	mean of $F_x(X)$	variance	worst $\lvert F_{\text{emp}}-u\rvert$
$\mathcal{N}(0,1)$	$0.500063$	$0.083365$	$0.001340$
Exponential$(1.5)$	$0.499978$	$0.083265$	$0.002682$
Gamma$(2,1)$	$0.500376$	$0.083258$	$0.002350$

against a uniform’s $0.5$ and $1/12 = 0.083333$.

All thirteen, in one table

#	What it is	Key result
6.1	joint table	marginals by summing an axis; conditionals by dividing by a row or column total
6.2	the bimodal mixture	mean $(4.0, 0.8)$ exactly; mode at $(10, 2)$; $p(x_1)$ bimodal, $p(x_2)$ not
6.3	Bernoulli conjugacy	$\operatorname{Beta}(\alpha+s,\ \beta+N-s)$
6.4	two bags	$P(\text{bag 2}\mid\text{mango}) = \tfrac13$
6.5	time series	the Kalman filter, from Equations 6.66/6.67/6.78/6.88
6.6	proof	$\mathbb{V}[x] = \mathbb{E}[x^2]-\mathbb{E}[x]^2$, by linearity
6.7	proof	pairwise sum is twice the raw-score variance; $N^2$ terms for $N$ terms’ worth
6.8	Bernoulli, exp family	$\theta = \log\frac{\mu}{1-\mu}$, $A(\theta) = \log(1+e^\theta)$
6.9	Binomial, Beta, product	natural parameters add, which is conjugacy
6.10	product of Gaussians	two routes, same $\mathbf{C}$ and $\mathbf{c}$; precisions add
6.11	iterated expectations	product rule then sum rule, three lines
6.12	Gaussian linear model	Bayesian linear regression’s posterior
6.13	proof	$F_x(x)$ is uniform, needing only strict monotonicity

The complete check

chapter_6_exercises.py

"""Every exercise at the end of Chapter 6, checked numerically."""
 
import math
 
import numpy as np
 
np.set_printoptions(precision=6, suppress=True, linewidth=150)
 
print("########## ex_6_1")
# Rows are y1, y2, y3; columns x1..x5.
P = np.array([
    [0.01, 0.02, 0.03, 0.10, 0.10],
    [0.05, 0.10, 0.05, 0.07, 0.20],
    [0.10, 0.05, 0.03, 0.05, 0.04],
])
print("joint p(x, y):")
print(P)
print(f"  sums to {P.sum():.10f}  -- a valid joint")
px = P.sum(axis=0)
py = P.sum(axis=1)
print(f"a. p(x) = {px}   sums to {px.sum():.10f}")
print(f"   p(y) = {py}   sums to {py.sum():.10f}")
# b. p(x | Y = y1) is the y1 ROW divided by its total.
cond_x_y1 = P[0] / py[0]
# p(y | X = x3) is the x3 COLUMN divided by its total.
cond_y_x3 = P[:, 2] / px[2]
print(f"b. p(x | Y=y1) = {cond_x_y1}")
print(f"      row total p(y1) = {py[0]:.4f},  conditional sums to {cond_x_y1.sum():.10f}")
print(f"   p(y | X=x3) = {cond_y_x3}")
print(f"      column total p(x3) = {px[2]:.4f},  conditional sums to {cond_y_x3.sum():.10f}")
print("   note which total divides which: conditioning on Y fixes a ROW.")
 
print()
print("########## ex_6_2")
W = np.array([0.4, 0.6])
MU = [np.array([10.0, 2.0]), np.array([0.0, 0.0])]
SIG = [np.array([[1.0, 0.0], [0.0, 1.0]]),
       np.array([[8.4, 2.0], [2.0, 1.7]])]
print("a. marginals: a Gaussian's marginal is Gaussian (Eq 6.68), and a mixture")
print("   marginalises component by component, so each marginal is a 1-D mixture:")
for d, name in ((0, "x1"), (1, "x2")):
    print(f"   p({name}) = {W[0]} N({MU[0][d]:.1f}, {SIG[0][d,d]:.1f})"
          f" + {W[1]} N({MU[1][d]:.1f}, {SIG[1][d,d]:.1f})")
 
grid = np.linspace(-14, 20, 400_001)
 
 
def norm1(t, m, v):
    return np.exp(-0.5 * (t - m) ** 2 / v) / np.sqrt(2 * np.pi * v)
 
 
print()
print("b. per-marginal mean, mode and median:")
print(f"   {'dim':>4}  {'mean':>10}  {'mode':>10}  {'median':>10}  {'modes found':>12}")
for d, name in ((0, "x1"), (1, "x2")):
    dens = W[0] * norm1(grid, MU[0][d], SIG[0][d, d]) \
        + W[1] * norm1(grid, MU[1][d], SIG[1][d, d])
    mass = np.trapezoid(dens, grid)
    mean = np.trapezoid(grid * dens, grid) / mass
    cdf = np.concatenate([[0.0], np.cumsum((dens[1:] + dens[:-1]) / 2 * np.diff(grid))])
    cdf /= cdf[-1]
    med = float(grid[int(np.searchsorted(cdf, 0.5))])
    peaks = [i for i in range(1, len(dens) - 1)
             if dens[i] > dens[i - 1] and dens[i] > dens[i + 1]
             and dens[i] > 0.02 * dens.max()]
    mode = float(grid[int(np.argmax(dens))])
    print(f"   {name:>4}  {mean:>10.6f}  {mode:>10.6f}  {med:>10.6f}  {len(peaks):>12}")
    # The exact mean, by linearity.
    print(f"        exact mean by Eq 6.34: "
          f"{W[0]*MU[0][d] + W[1]*MU[1][d]:.6f}")
 
print()
print("c. the two-dimensional mean and mode:")
mean2 = W[0] * MU[0] + W[1] * MU[1]
print(f"   mean, exactly by linearity: {mean2}")
# The mode of the joint: search the density on a grid.
g1 = np.linspace(-9, 15, 1200)
g2 = np.linspace(-7, 9, 1200)
G1, G2 = np.meshgrid(g1, g2)
pts = np.stack([G1.ravel(), G2.ravel()], axis=1)
dens2 = np.zeros(pts.shape[0])
for w, m, S in zip(W, MU, SIG):
    d = pts - m
    Si = np.linalg.inv(S)
    dens2 += w * np.exp(-0.5 * np.einsum("ij,jk,ik->i", d, Si, d)) \
        / (2 * np.pi * np.sqrt(np.linalg.det(S)))
i = int(np.argmax(dens2))
print(f"   mode, by grid search: {pts[i]}   density {dens2[i]:.6f}")
# Compare the density at the two component centres.
for w, m, S in zip(W, MU, SIG):
    dd = np.zeros(1)
    for w2, m2, S2 in zip(W, MU, SIG):
        v = (m - m2)[None, :]
        dd += w2 * np.exp(-0.5 * v @ np.linalg.inv(S2) @ v.T).ravel() \
            / (2 * np.pi * np.sqrt(np.linalg.det(S2)))
    print(f"   density at component centre {m}: {float(dd[0]):.6f}")
print("   the mode is NOT the mean, and the mean sits between the two humps.")
 
print()
print("########## ex_6_3")
# Conjugate prior for a Bernoulli likelihood is the Beta (Table 6.2).
print("  conjugate prior: Beta(alpha, beta).  With N observations summing to s:")
print("  p(mu | x_1..x_N) = Beta(alpha + s, beta + N - s)")
rng = np.random.default_rng(3)
mu_true = 0.7
xs = (rng.random(50) < mu_true).astype(int)
a0, b0 = 2.0, 2.0
s = int(xs.sum())
N = xs.size
print(f"  simulated: N = {N}, successes s = {s}, prior Beta({a0}, {b0})")
print(f"  posterior Beta({a0 + s}, {b0 + N - s})   mean "
      f"{(a0+s)/(a0+b0+N):.6f}   (truth {mu_true})")
# Check against a grid.
gg = np.linspace(1e-6, 1 - 1e-6, 200_001)
lik = gg ** s * (1 - gg) ** (N - s)
pri = gg ** (a0 - 1) * (1 - gg) ** (b0 - 1)
post = lik * pri
post /= np.trapezoid(post, gg)
lc = math.lgamma(a0 + s + b0 + N - s) - math.lgamma(a0 + s) - math.lgamma(b0 + N - s)
closed = np.exp(lc + (a0 + s - 1) * np.log(gg) + (b0 + N - s - 1) * np.log1p(-gg))
print(f"  grid posterior vs closed form: worst gap {np.abs(post - closed).max():.2e}")
print("  the compiler either works or it does not, and one Beta absorbs every run.")
 
print()
print("########## ex_6_4")
# Two bags, a biased coin, Bayes.
p_h = 0.6
p_mango_b1 = 4 / 6
p_mango_b2 = 4 / 8
p_mango = p_h * p_mango_b1 + (1 - p_h) * p_mango_b2
post_b2 = (1 - p_h) * p_mango_b2 / p_mango
print(f"  P(heads) = {p_h}, so bag 1 with prob {p_h} and bag 2 with prob {1-p_h}")
print(f"  P(mango | bag 1) = 4/6 = {p_mango_b1:.6f}")
print(f"  P(mango | bag 2) = 4/8 = {p_mango_b2:.6f}")
print(f"  evidence P(mango) = {p_h}*{p_mango_b1:.4f} + {1-p_h}*{p_mango_b2:.4f}"
      f" = {p_mango:.6f}")
print(f"  P(bag 2 | mango) = {1-p_h}*{p_mango_b2:.4f} / {p_mango:.6f} = "
      f"{post_b2:.10f}  = 1/3")
print(f"  P(bag 1 | mango) = {1 - post_b2:.10f}  = 2/3   (they sum to "
      f"{post_b2 + (1-post_b2):.10f})")
# Simulate to confirm.
M = 4_000_000
heads = rng.random(M) < p_h
draw = rng.random(M)
is_mango = np.where(heads, draw < p_mango_b1, draw < p_mango_b2)
print(f"  simulated P(bag 2 | mango) = "
      f"{float((~heads & is_mango).sum() / is_mango.sum()):.6f}")
print("  note the prior favoured bag 1 (0.6) AND bag 1 is richer in mangos, so")
print("  seeing a mango pushes further toward bag 1: 0.4 prior -> 0.333 posterior.")
 
print()
print("########## ex_6_5_kalman")
# The time-series model. Everything stays Gaussian, so only means and
# covariances need propagating.
A = np.array([[0.9, 0.2], [-0.1, 0.95]])
Q = np.array([[0.30, 0.05], [0.05, 0.20]])
C = np.array([[1.0, 0.5]])
R = np.array([[0.25]])
mu_t = np.array([1.0, -0.5])
S_t = np.array([[0.40, 0.10], [0.10, 0.30]])
print("a. p(x_0, ..., x_T) is GAUSSIAN: x_{t+1} = A x_t + w is an affine map of a")
print("   Gaussian plus independent Gaussian noise, and Eq 6.78 and 6.88 say both")
print("   operations preserve Gaussianity. The joint of jointly-Gaussian variables")
print("   is Gaussian, so the whole trajectory is one big Gaussian.")
print()
# b1: the predictive.
mu_pred = A @ mu_t
S_pred = A @ S_t @ A.T + Q
print(f"b1. p(x_{{t+1}} | y_1..y_t) = N(A mu_t, A Sigma_t A^T + Q)")
print(f"    mean {mu_pred}")
print(f"    cov\n{S_pred}")
# b2: the joint of x_{t+1} and y_{t+1}.
mu_y = C @ mu_pred
S_yy = C @ S_pred @ C.T + R
S_xy = S_pred @ C.T
print(f"b2. the joint p(x_{{t+1}}, y_{{t+1}} | y_1..y_t) is Gaussian with")
print(f"    mean_x {mu_pred}   mean_y {mu_y}")
print(f"    Sigma_xx\n{S_pred}")
print(f"    Sigma_xy {S_xy.ravel()}   Sigma_yy {S_yy.ravel()}")
# b3: condition on the observation. This is the Kalman update.
y_hat = np.array([0.8])
K = S_xy @ np.linalg.inv(S_yy)
mu_post = mu_pred + (K @ (y_hat - mu_y)).ravel()
S_post = S_pred - K @ S_xy.T
print(f"b3. observing y = {y_hat[0]}, Eq 6.66 and 6.67 give")
print(f"    gain K = {K.ravel()}")
print(f"    posterior mean {mu_post}")
print(f"    posterior cov\n{S_post}")
print(f"    trace fell {np.trace(S_pred):.6f} -> {np.trace(S_post):.6f}, and")
print(f"    Sigma_pred - Sigma_post has eigenvalues "
      f"{np.linalg.eigvalsh(S_pred - S_post)}  -> PSD, as conditioning must be")
print("    that is exactly the Kalman filter: predict with 6.78/6.88, update with")
print("    6.66/6.67. No new mathematics is needed.")
 
print()
print("########## ex_6_6_and_6_7")
d = np.array([2.0, 4.0, 4.0, 4.0, 5.0, 5.0, 7.0, 9.0])
n = d.size
mu_d = d.mean()
lhs66 = float(np.mean((d - mu_d) ** 2))
rhs66 = float(np.mean(d ** 2) - mu_d ** 2)
print("6.6  E[(x-mu)^2] = E[x^2] - 2 mu E[x] + mu^2 = E[x^2] - mu^2")
print(f"     measured: {lhs66:.12f} against {rhs66:.12f}   gap {abs(lhs66-rhs66):.1e}")
lhs67 = float(np.sum((d[:, None] - d[None, :]) ** 2) / n ** 2)
print("6.7  expanding (x_i - x_j)^2 = x_i^2 - 2 x_i x_j + x_j^2 and summing over")
print("     both indices gives 2 N sum x_i^2 - 2 (sum x_i)^2, hence 2[E[x^2]-E[x]^2]")
print(f"     measured: {lhs67:.12f} against 2 * {rhs66:.12f} = {2*rhs66:.12f}")
print(f"     gap {abs(lhs67 - 2*rhs66):.1e}   -- N^2 = {n**2} terms for what N = {n} give")
 
print()
print("########## ex_6_8_and_6_9")
print("6.8  Bernoulli: p(x|mu) = exp[x log(mu/(1-mu)) + log(1-mu)], so")
print("     phi(x) = x,  theta = log(mu/(1-mu)),  A(theta) = log(1+e^theta), h(x) = 1")
for m in (0.25, 0.6):
    th = math.log(m / (1 - m))
    for xv in (0, 1):
        direct = m ** xv * (1 - m) ** (1 - xv)
        expfam = math.exp(xv * th - math.log1p(math.exp(th)))
        print(f"     mu={m}, x={xv}: direct {direct:.10f}  exp-family {expfam:.10f}"
              f"  gap {abs(direct-expfam):.1e}")
print()
print("6.9  Binomial(N, mu) at fixed N: phi(m) = m, theta = log(mu/(1-mu)),")
print("     A(theta) = N log(1+e^theta), h(m) = C(N, m)")
Nb, mub = 12, 0.35
thb = math.log(mub / (1 - mub))
worst = 0.0
for k in range(Nb + 1):
    direct = math.comb(Nb, k) * mub ** k * (1 - mub) ** (Nb - k)
    ef = math.comb(Nb, k) * math.exp(k * thb - Nb * math.log1p(math.exp(thb)))
    worst = max(worst, abs(direct - ef))
print(f"     worst gap over all {Nb+1} outcomes: {worst:.1e}")
print("     Beta(a,b): phi(mu) = [log mu, log(1-mu)], theta = [a-1, b-1],")
print("     A(theta) = log B(a,b), h(mu) = 1")
aa, bb = 3.0, 4.0
gm = np.linspace(1e-6, 1 - 1e-6, 20001)
lb = math.lgamma(aa) + math.lgamma(bb) - math.lgamma(aa + bb)
direct = np.exp(-lb + (aa - 1) * np.log(gm) + (bb - 1) * np.log1p(-gm))
phi_mu = np.stack([np.log(gm), np.log1p(-gm)], axis=1)
ef = np.exp(phi_mu @ np.array([aa - 1, bb - 1]) - lb)
print(f"     worst gap over the grid: {np.abs(direct - ef).max():.1e}")
print("     and the PRODUCT of the two: adding exponents means adding natural")
print("     parameters, so the product is exponential-family with theta summed --")
print("     which is exactly why the Beta is conjugate for the Binomial.")
 
print()
print("########## ex_6_10_product_of_gaussians")
a_, A_ = np.array([1.0, -0.5]), np.array([[1.2, 0.3], [0.3, 0.8]])
b_, B_ = np.array([-0.4, 1.1]), np.array([[0.6, -0.2], [-0.2, 1.5]])
C_ = np.linalg.inv(np.linalg.inv(A_) + np.linalg.inv(B_))
c_ = C_ @ (np.linalg.solve(A_, a_) + np.linalg.solve(B_, b_))
sign, ld = np.linalg.slogdet(A_ + B_)
diff = a_ - b_
log_scale = -np.log(2 * np.pi) - 0.5 * ld - 0.5 * diff @ np.linalg.solve(A_ + B_, diff)
print("a. completing the square: the exponents are two quadratics in x; adding them")
print("   gives a single quadratic whose Hessian is A^-1 + B^-1, so C = (A^-1+B^-1)^-1")
print("   and whose linear term gives c = C(A^-1 a + B^-1 b). What is left over does")
print("   not involve x, and is the scaling constant.")
print("b. in exponential-family form: a Gaussian's natural parameters are")
print("   [Sigma^-1 mu, -Sigma^-1/2]. Multiplying densities ADDS natural parameters,")
print("   so Sigma_new^-1 = A^-1 + B^-1 and Sigma_new^-1 mu_new = A^-1 a + B^-1 b.")
print("   Both routes give the same C and c, which is the point of the exercise.")
print(f"   C =\n{C_}")
print(f"   c = {c_}")
print(f"   log scaling constant = {log_scale:.10f}")
# Verify by integration.
gx = np.linspace(-8, 8, 1600)
X1, X2 = np.meshgrid(gx, gx)
pp = np.stack([X1.ravel(), X2.ravel()], axis=1)
 
 
def gauss(pts, m, S):
    d0 = pts - m
    Si = np.linalg.inv(S)
    return np.exp(-0.5 * np.einsum("ij,jk,ik->i", d0, Si, d0)) \
        / (2 * np.pi * np.sqrt(np.linalg.det(S)))
 
 
prod = gauss(pp, a_, A_) * gauss(pp, b_, B_)
dxg = gx[1] - gx[0]
mass = prod.sum() * dxg * dxg
print(f"   integral of the product: {mass:.10e}   predicted {np.exp(log_scale):.10e}")
print(f"   relative gap {abs(mass - np.exp(log_scale))/np.exp(log_scale):.2e}")
 
print()
print("########## ex_6_11_iterated_expectations")
# E_X[x] = E_Y[E_X[x|y]].  Check on a discrete joint.
J = np.array([[0.05, 0.10, 0.02],
              [0.15, 0.08, 0.20],
              [0.10, 0.20, 0.10]])
J = J / J.sum()
xv = np.array([-1.0, 0.5, 2.0])       # states of X, indexing columns
yv = np.array([0.0, 1.0, 2.0])        # states of Y, indexing rows
p_y = J.sum(axis=1)
E_x = float((J.sum(axis=0) * xv).sum())
E_x_given_y = (J / p_y[:, None]) @ xv
outer = float((p_y * E_x_given_y).sum())
print("  the proof: E_Y[E_X[x|y]] = sum_y p(y) sum_x x p(x|y) = sum_y sum_x x p(x,y)")
print("  = sum_x x sum_y p(x,y) = sum_x x p(x) = E_X[x]. The product rule is the")
print("  only step, and the sum rule collapses the inner sum.")
print(f"  E[x] directly              {E_x:.12f}")
print(f"  E_Y[E_X[x|y]]              {outer:.12f}")
print(f"  gap                        {abs(E_x - outer):.1e}")
print(f"  the inner conditional means are {E_x_given_y}, which differ a lot -- the")
print("  identity is about their weighted average, not about them individually.")
 
print()
print("########## ex_6_12_gaussian_manipulation")
D_, E_, F_ = 3, 2, 2
rngz = np.random.default_rng(12)
Araw = rngz.normal(size=(D_, D_))
Sx = Araw @ Araw.T + D_ * np.eye(D_)
mux = rngz.normal(size=D_)
Am = rngz.normal(size=(E_, D_))
bm = rngz.normal(size=E_)
Qm = np.diag(np.abs(rngz.normal(size=E_)) + 0.4)
print(f"a. p(y | x) = N(y | A x + b, Q)   -- shapes A {Am.shape}, Q {Qm.shape}")
mu_y = Am @ mux + bm
S_y = Am @ Sx @ Am.T + Qm
print(f"b. mu_y = A mu_x + b = {mu_y}")
print(f"   Sigma_y = A Sigma_x A^T + Q =\n{S_y}")
Cm = rngz.normal(size=(F_, E_))
Rm = np.diag(np.abs(rngz.normal(size=F_)) + 0.3)
mu_z = Cm @ mu_y
S_z = Cm @ S_y @ Cm.T + Rm
print(f"c. p(z | y) = N(z | C y, R);  mu_z = C mu_y = {mu_z}")
print(f"   Sigma_z = C Sigma_y C^T + R =\n{S_z}")
# d. the posterior p(x | y_hat), via the joint and Gaussian conditioning.
S_xy = Sx @ Am.T                       # Cov[x, y] = Sigma_x A^T
y_obs = mu_y + np.array([0.7, -0.4])
Kg = S_xy @ np.linalg.inv(S_y)
mu_post = mux + Kg @ (y_obs - mu_y)
S_post = Sx - Kg @ S_xy.T
print(f"d. Cov[x,y] = Sigma_x A^T, shape {S_xy.shape}")
print(f"   posterior mean {mu_post}")
print(f"   posterior cov\n{S_post}")
print(f"   trace {np.trace(Sx):.6f} -> {np.trace(S_post):.6f}")
print(f"   Sigma_x - Sigma_post eigenvalues {np.linalg.eigvalsh(Sx - S_post)}  -> PSD")
# Confirm b and c by simulation.
M2 = 2_000_000
xsim = rngz.multivariate_normal(mux, Sx, M2)
ysim = xsim @ Am.T + bm + rngz.multivariate_normal(np.zeros(E_), Qm, M2)
zsim = ysim @ Cm.T + rngz.multivariate_normal(np.zeros(F_), Rm, M2)
print(f"   simulated: mu_y gap {np.abs(ysim.mean(axis=0)-mu_y).max():.4f}, "
      f"Sigma_y gap {np.abs(np.cov(ysim.T, bias=True)-S_y).max():.4f}")
print(f"              mu_z gap {np.abs(zsim.mean(axis=0)-mu_z).max():.4f}, "
      f"Sigma_z gap {np.abs(np.cov(zsim.T, bias=True)-S_z).max():.4f}")
print(f"              Cov[x,y] gap "
      f"{np.abs(np.cov(np.hstack([xsim, ysim]).T, bias=True)[:D_, D_:] - S_xy).max():.4f}")
 
print()
print("########## ex_6_13_probability_integral_transform")
print("  proof: let Y = F_x(X) and take u in [0,1]. Then")
print("    P(Y <= u) = P(F_x(X) <= u) = P(X <= F_x^-1(u)) = F_x(F_x^-1(u)) = u,")
print("  using strict monotonicity to invert. A cdf equal to u on [0,1] IS the")
print("  uniform, so Y ~ U[0,1].")
rngp = np.random.default_rng(13)
for name, smp, cdf in (
        ("N(0,1)", rngp.standard_normal(2_000_000),
         lambda t: 0.5 * (1 + np.vectorize(math.erf)(t / math.sqrt(2)))),
        ("Exponential(1.5)", rngp.exponential(1 / 1.5, 2_000_000),
         lambda t: 1 - np.exp(-1.5 * t)),
        ("Gamma(2, 1)", rngp.gamma(2.0, 1.0, 2_000_000),
         lambda t: 1 - np.exp(-t) * (1 + t))):
    u = cdf(smp)
    srt = np.sort(u[:200_000])
    ks = float(np.abs(srt - np.linspace(0, 1, srt.size)).max())
    print(f"  {name:>18}  mean {u.mean():.6f}  var {u.var():.6f}  "
          f"worst |F_emp - u| {ks:.6f}")
print("  a uniform has mean 0.5 and variance 0.083333.")

output

########## ex_6_1
joint p(x, y):
[[0.01 0.02 0.03 0.1  0.1 ]
 [0.05 0.1  0.05 0.07 0.2 ]
 [0.1  0.05 0.03 0.05 0.04]]
  sums to 1.0000000000  -- a valid joint
a. p(x) = [0.16 0.17 0.11 0.22 0.34]   sums to 1.0000000000
   p(y) = [0.26 0.47 0.27]   sums to 1.0000000000
b. p(x | Y=y1) = [0.038462 0.076923 0.115385 0.384615 0.384615]
      row total p(y1) = 0.2600,  conditional sums to 1.0000000000
   p(y | X=x3) = [0.272727 0.454545 0.272727]
      column total p(x3) = 0.1100,  conditional sums to 1.0000000000
   note which total divides which: conditioning on Y fixes a ROW.
 
########## ex_6_2
a. marginals: a Gaussian's marginal is Gaussian (Eq 6.68), and a mixture
   marginalises component by component, so each marginal is a 1-D mixture:
   p(x1) = 0.4 N(10.0, 1.0) + 0.6 N(0.0, 8.4)
   p(x2) = 0.4 N(2.0, 1.0) + 0.6 N(0.0, 1.7)
 
b. per-marginal mean, mode and median:
    dim        mean        mode      median   modes found
     x1    4.000008    9.998390    2.803905             2
        exact mean by Eq 6.34: 4.000000
     x2    0.800000    1.330770    0.868710             1
        exact mean by Eq 6.34: 0.800000
 
c. the two-dimensional mean and mode:
   mean, exactly by linearity: [4.  0.8]
   mode, by grid search: [9.99583  1.994162]   density 0.063734
   density at component centre [10.  2.]: 0.063735
   density at component centre [0. 0.]: 0.029783
   the mode is NOT the mean, and the mean sits between the two humps.
 
########## ex_6_3
  conjugate prior: Beta(alpha, beta).  With N observations summing to s:
  p(mu | x_1..x_N) = Beta(alpha + s, beta + N - s)
  simulated: N = 50, successes s = 34, prior Beta(2.0, 2.0)
  posterior Beta(36.0, 18.0)   mean 0.666667   (truth 0.7)
  grid posterior vs closed form: worst gap 8.79e-14
  the compiler either works or it does not, and one Beta absorbs every run.
 
########## ex_6_4
  P(heads) = 0.6, so bag 1 with prob 0.6 and bag 2 with prob 0.4
  P(mango | bag 1) = 4/6 = 0.666667
  P(mango | bag 2) = 4/8 = 0.500000
  evidence P(mango) = 0.6*0.6667 + 0.4*0.5000 = 0.600000
  P(bag 2 | mango) = 0.4*0.5000 / 0.600000 = 0.3333333333  = 1/3
  P(bag 1 | mango) = 0.6666666667  = 2/3   (they sum to 1.0000000000)
  simulated P(bag 2 | mango) = 0.333517
  note the prior favoured bag 1 (0.6) AND bag 1 is richer in mangos, so
  seeing a mango pushes further toward bag 1: 0.4 prior -> 0.333 posterior.
 
########## ex_6_5_kalman
a. p(x_0, ..., x_T) is GAUSSIAN: x_{t+1} = A x_t + w is an affine map of a
   Gaussian plus independent Gaussian noise, and Eq 6.78 and 6.88 say both
   operations preserve Gaussianity. The joint of jointly-Gaussian variables
   is Gaussian, so the whole trajectory is one big Gaussian.
 
b1. p(x_{t+1} | y_1..y_t) = N(A mu_t, A Sigma_t A^T + Q)
    mean [ 0.8   -0.575]
    cov
[[0.672   0.1545 ]
 [0.1545  0.45575]]
b2. the joint p(x_{t+1}, y_{t+1} | y_1..y_t) is Gaussian with
    mean_x [ 0.8   -0.575]   mean_y [0.5125]
    Sigma_xx
[[0.672   0.1545 ]
 [0.1545  0.45575]]
    Sigma_xy [0.74925  0.382375]   Sigma_yy [1.190438]
b3. observing y = 0.8, Eq 6.66 and 6.67 give
    gain K = [0.62939  0.321205]
    posterior mean [ 0.98095  -0.482653]
    posterior cov
[[ 0.200429 -0.086163]
 [-0.086163  0.332929]]
    trace fell 1.127750 -> 0.533358, and
    Sigma_pred - Sigma_post has eigenvalues [0.       0.594392]  -> PSD, as conditioning must be
    that is exactly the Kalman filter: predict with 6.78/6.88, update with
    6.66/6.67. No new mathematics is needed.
 
########## ex_6_6_and_6_7
6.6  E[(x-mu)^2] = E[x^2] - 2 mu E[x] + mu^2 = E[x^2] - mu^2
     measured: 4.000000000000 against 4.000000000000   gap 0.0e+00
6.7  expanding (x_i - x_j)^2 = x_i^2 - 2 x_i x_j + x_j^2 and summing over
     both indices gives 2 N sum x_i^2 - 2 (sum x_i)^2, hence 2[E[x^2]-E[x]^2]
     measured: 8.000000000000 against 2 * 4.000000000000 = 8.000000000000
     gap 0.0e+00   -- N^2 = 64 terms for what N = 8 give
 
########## ex_6_8_and_6_9
6.8  Bernoulli: p(x|mu) = exp[x log(mu/(1-mu)) + log(1-mu)], so
     phi(x) = x,  theta = log(mu/(1-mu)),  A(theta) = log(1+e^theta), h(x) = 1
     mu=0.25, x=0: direct 0.7500000000  exp-family 0.7500000000  gap 0.0e+00
     mu=0.25, x=1: direct 0.2500000000  exp-family 0.2500000000  gap 5.6e-17
     mu=0.6, x=0: direct 0.4000000000  exp-family 0.4000000000  gap 0.0e+00
     mu=0.6, x=1: direct 0.6000000000  exp-family 0.6000000000  gap 0.0e+00
 
6.9  Binomial(N, mu) at fixed N: phi(m) = m, theta = log(mu/(1-mu)),
     A(theta) = N log(1+e^theta), h(m) = C(N, m)
     worst gap over all 13 outcomes: 1.7e-16
     Beta(a,b): phi(mu) = [log mu, log(1-mu)], theta = [a-1, b-1],
     A(theta) = log B(a,b), h(mu) = 1
     worst gap over the grid: 1.3e-15
     and the PRODUCT of the two: adding exponents means adding natural
     parameters, so the product is exponential-family with theta summed --
     which is exactly why the Beta is conjugate for the Binomial.
 
########## ex_6_10_product_of_gaussians
a. completing the square: the exponents are two quadratics in x; adding them
   gives a single quadratic whose Hessian is A^-1 + B^-1, so C = (A^-1+B^-1)^-1
   and whose linear term gives c = C(A^-1 a + B^-1 b). What is left over does
   not involve x, and is the scaling constant.
b. in exponential-family form: a Gaussian's natural parameters are
   [Sigma^-1 mu, -Sigma^-1/2]. Multiplying densities ADDS natural parameters,
   so Sigma_new^-1 = A^-1 + B^-1 and Sigma_new^-1 mu_new = A^-1 a + B^-1 b.
   Both routes give the same C and c, which is the point of the exercise.
   C =
[[0.376271 0.020339]
 [0.020339 0.482567]]
   c = [ 0.237288 -0.160533]
   log scaling constant = -3.7048850193
   integral of the product: 2.4603046081e-02   predicted 2.4603046081e-02
   relative gap 2.82e-15
 
########## ex_6_11_iterated_expectations
  the proof: E_Y[E_X[x|y]] = sum_y p(y) sum_x x p(x|y) = sum_y sum_x x p(x,y)
  = sum_x x sum_y p(x,y) = sum_x x p(x) = E_X[x]. The product rule is the
  only step, and the sum rule collapses the inner sum.
  E[x] directly              0.530000000000
  E_Y[E_X[x|y]]              0.530000000000
  gap                        1.1e-16
  the inner conditional means are [0.235294 0.674419 0.5     ], which differ a lot -- the
  identity is about their weighted average, not about them individually.
 
########## ex_6_12_gaussian_manipulation
a. p(y | x) = N(y | A x + b, Q)   -- shapes A (2, 3), Q (2, 2)
b. mu_y = A mu_x + b = [-2.419224  2.837356]
   Sigma_y = A Sigma_x A^T + Q =
[[ 17.57386  -13.710442]
 [-13.710442  19.444137]]
c. p(z | y) = N(z | C y, R);  mu_z = C mu_y = [ 0.877603 -0.221628]
   Sigma_z = C Sigma_y C^T + R =
[[3.621597 1.61282 ]
 [1.61282  4.471699]]
d. Cov[x,y] = Sigma_x A^T, shape (3, 2)
   posterior mean [0.732067 0.171968 0.379045]
   posterior cov
[[ 1.141729  0.108221 -2.077995]
 [ 0.108221  1.413374 -0.873665]
 [-2.077995 -0.873665  4.555614]]
   trace 17.391171 -> 7.110717
   Sigma_x - Sigma_post eigenvalues [0.       3.593818 6.686635]  -> PSD
   simulated: mu_y gap 0.0060, Sigma_y gap 0.0143
              mu_z gap 0.0027, Sigma_z gap 0.0056
              Cov[x,y] gap 0.0180
 
########## ex_6_13_probability_integral_transform
  proof: let Y = F_x(X) and take u in [0,1]. Then
    P(Y <= u) = P(F_x(X) <= u) = P(X <= F_x^-1(u)) = F_x(F_x^-1(u)) = u,
  using strict monotonicity to invert. A cdf equal to u on [0,1] IS the
  uniform, so Y ~ U[0,1].
              N(0,1)  mean 0.500063  var 0.083365  worst |F_emp - u| 0.001340
    Exponential(1.5)  mean 0.499978  var 0.083265  worst |F_emp - u| 0.002682
         Gamma(2, 1)  mean 0.500376  var 0.083258  worst |F_emp - u| 0.002350
  a uniform has mean 0.5 and variance 0.083333.

pch.quizTag Check your understanding

1. In 6.2, the mixture puts weight 0.4 on the component at (10, 2) and 0.6 on the one at (0, 0). Where is the mode?

   A At (0, 0), the heavier component B At (10, 2) — the lighter component, because its covariance is the identity while the other's is much broader, so its mass concentrates into a taller peak. Measured densities: 0.063735 against 0.029783 C At the mean, (4.0, 0.8) D There are two equal modes
   pch.quizShowAnswer

   B — At (10, 2) — the lighter component, because its covariance is the identity while the other's is much broader, so its mass concentrates into a taller peak. Measured densities: 0.063735 against 0.029783 — Weight is not height. And the mean at (4.0, 0.8) sits in the low-density gap between the clusters — a point the distribution almost never produces, which is why a mean is a poor summary of a multimodal density.

2. Exercise 6.5 asks for the form of the joint over a whole trajectory. What is the argument?

   A It must be computed component by component B It is Gaussian, because the dynamics apply only two operations — an affine map (Eq 6.88) and adding independent Gaussian noise (Eq 6.78) — and neither leaves the Gaussian family C It is Gaussian only if A is invertible D It is a mixture of Gaussians
   pch.quizShowAnswer

   B — It is Gaussian, because the dynamics apply only two operations — an affine map (Eq 6.88) and adding independent Gaussian noise (Eq 6.78) — and neither leaves the Gaussian family — That is why the exercise says you need not compute the joint explicitly: the closure properties settle the form. Parts b1 to b3 then produce the Kalman filter from Equations 6.66 and 6.67 with no new mathematics.

3. In 6.4 the prior on bag 2 is 0.4 and the posterior after seeing a mango is 1/3. Why did it go down?

   A Because the evidence P(mango) is less than one B Because bag 1 is the mango-richer bag — 2/3 against 1/2 — so a mango is evidence FOR bag 1. Prior and likelihood point the same way and reinforce each other C Because the coin is biased D It should have gone up
   pch.quizShowAnswer

   B — Because bag 1 is the mango-richer bag — 2/3 against 1/2 — so a mango is evidence FOR bag 1. Prior and likelihood point the same way and reinforce each other — Worth checking the direction of every Bayes update. Here P(mango) = 0.6 exactly, and the posterior 0.2/0.6 = 1/3 confirms it — simulated at 0.333517.

4. Exercise 6.9 asks you to show the product of a Beta and a Binomial is in the exponential family. What is the one-line reason?

   A Because both are bounded on [0,1] B Because both are exponentials of linear functions of the same sufficient statistics, so multiplying ADDS the exponents — which means adding the natural parameters, leaving the form unchanged. That is exactly conjugacy C Because the Beta is the conjugate prior D Because the product is normalisable
   pch.quizShowAnswer

   B — Because both are exponentials of linear functions of the same sufficient statistics, so multiplying ADDS the exponents — which means adding the natural parameters, leaving the form unchanged. That is exactly conjugacy — It also explains WHY Beta-Binomial conjugacy holds, rather than just verifying it: theta_prior + theta_likelihood gives Beta(alpha+m, beta+N-m), which is Example 6.11's answer arrived at structurally.

5. In both 6.5 and 6.12 the covariance difference has a zero eigenvalue. What does that mean?

   A The computation is numerically unstable B One direction in the state learned nothing, because the observation has fewer components than the state — a scalar y against a 2-D state in 6.5, a 2-D y against a 3-D x in 6.12 C The posterior is singular D The Kalman gain is wrong
   pch.quizShowAnswer

   B — One direction in the state learned nothing, because the observation has fewer components than the state — a scalar y against a 2-D state in 6.5, a 2-D y against a 3-D x in 6.12 — Equation 6.67 subtracts a positive-semidefinite term whose rank is at most the observation's dimension. A zero eigenvalue is the honest report that the unobserved direction is exactly as uncertain as it was before.

🧪 Try It Yourself

Exercise 1 – 6.1, marginals and conditionals

pch.exerciseTag pch.exerciseTitle

Exercise 2 – 6.4, the mango

pch.exerciseTag pch.exerciseTitle

Exercise 3 – 6.5, one Kalman step

pch.exerciseTag pch.exerciseTitle

Exercise 4 – 6.6 and 6.7, the two variance identities

pch.exerciseTag pch.exerciseTitle

Exercise 5 – 6.13, the probability integral transform

pch.exerciseTag pch.exerciseTitle

Recall card

• 6.1: marginals sum an axis, conditionals divide by a total. Conditioning on Y fixes a ROW and divides by p(y); conditioning on X fixes a COLUMN. Each conditional must sum to 1 — the check that catches a transpose.
• 6.2: the mixture’s mean is exact by linearity, (4.0, 0.8), but its mode is at (10, 2) — the LIGHTER component, because its covariance is tighter and so its peak is taller. Measured 0.063735 against 0.029783.
• 6.2: marginalising can lose a mode. p(x1) is bimodal with peaks near 0 and 10; p(x2) is unimodal with its peak at 1.3308, at neither component mean.
• 6.3: the Beta is conjugate for a Bernoulli, so the posterior is Beta(alpha + s, beta + N - s). The whole run history collapses to two numbers.
• 6.4: P(bag 2 | mango) = 1/3. The prior on bag 2 was 0.4 and it FELL, because bag 1 is the mango-richer bag, so prior and likelihood reinforce.
• 6.5 is the Kalman filter, and it needs no new mathematics: predict with Eq 6.78 and 6.88, update with Eq 6.66 and 6.67. Measured: trace 1.1278 down to 0.5334 after one scalar observation.
• A scalar observation can only inform one direction. The covariance reduction had eigenvalues [0, 0.594392] — the zero is the unobserved direction, and it is honest rather than an error.
• 6.6: V[x] = E[x^2] - E[x]^2 follows from linearity of expectation plus the fact that mu is a constant.
• 6.7: the pairwise sum over N^2 terms equals TWICE the raw-score variance. The cross term factorises into the square of a sum, which is what collapses N^2 into N.
• 6.8 and 6.9: Bernoulli, Binomial and Beta are all exponential-family. The Binomial’s A(theta) is the Bernoulli’s scaled by N; the Beta’s sufficient statistics are (log mu, log(1-mu)).
• 6.9’s real content: multiplying exponential-family densities ADDS natural parameters. That is why the Beta is conjugate for the Binomial — conjugacy is structural, not lucky.
• 6.10: two derivations, one answer. Completing the square and adding natural parameters both give C = (A-inverse + B-inverse)-inverse. Precisions add.
• 6.11: iterated expectations is the product rule then the sum rule, three lines. The conditional means were 0.235, 0.674, 0.500 — the identity is about their weighted average.
• 6.12 is Bayesian linear regression. Cov[x,y] = Sigma_x A-transpose, then Gaussian conditioning. Verified against two million simulated draws to within 0.018.
• 6.13: F_x(x) is uniform, and the proof needs only strict monotonicity — enough to invert F and preserve the inequality. Verified on a Gaussian, an exponential and a Gamma.

Next: Chapter 6 Formula Sheet — every definition, theorem and identity in the chapter on one page.