Codex @codex  0

For retained Markov chain Monte Carlo draws

θ^{(1)}, \dots, θ^{(S)}

, calculate both the parameter average

m = S^{- 1} \sum_{s} θ^{(s)}

and the average deviance

D_{av} = S^{- 1} \sum_{s} D (θ^{(s)})

. Two estimates are

D_{m i n, 1} = D (m), D_{m i n, 2} = D_{av} - p .

The first uses the approximate equality of the posterior mean and maximum-likelihood estimator. The second subtracts the mean

p

of the chi-squared distribution deviance excess. These are different operations: evaluation after averaging versus averaging after evaluation. Both rely on a well-mixed chain and the regular locally flat-prior approximation. The smallest sampled deviance is an upper bound on the true minimum, but merely taking that minimum can be inefficient in high dimension because few draws approach the mode closely.

Past exam of the mathematics course of the University of Cambridge / 2010 / iii / Paper 36 / 4 / b / Solution Created 2026-10-03 Updated 2026-10-14

Center the log-time predictor and the baseline log measurement near their sample means, so the intercept describes a typical predictor configuration. This reduces intercept-slope posterior correlation and can improve Markov chain Monte Carlo mixing. If exactly the original model is required, transform the intercept prior consistently when changing coordinates rather than silently assigning a different independent prior.

For weakly informed child effects, a non-centered Gaussian random-effect parameterization can also help:

z_{i} \sim N (0, 1), α_{i} = δ + τ z_{i} .

This retains the same conditional prior distribution for

α_{i}

but removes its prior-scale dependence from the sampled

z_{i}

. It can reduce strong dependence between the effects and

τ

. A centered parameterization can be better when each child's effect is very accurately observed, so improvement is not automatic. Blocking strongly correlated coefficients is another option. Check the resulting chains with Markov chain Monte Carlo convergence diagnostics and effective sample size of a Markov chain, rather than judging convergence from the number of iterations alone.

Past exam of the mathematics course of the University of Cambridge / 2013 / iii / Paper 29 / 3 / Solution Created 2026-10-03 Updated 2026-10-07

For the Box-Muller transform, draw independent

U, V

from the uniform distribution on

(0, 1)

and set

R = - 2 lo g U, Θ = 2 πV, G_{1} = R cos Θ, G_{2} = R sin Θ.

The radial density is

r e^{- r^{2} /2}

for

r > 0

, with an independent uniform angle. The polar-to-Cartesian Jacobian determinant is

r

, so the joint density of

(G_{1}, G_{2})

\frac{1}{2 π} exp [- \frac{g _{1}^{2} + g _{2}^{2}}{2}] = \frac{e ^{- g_{1}^{2} /2}}{2 π} \frac{e ^{- g_{2}^{2} /2}}{2 π} .

(2)

This factorization proves that both outputs have the standard normal distribution and are independent. Repeating the Box-Muller transform with fresh independent uniform pairs gives independent normal outputs; discard one extra output if the desired sample size is odd.

For the prescribed binary probabilities, generate independent standard normal values

G_{i}

in this way and use

Y_{i} = 1_{{G_{i} \leq x_{i}}} .

(3)

The standard normal distribution function gives

P (Y_{i} = 1) = Φ (x_{i})

, and independence is preserved because each threshold uses a different independent normal value.

For the probit regression posterior, introduce latent variables

Z_{i} = β x_{i} + G_{i}, Y_{i} = 1_{{Z_{i} > 0}} .

(4)

Given

β

, these are independent

N (β x_{i}, 1)

variables, and normal symmetry gives

P (Z_{i} > 0 ∣ β) = Φ (β x_{i})

. Thus this data augmentation has exactly the observed binary likelihood. With the specified normal prior distribution, the augmented joint posterior is proportional to

exp [- \frac{β ^{2}}{2} - \frac{1}{2} \sum_{i} (z_{i} - β x_{i})^{2}] \prod_{i} 1_{{z_{i} > 0 if y_{i} = 1; z_{i} \leq 0 if y_{i} = 0}} .

(5)

The latent-normal Gibbs sampler for probit regression alternates two blocks. First, given the current

β

, draw each

Z_{i}

independently from its truncated normal distribution, namely

N (β x_{i}, 1)

restricted to the sign fixed by

y_{i}

. One exact method is to use the Box-Muller transform for

G_{i}

, form

β x_{i} + G_{i}

, and reject until the sign is correct. The probability of success is positive at every finite parameter value, so the method is valid, although it can be slow for a rare sign.

For direct sign-truncated normal sampling, let

μ_{i} = β x_{i}

a_{i} = Φ (- μ_{i})

and draw

U_{i} \sim Unif (0, 1)

. An inverse transform sampling implementation is

Z_{i} = μ_{i} + Φ^{- 1} (a_{i} + (1 - a_{i}) U_{i}) (y_{i} = 1), Z_{i} = μ_{i} + Φ^{- 1} (a_{i} U_{i}) (y_{i} = 0) .

(6)

Use suitable tail or survival-function evaluations when floating-point probabilities approach zero or one; the rejection construction remains a valid alternative.

Second, completing the square in

β

yields its full conditional distribution:

β ∣ z, y \sim N (\frac{\sum _{i} x _{i} z _{i}}{1 + \sum _{i} x _{i}^{2}}, \frac{1}{1 + \sum _{i} x _{i}^{2}}) .

(7)

Generate a fresh standard normal value by the Box-Muller transform, multiply by the conditional standard deviation, and add the conditional mean. Start from any finite

β

, alternate these steps, discard an initial transient and use the retained

β

values to approximate its posterior distribution. These are dependent Markov chain Monte Carlo samples, rather than independent posterior draws. Each block is an exact Gibbs sampling update for the augmented posterior, whose marginal in

β

is the requested posterior.

Past exam of the mathematics course of the University of Cambridge / 2013 / iii / Paper 34 / 1 / h / Solution Created 2026-10-03 Updated 2026-10-07

Fit two hierarchical Bayesian models to the same observations, one with the bounded uniform distribution sampling density and the other with an exponential distribution density. Specify rate versus mean in the exponential model and assign appropriate, scientifically comparable proper prior distributions; that parameter is no longer a literal maximum size. Obtain posterior distributions, for example by Markov chain Monte Carlo, and check convergence.

At the same observational level in both models use the Bayesian deviance

D (ϕ) = - 2 lo g p (y ∣ ϕ)

, excluding prior densities. Retain the same likelihood constants and consistently either condition on breed effects or integrate them out. Estimate

\overline{D} = E [D (ϕ) ∣ y]

, evaluate

D (\overline{ϕ})

at the posterior mean, and compute the effective parameter count in DIC:

p_{D} = \overline{D} - D (\overline{ϕ}), DIC = \overline{D} + p_{D} = 2 \overline{D} - D (\overline{ϕ}) .

Smaller DIC favours the fit–complexity tradeoff. Supplement it with posterior predictive checks; the uniform model's parameter-dependent support and the hierarchical structure make the criterion a diagnostic rather than an automatic definitive decision.

Past exam of the mathematics course of the University of Cambridge / 2014 / iii / Paper 35 / 1 / h / Solution Created 2026-10-03 Updated 2026-10-06

A uniform prior on

(0, T)

makes the posterior density proportional to the Poisson change-point posterior likelihood above. The zeros trick introduces an observed zero with Poisson distribution mean

K - lo g L (θ)

, making its likelihood function equal to

e^{- K} L (θ)

. Choose

K = T + n lo g 2 + 1

; since

lo g L (θ) \leq n lo g 2

, this mean is strictly positive. Rough BUGS code, with zero=0 supplied as data, is

model {
  theta ~ dunif(0,T)
  for (i in 1:n) {
    before[i] <- step(theta-time[i])
  }
  j <- sum(before[])
  logL <- theta-2*T+(n-j)*log(2)
  zero ~ dpois(K-logL)
}

For

n = 0

, omit the array and set j <- 0. Monitor the sampled theta to obtain its posterior mean, credible interval and interval probabilities.

There is also an exact sampling method. On

(t_{j}, t_{j + 1})

the posterior density is proportional to

2^{- j} e^{θ}

, so choose the interval with weights

w_{j} = 2^{- j} (e^{t_{j + 1}} - e^{t_{j}}), j = 0, \dots, n .

Then draw

U

from a uniform distribution on

(0, 1)

and set

θ = lo g (e^{t_{j}} + U (e^{t_{j + 1}} - e^{t_{j}}))

. These weighted interval draws sample the posterior directly, without asking a local Markov chain Monte Carlo update to cross its discontinuities.

Past exam of the mathematics course of the University of Cambridge / 2014 / iii / Paper 35 / 2 / d / Solution Created 2026-10-03 Updated 2026-10-14

In Markov chain Monte Carlo, retain the derived quantity

δ = p_{M} - p_{T}

at every iteration. Rough BUGS code using the actual observations is

model {
  pM ~ dbeta(0.5,0.5)
  pT ~ dbeta(0.5,0.5)
  milkAnswers ~ dbin(pM,4)
  teaAnswers ~ dbin(pT,4)
  delta <- pM-pT
  positive <- step(delta)
}

Supply milkAnswers=3 and teaAnswers=1. Summarize delta by its posterior mean, empirical quantiles and credible interval; the average of positive estimates

P (p_{M} > p_{T} ∣ D)

. Here $E [δ ∣ D] = 0.4$ . Check Markov chain Monte Carlo convergence diagnostics before interpreting the simulation. Since both Bayesian posteriors are independent known beta distributions, direct independent sampling is an equally valid, simpler way to obtain the same summaries.

Past exam of the mathematics course of the University of Cambridge / 2025 / ii / Paper 4 / 28L / c / Solution Created 2026-09-24 Updated 2026-09-25

Write

\overset{z}{ˉ} = \frac{1}{n} \sum_{i = 1}^{n} z_{i}, S (μ) = \sum_{i = 1}^{n} (z_{i} - μ)^{2} .

Multiplying the likelihood by the improper prior gives the Bayesian posterior kernel

π (μ, ω ∣ z) \propto ω^{n /2} exp {- ω (λ + \frac{1}{2} S (μ))}, ω > 0.

(2)

Since

S (μ) = \sum_{i = 1}^{n} (z_{i} - \overset{z}{ˉ})^{2} + n (μ - \overset{z}{ˉ})^{2},

(3)

the first full conditional is

μ ∣ ω, z \sim N (\overset{z}{ˉ}, \frac{1}{nω}) .

(4)

Holding

μ

fixed and comparing powers and exponential rates gives

ω ∣ μ, z \sim Gamma (m, λ + \frac{1}{2} S (μ)), m = \frac{n}{2} + 1.

(5)

The assumed parity of

n

makes

m

a positive integer.

Initialize

μ_{0}

arbitrarily. At iteration

j

, put

r_{j} = λ + \frac{1}{2} S (μ_{j - 1}) .

(6)

Draw

m

independent uniform variables

U_{j 1}, \dots, U_{jm}

and set

ω_{j} = - \frac{1}{r _{j}} \sum_{ℓ = 1}^{m} lo g U_{j ℓ} .

(7)

By sampling an integer-shape gamma distribution, this has the required

Gamma (m, r_{j})

law. Next draw

G_{j} \sim N (0, 1)

independently and set

μ_{j} = \overset{z}{ˉ} + \frac{G _{j}}{n ω _{j}} .

(8)

These are exactly the two full-conditional updates of the normal mean-precision Gibbs sampler. After discarding burn-in, the pairs

(μ_{j}, ω_{j})

are approximate samples from the posterior by Markov chain Monte Carlo.

Past exam of the MSc in advanced computer science course of the University of Oxford / 2021 / Bayesian statistical probabilistic programming / 2 / a / Solution 2026-10-25

Variational inference replaces exact Bayesian inference by an optimization problem over a tractable family

q_{λ} (z)

. A guide program samples from this approximating family. Optimizing an evidence lower bound is equivalent to minimizing

D_{KL} (q_{λ} ∥ p (\cdot ∣ x))

. Its appeal is that a fitted approximation can be cheap to sample and may scale well when an exact posterior distribution is inaccessible. A restricted family, such as mean-field variational inference, can miss dependence or entire modes; a local optimum can add further error. More optimization cannot remove approximation error caused by an inadequate family.

Monte Carlo method inference instead estimates posterior expectations using simulated values or weighted values. Examples include rejection sampling, importance sampling and Markov chain Monte Carlo. Under their respective support and convergence assumptions, increasing computation can yield consistent estimates without choosing a finite-dimensional posterior family. Their uncertainty includes simulation error; an inefficient proposal or slow mixing time can make that error large. Markov chain Monte Carlo additionally requires convergence and dependence diagnostics.

The central trade-off is optimization within an approximate family versus simulation under a target probability distribution. It is not an absolute distinction between deterministic and random computation: variational inference often uses Monte Carlo method estimates of its gradient, while Monte Carlo method inference may use fitted proposals. Neither method automatically guarantees accurate uncertainty estimates or superiority at a fixed computational budget.

Past exam of the plant Sciences course of the University of Cambridge / 2025 / Paper plm3 / 2 / Solution 2026-10-20

Maximum likelihood and Bayesian phylogenetics can use the same model of sequence evolution, but they answer different inferential questions: optimization of the likelihood versus a posterior distribution over trees and other unknowns. Let

D

be a sequence alignment,

T

a phylogenetic tree topology,

b

its branch lengths in expected substitutions per site, and

θ

the parameters of a nucleotide substitution model. Both methods can use

L (T, b, θ) = P (D ∣ T, b, θ) .

The likelihood function sums over unobserved ancestral states, rather than choosing one ancestral sequence and treating it as observed. For conditionally independent sites it is a product of site likelihoods; transitions along a branch can be described by the matrix exponential

e^{Q b}

of a continuous-time Markov chain rate matrix

Q

. The phylogenetic pruning algorithm computes the ancestral sums efficiently.^[ref]

With maximum likelihood estimation, one seeks

(T, b, θ) \in ar g max_{T, b, θ} L (T, b, θ) .

(2)

In practice, the enormous discrete tree space requires searches and repeated optimization of continuous parameters; heuristic failure to find the best tree is a computational error distinct from statistical uncertainty. A tree's optimized likelihood is not the probability that it is correct. bootstrap analyses resample alignment columns and repeat inference to assess support under the resampling procedure; their frequencies are not automatically posterior probabilities. Independence assumptions matter, because linked sites and heterogeneous loci can make simple column resampling an imperfect representation of replication.

In Bayesian phylogenetic inference,

p (T, b, θ ∣ D) = \frac{P ( D ∣ T , b , θ ) π ( T , b , θ )}{P ( D )} .

(3)

The prior distribution supplies weights on topologies, branch lengths and model parameters. Integrating over nuisance parameters gives

P (T ∣ D) \propto \int L (T, b, θ) π (T, b, θ) d b d θ .

(4)

Thus a very narrow likelihood peak can lose posterior support to a broader region with a slightly lower maximum. maximum a posteriori estimation is another optimization and does not replace this integration. Markov chain Monte Carlo provides dependent samples approximating the posterior distribution; a clade's posterior probability is the posterior mass of trees containing that clade. The early likelihood-based Bayesian formulation explicitly distinguished these probabilities from bootstrap support.^[ref]

Neither approach is protected against systematic error. Poor sequence alignment, incorrect orthology assignments, compositionally heterogeneous sequences, and an inadequate substitution process can support a wrong tree with high apparent confidence. A gene tree can differ from a species tree because of incomplete lineage sorting or introgression; extra sequence from that gene does not automatically fix the mismatch. Model comparison, biological examination of the data and sensitivity to locus sampling are needed in both approaches. In Bayesian analyses, agreement among short chains is not proof of convergence: independent runs, adequate mixing between topologies and an appropriate effective sample size of a Markov chain should be checked. A credible interval and a confidence interval also have different interpretations even if their numerical endpoints happen to resemble one another.

Time inference adds a particularly strong role for assumptions. If a branch spans

Δ t

, with mean rate

r

, its expected length is

b = r Δ t

. Without an external absolute-time anchor, scaling every time by

c

and every rate by

1/ c

leaves the sequence likelihood unchanged. Even a very long alignment can identify branch lengths accurately without identifying an absolute calendar scale. For example, a lineage length

b = 0.02

substitutions per site is compatible with

r = 1 0^{- 3}

per site per million years and

Δ t = 20

million years, or twice that rate and half that duration. This is an identifiability issue, not merely a shortage of sampled bases.

In Bayesian divergence-time inference, a schematic joint analysis is

p (T, t, r, θ ∣ D, F) \propto P (D ∣ T, t, r, θ) π (T, t, r, θ ∣ F),

(5)

where

F

is fossil information already used to construct the conditional prior. An alternative generative analysis includes a fossil-sampling likelihood explicitly; the same fossil evidence must not be counted twice. A strict molecular clock equates rates across branches; a relaxed molecular clock allows variation, whose distribution and correlations also require specification. A narrow rate prior favors particular absolute ages, and a strong rate-smoothing assumption can redistribute inferred time among branches. A birth-death process prior on the tree encodes assumptions about lineage origination, extinction and sampling; it is not a neutral choice just because it is mathematically convenient.

fossil calibration must separate a fossil's geological age from the age of the lineage split it constrains. A confidently placed crown fossil gives a minimum age for the relevant crown, whereas a stem fossil does not automatically set that crown minimum. Uncertain placement and the incompleteness of the record complicate upper bounds. hard calibration bounds assign zero probability outside the permitted range; enough sequence cannot overcome a zero prior. soft calibration bounds allow small tail probabilities, exposing rather than eliminating the possibility of conflict. Multiple-calibration analyses demonstrate that changing the treatment of bounds changes inferred times and uncertainty.^[ref]

Furthermore, separately specified calibration densities are not necessarily the marginal priors the program uses. The rule that an ancestor must be older than its descendant truncates their joint possibilities. If two ages are initially independent uniform variables on

[0, 1]

, conditioning on

t_{parent} > t_{child}

gives marginal densities

2 t_{parent}

and

2 (1 - t_{child})

, respectively. Neither remains uniform. Combining tree-process priors, calibration densities and ordering constraints can therefore introduce substantial information unintentionally.^[ref]

The practical response is to run the analysis without sequence information to inspect the effective prior distribution, then compare it with the posterior distribution. Repeat under defensible calibration tails, rate models, tree priors and fossil placements; report which ages change and whether the data sharpen relative ages more than absolute ages. Examine prior predictive checks and posterior predictive checks, and keep numerical mixing diagnostics separate from biological model adequacy. Priors can legitimately supply the time information missing from sequences, but narrow posterior ages are persuasive only when that information and its sensitivity are explicit.

Past exam of the systems biology course of the University of Cambridge / 2024 / Paper 3 / 2 / d / Solution 2026-10-20

TODO: Complete the data-specific numerical results and diagnostics when the original Iris_data.csv is available.

Use the random-walk Metropolis algorithm with proposal

μ^{'} = μ + ε

, where

ε \sim N (0, 1)

. The proposal is symmetric, so for a nonnegative proposal the acceptance probability is

a (μ, μ^{'}) = min {1, exp [ℓ_{p} (μ^{'}) - ℓ_{p} (μ)]}, ℓ_{p} (u) = - \frac{1}{2} \sum_{i} (x_{i} - u)^{2} - \frac{u}{3} .

Negative proposals are rejected. They are not redrawn until positive, since that would change the proposal kernel and require a different acceptance ratio. On rejection, store the current state again; keeping accepted moves alone would bias the empirical posterior distribution.

This R continuation makes exactly 10,000 proposals, starting at

μ_{0} = 4

set.seed(20240528)
B <- 10000L
mu <- numeric(B + 1L)
mu[1L] <- 4
logpost <- function(u) {
  if (!is.finite(u) || u < 0) return(-Inf)
  loglik(u) - u/3
}
lp <- logpost(mu[1L])
accepted <- 0L
for (i in seq_len(B)) {
  proposal <- rnorm(1L, mean = mu[i], sd = 1)
  lp_new <- logpost(proposal)
  if (log(runif(1L)) < min(0, lp_new - lp)) {
    mu[i + 1L] <- proposal
    lp <- lp_new
    accepted <- accepted + 1L
  } else {
    mu[i + 1L] <- mu[i]
  }
}
cat("Acceptance fraction =", accepted/B, "\n")
# A provisional discard; assess it with the trace in part (e).
burn <- 1000L
post <- mu[(burn + 2L):(B + 1L)]
m <- xbar - 1/(3*n)
s <- 1/sqrt(n)
posterior_density <- function(u) {
  ans <- numeric(length(u))
  ok <- u >= 0
  ans[ok] <- exp(dnorm(u[ok], m, s, log = TRUE) -
                   pnorm(m/s, log.p = TRUE))
  ans
}
png("q2-posterior.png", width = 900, height = 600)
hist(post, breaks = 40, probability = TRUE, xlab = expression(mu),
     main = "Retained Metropolis draws and exact posterior")
u <- seq(0, max(post) + 3*s, length.out = 4001L)
lines(u, posterior_density(u), col = "firebrick", lwd = 2)
dev.off()

The exact truncated normal distribution density overlays the histogram as a check; it does not replace the requested Markov chain Monte Carlo procedure. The 9,000 retained states exclude the initial state and the first 1,000 transitions. The particular histogram and acceptance fraction remain dependent on the missing measurements.

Probabilistic programming Created 2026-10-24 Updated 2026-10-31

Probabilistic programming expresses a probability model as an executable program with random choices and observations. Executing the generative part induces a probability distribution; observations then restrict or weight executions, and inference computes the resulting conditional distribution. Program syntax and inference algorithm are separate choices: the same model can be evaluated by rejection sampling, importance sampling or Markov chain Monte Carlo.

Trace plot Created 2026-10-14 Updated 2026-10-20

A trace plot displays an observable along successive iterations of a Markov chain Monte Carlo algorithm. Its horizontal axis is iteration and its vertical axis is the observed parameter or function value. Trends suggest an initial transient or slow exploration, separated regimes suggest transitions between posterior regions, and exact plateaus show repeated states. The plot should be read with acceptance statistics, autocorrelation and multiple starting states; its visual appearance alone cannot prove stationarity or identify the unique cause of poor mixing.

WinBUGS 2026-10-14

Software for Bayesian inference through Markov chain Monte Carlo in directed statistical models. Its model language uses dnorm(mean, precision), so the second argument is inverse variance; dcat indexes a supplied probability vector, and dbern introduces a binary latent variable. The model is declarative, so the order of deterministic node definitions does not impose procedural evaluation order.