Past exam of the mathematical and Computational Biology course of the University of Cambridge / 2024 / Paper 2 / B4 / b / Solution 2026-10-20
View more
A single randomly chosen internal-site flip is a symmetric proposal: the reverse flip has the same probability . Therefore the Metropolis–Hastings algorithm accepts with
(1)
Comparing logarithms avoids overflow. On rejection the old state must be retained, including a repeated copy in the stored Markov chain; retaining only accepted states would sample a different distribution. This implementation also stores the actual initial state at row zero:
def sample_methylation(observed, nsteps=5000, seed=2024, score=methylation_score):
    rng=np.random.default_rng(seed)
    state=rng.integers(0,2,size=len(observed));state[[0,-1]]=1
    states=np.empty((nsteps,len(state)),dtype=int);values=np.empty(nsteps)
    states[0]=state;values[0]=score(state,observed)
    accepted=0
    for t in range(1,nsteps):
        candidate=state.copy();site=rng.integers(1,len(state)-1)
        candidate[site]=1-candidate[site]
        trial=score(candidate,observed)
        if np.log(rng.random())<min(0.,trial-values[t-1]):
            state=candidate;values[t]=trial;accepted+=1
        else:values[t]=values[t-1]
        states[t]=state
    return values,states,accepted/(nsteps-1)

lls, meths, acceptance = sample_methylation(obs, nsteps=5000, seed=2024)
import matplotlib.pyplot as plt
plt.plot(lls)
plt.xlabel("Stored state index")
plt.ylabel("Conditional-product log score")
plt.show()
For this random-number convention, the acceptance fraction is approximately . The trace climbs substantially during its first several thousand updates, so a small fixed burn-in such as 100 iterations would discard too little of the initial transient. Five thousand single-site proposals are only about fifty updates per site on average.
Figure 1.
Conditional-product sampling trace and retained-sample methylation profile, with a provisional burn-in of 3000 stored states
. Details.
Read the full article
Past exam of the mathematical and Computational Biology course of the University of Cambridge / 2024 / Paper 2 / B4 / c / Solution 2026-10-20
View more
For a binary state, its expected value is its probability of being one. Hence the Monte Carlo method estimate of the posterior methylation profile under the chosen target is
(1)
The averaging axis must be zero: average over sampled patterns while keeping the site index. A provisional choice for this particular trace is , retaining 2000 states:
def mean_methylation(meths, burnin):
    if not 0 <= burnin < len(meths):
        raise ValueError("burn-in must leave at least one state")
    return np.asarray(meths)[burnin:].mean(axis=0)

burnin = 3000
profile_estimate = mean_methylation(meths, burnin)
print(lls[burnin:].mean())
plt.plot(profile_estimate, label="retained-sample mean")
plt.scatter(np.arange(len(obs)), .2*obs, s=10, label="0.2 × observation")
plt.xlabel("Site index")
plt.ylabel("Methylation probability")
plt.legend()
plt.show()
The retained mean log score is , and the profile is shown below the trace in part (b). The fixed endpoints have estimate exactly one; adjacent sites are often smoothed into methylated or unmethylated stretches despite individual conflicting observations. The rescaled observations are only a visual comparison and are not another estimate of the posterior probability.
This is a finite-chain estimate for the conditional-product target defined in part (a). A burn-in chosen by inspection does not establish convergence, and the long excursions in this short Markov chain leave substantial uncertainty in individual-site estimates; several longer runs would be needed for precise inference. To sample the coherent joint alternative instead, pass score=joint_methylation_score and inspect that new trace before choosing its own burn-in.
Read the full article
Past exam of the systems biology course of the University of Cambridge / 2024 / Paper 3 / 2 / e / Solution 2026-10-20
View more
TODO: Complete the data-specific numerical results and diagnostics when the original Iris_data.csv is available.
Plot all 10,001 stored states, including the initial value, and indicate the provisional burn-in boundary. Also display the retained chain's autocorrelation function:
png("q2-trace.png", width = 1000, height = 700)
par(mfrow = c(2, 1), mar = c(4, 4.5, 2.5, 1))
plot(0:B, mu, type = "l", xlab = "Iteration", ylab = expression(mu),
     main = "Complete Metropolis trace")
abline(v = burn, col = "firebrick", lty = 2)
abline(h = m, col = "steelblue", lty = 3)
acf(post, main = "Autocorrelation after provisional burn-in")
dev.off()
# Check sensitivity to the chosen discard and to the retained time window.
for (discard in c(500L, 1000L, 2000L)) {
  kept <- mu[(discard + 2L):(B + 1L)]
  cat("Discard", discard, ": posterior mean", mean(kept), "\n")
}
cat("Means of retained halves:", mean(post[1:4500]),
    mean(post[4501:9000]), "\n")
A satisfactory trace plot should show an initial transient, if any, followed by repeated exploration of one stable range with no continuing drift. Repeated values are expected because rejections are stored. If is large, the prescribed proposal standard deviation one can be much wider than the posterior scale , giving low acceptance and conspicuous plateaus. A flat-looking trace is therefore not by itself evidence of good exploration.
Compare the retained histogram with the exact posterior density, inspect autocorrelation, compare separate retained windows, and assess sensitivity to the discarded count. Agreement supports the adequacy of this run; neither 10,000 iterations nor a chosen burn-in of 1,000 proves convergence. The effective sample size of a Markov chain can be much smaller than the number of stored states.
The actual trace cannot be diagnosed without Iris_data.csv. The code produces the required trace when the measurements are available; the comments here specify what to examine rather than asserting an unseen chain has converged.
Read the full article