Python - 如何在大型数据集上计算多项式概率密度函数？

Question

Python - 如何在大型数据集上计算多项式概率密度函数？

9

我最初打算使用MATLAB解决这个问题，但内置函数有一些限制，不符合我的目标。NumPy中也存在同样的限制。

我有两个制表符分隔的文件。第一个文件显示了一个蛋白质结构内部数据库的氨基酸残基、频率和计数。

A    0.25    1
S    0.25    1
T    0.25    1
P    0.25    1

第二个文件由氨基酸四元组及其出现次数组成，即：

ASTP    1

请注意，有超过8,000个这样的四重组。

基于每种氨基酸的背景出现频率和四重体的计数，我旨在计算每个四重体的多项式概率密度函数，并随后将其用作最大似然计算中的期望值。

多项式分布如下：

f(x|n, p) = n!/(x1!*x2!*...*xk!)*((p1^x1)*(p2^x2)*...*(pk^xk))

其中x是在n次试验中，每个k个结果的数量。在我的计算中，n在所有情况下均为4。

我创建了四个函数来计算此分布。

# functions for multinomial distribution


def expected_quadruplets(x, y):
    expected = x*y
    return expected

# calculates the probabilities of occurence raised to the number of occurrences

def prod_prob(p1, a, p2, b, p3, c, p4, d):
    prob_prod = (pow(p1, a))*(pow(p2, b))*(pow(p3, c))*(pow(p4, d))
    return prob_prod 


# factorial() and multinomial_coefficient() work in tandem to calculate C, the multinomial coefficient

def factorial(n):
    if n <= 1:
        return 1
    return n*factorial(n-1)


def multinomial_coefficient(a, b, c, d):
    n = 24.0
    multi_coeff =  (n/(factorial(a) * factorial(b) * factorial(c) * factorial(d)))
    return multi_coeff

问题在于如何最好地组织数据，以便以最高效的方式处理计算，以一种我可以阅读的方式（你们写的某些代码很神秘 :-))，并且不会创建溢出或运行时错误。

到目前为止，我的数据表示为嵌套列表。

amino_acids = [['A', '0.25', '1'], ['S', '0.25', '1'], ['T', '0.25', '1'], ['P', '0.25', '1']]

quadruplets = [['ASTP', '1']]

我最初打算在嵌套的for循环中调用这些函数，但这会导致运行时错误或溢出错误。我知道可以重置递归限制，但我更愿意以更优雅的方式解决这个问题。

我有以下代码：

for i in quadruplets:
    quad = i[0].split(' ')
    for j in amino_acids:
        for k in quadruplets:
            for v in k:
                if j[0] == v:
                    multinomial_coefficient(int(j[2]), int(j[2]), int(j[2]), int(j[2]))

我还没有涉及如何整合其他功能。我认为我的当前嵌套列表排列是次优的。

我希望将字符串“ASTP”中的每个字母与氨基酸中每个子列表的第一个组件进行比较。如果存在匹配，我希望使用索引将适当的数字值传递给函数。

有更好的方法吗？我能否在循环内将每个氨基酸和四联体的适当数字附加到临时数据结构中，将其传递给函数，并清除下一次迭代的内容？

谢谢，S :-)

- Darren J. Fitzpatrick

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Tamás · Accepted Answer

这可能与您最初的问题无关，但我强烈建议不要显式计算阶乘以避免溢出。相反，利用 factorial(n) = gamma(n+1) 这个事实，使用 gamma 函数的对数，并使用加法而不是乘法，减法而不是除法。 scipy.special 包含一个名为 gammaln 的函数，它可以给出 gamma 函数的对数。

from itertools import izip
from numpy import array, log, exp
from scipy.special import gammaln

def log_factorial(x):
    """Returns the logarithm of x!
    Also accepts lists and NumPy arrays in place of x."""
    return gammaln(array(x)+1)

def multinomial(xs, ps):
    n = sum(xs)
    xs, ps = array(xs), array(ps)
    result = log_factorial(n) - sum(log_factorial(xs)) + sum(xs * log(ps))
    return exp(result)

如果你不想仅仅为了使用 gammaln 而安装 SciPy，那么这里有一个纯 Python 的替代方案（当然它比 SciPy 中的版本速度慢且没有向量化）:

def gammaln(n):
    """Logarithm of Euler's gamma function for discrete values."""
    if n < 1:
        return float('inf')
    if n < 3:
        return 0.0
    c = [76.18009172947146, -86.50532032941677, \
         24.01409824083091, -1.231739572450155, \
         0.001208650973866179, -0.5395239384953 * 0.00001]
    x, y = float(n), float(n)
    tm = x + 5.5
    tm -= (x + 0.5) * log(tm)
    se = 1.0000000000000190015
    for j in range(6):
        y += 1.0
        se += c[j] / y
    return -tm + log(2.5066282746310005 * se / x)

另一个简单的技巧是使用一个由氨基酸本身作为索引的dict，用于amino_acids。在给定原始的amino_acids结构的情况下，您可以这样做：

amino_acid_dict = dict((amino_acid[0], amino_acid) for amino_acid in amino_acids)
print amino_acid_dict
{"A": ["A", 0.25, 1], "S": ["S", 0.25, 1], "T": ["T", 0.25, 1], "P": ["P", 0.25, 1]}

您可以更轻松地按残基查找频率或计数：

freq_A = amino_acid_dict["A"][1]
count_A = amino_acid_dict["A"][2]

这将为您节省主循环中的一些时间：

for quadruplet in quadruplets:
    probs = [amino_acid_dict[aa][1] for aa in quadruplet]
    counts = [amino_acid_dict[aa][2] for aa in quadruplet]
    print quadruplet, multinomial(counts, probs)