序列熵

Question

有人能帮我纠正我的代码吗？ 我有3个序列

 seq1='AATC'
 seq2='GCCT'
 seq3='ATCA'

我要计算三个序列的熵，它的公式是：

^{pr2}$

式中，p峎ia是i列中字母a的概率，C峎ia是a列中字母a的数目[：][i]

示例：从以上3个序列中：第一个列：a'取值C，A和C，这里C_ia'=3，如果A='C'C_ia=2，那么p_ia=2/3

如果a='a'C_ia=1，则p_ia=1/3

Entropy (A[:][0])=-(2/3 log (2/3)+1/3 log(1/3))

其他列也是一样。 我试过了： 此函数计算给定字符串中每个字母的出现（频率），l是字符串的长度，表示l=∑ C_ia'.

def occ(source):
    occ={}
    l=0
    for e in source:
        l +=1
        if e not in occ:
            occ[e]=0
        occ[e]+=1
    return (l,occ)

示例：

>>> source='aaazzzedc'
>>> occ(source)
(9, {'z': 3, 'c': 1, 'a': 3, 'd': 1, 'e': 1})

那么，这就是熵的函数：

def entropy(occ,l):
    elist=[]
    for v in occ.values():
        c=v/l
        elist.append(-c *math.log(c, 2))
    return sum(elist)

示例：

>>> (l,h)=occ(source)
>>> entropy(h,l)
2.113283334294875

好吧。 这是一个简单的例子，但对我来说，我想把它应用到序列的列上。在

>>> source=[occ(t) for t in zip(seq1,seq2,seq3)]
>>> source
[(3, {'A': 2, 'G': 1}), (3, {'C': 1, 'A': 1, 'T': 1}), (3, {'C': 2, 'T': 1}), (3, {'C': 1, 'A': 1, 'T': 1})]

如何计算这些序列的熵？ 第一列给出：

(3, {'A': 2, 'G': 1})
2/3 * log(2/3)+ 1/3 * log (1/3)

如何写出所有列的熵？在