= Solution
Write $m_i\in\{0,1\}$ for the true <DNA methylation> status and $o_i$ for its observation. The observation <likelihood function> is a <binary symmetric channel>, so
$$
\log p(o\mid m)=\sum_{i=0}^{99}\begin{cases}\log(0.8),&o_i=m_i,\\\log(0.2),&o_i\ne m_i.\end{cases}
$$
For the loop suggested in the practical, set $p_i=(0.1,0.5,0.9)_{m_{i-1}+m_{i+1}}$ and add $\log p_i$ when $m_i=1$, or $\log(1-p_i)$ otherwise. The following independent implementation combines these terms without multiplying many small numbers:
``
import numpy as np
def methylation_score(pattern, observed):
"""Conditional-product score used by the practical's natural loop filling."""
pattern=np.asarray(pattern,dtype=int);observed=np.asarray(observed,dtype=int)
if pattern[0]!=1 or pattern[-1]!=1:return -np.inf
neighbour_sum=pattern[:-2]+pattern[2:]
p=np.array([.1,.5,.9])[neighbour_sum]
prior_terms=np.where(pattern[1:-1]==1,p,1-p)
observation_terms=np.where(pattern==observed,.8,.2)
return float(np.log(prior_terms).sum()+np.log(observation_terms).sum())
log_likelihood = methylation_score
print(f"{log_likelihood(obs, obs):.2f}")
``
\b[The conditional-product log score at the observed pattern is $-102.14$.] The endpoint observations are included once each, while the endpoint statuses are fixed rather than assigned an additional prior probability.
There is a modeling qualification: multiplying the specified neighbour <conditional probabilities> gives a <pseudolikelihood>, not generally the <joint probability distribution> having those <conditional probabilities>. Thus the function above defines the practical's conditional-product target; the pure observation <log-likelihood> at $m=o$ would instead be $100\log(0.8)=-22.3144$.
If the neighbour probabilities are intended literally as full conditionals of a <Markov random field>, a coherent alternative has edge weights $T_{aa}=3/4$ and $T_{ab}=1/4$ for $a\ne b$, with both endpoints fixed. Given two equal neighbours, the matching centre has conditional odds $(3/4)^2:(1/4)^2=9:1$; mixed neighbours give equal weights. Consequently its unnormalized <posterior probability> has log score
$$
\ell_J(m;o)=\sum_{i=0}^{98}\log T_{m_i m_{i+1}}+\log p(o\mid m).
$$
For this alternative, $\ell_J(o;o)=-94.7394$ before the constant associated with conditioning the endpoints. It can be implemented as follows, but its numerical results must not be mixed with the conditional-product results in the remaining parts:
``
def joint_methylation_score(pattern,observed):
"""Coherent nearest-neighbour model with the specified full conditionals."""
pattern=np.asarray(pattern,dtype=int);observed=np.asarray(observed,dtype=int)
if pattern[0]!=1 or pattern[-1]!=1:return -np.inf
edges=np.where(pattern[1:]==pattern[:-1],.75,.25)
emissions=np.where(pattern==observed,.8,.2)
return float(np.log(edges).sum()+np.log(emissions).sum())
``
Back to article page
