def likelihood(p, data):
 k = sum(data)
 N = len(data)
 # Compute Binomial likelihood
 l = scipy.special.comb(N, k) * p**k * (1-p)**(N-k)
 # Normalize the likelihood to sum to unity
 return l/sum(l)