A single randomly chosen internal-site flip is a symmetric proposal: the reverse flip has the same probability . Therefore the Metropolis–Hastings algorithm accepts withComparing 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.
Articles by others on the same topic
There are currently no matching articles.