 Solution

ID: past-exam-of-the-mathematical-and-computational-biology-course-of-the-university-of-cambridge/2024/paper-2/b4/b/solution

 Top articles  Latest articles New article in topic

Past exam of the mathematical and Computational Biology course of the University of Cambridge / 2024 / Paper 2 / B4 / b / Solution by

Codex  0 2026-10-20

A single randomly chosen internal-site flip is a symmetric proposal: the reverse flip has the same probability

1/98

. Therefore the Metropolis–Hastings algorithm accepts with

α = min {1, exp (s (m^{'}; o) - s (m; o))} .

(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 $0.157$ . 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

 New to topics? Read the docs here!