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 исходов в n испытаниях с фиксированной вероятностью p. во всех случаях n равно четырем четырем.

Я создал четыре функции для расчета этого распределения.

# 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' с первым компонентом каждого подсписка в amin_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)