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)) 

где х число каждого из исходов к в п испытаниях с фиксированными вероятностями р. 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 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)