Python - 計算大數據集上的多項式概率密度函數？

Question

我原本打算使用MATLAB來解決這個問題，但內置函數的侷限性不符合我的目標。 NumPy中出現相同的限制。

我有兩個製表符分隔的文件。首先是一個文件表示的氨基酸殘基，頻率和計數蛋白質結構的一個內部數據庫，即

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

第二個文件由氨基酸四胞胎和它們發生的次數，即，

的

ASTP 1 

請注意，有> 8000這樣的四聯組。

根據每個氨基酸出現的背景頻率和四聯體計數，我的目標是計算每個四聯體的多項式概率密度函數，並隨後將其用作最大似然計算中的期望值。

的多項分佈如下：

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

其中x是每個k成果的與固定概率p n次試驗的數目。在我的計算中，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])) 

我haven'te真正得到如何整合其他功能呢。我認爲我目前的嵌套列表安排是次優的。

我希望將字符串'ASTP'中的每個字母與amino_acids中每個子列表的第一個組件進行比較。如果存在匹配，我希望將適當的數值傳遞給使用索引的函數。

他們是更好的方法嗎？我可以將每個氨基酸和四聯體的適當數字附加到循環中的臨時數據結構中，將其傳遞給函數並將其清除以用於下一次迭代？

感謝，S :-)

Answer 1

這可能是切到你原來的問題，但我強烈建議不要因溢出計算階乘明確。相反，利用factorial(n) = gamma(n+1)這一事實，使用伽瑪函數的對數，並使用加法代替乘法，減法而不是除法。 scipy.special包含一個名爲gammaln的函數，它給出伽馬函數的對數。

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) 

如果你不想安裝SciPy的只是爲gammaln緣故，這裏是純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)