Write for the true DNA methylation status and for its observation. The observation likelihood function is a binary symmetric channel, soFor the loop suggested in the practical, set and add when , or 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}")The conditional-product log score at the observed pattern is . 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 would instead be .
If the neighbour probabilities are intended literally as full conditionals of a Markov random field, a coherent alternative has edge weights and for , with both endpoints fixed. Given two equal neighbours, the matching centre has conditional odds ; mixed neighbours give equal weights. Consequently its unnormalized posterior probability has log scoreFor this alternative, 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()) Articles by others on the same topic
There are currently no matching articles.