= Solution
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
$$
\boxed{\alpha=\min\{1,\exp(s(m';o)-s(m;o))\}}.
$$
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, \b[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.
\Image[/past-exam-of-the-mathematical-and-computational-biology-course-of-the-university-of-cambridge/2024/paper-2-methylation.png]
{title=Conditional-product sampling trace and retained-sample methylation profile, with a provisional burn-in of 3000 stored states}
Back to article page
