Python - tính toán các hàm mật độ xác suất đa thức trên tập dữ liệu lớn?

Question

Tôi ban đầu dự định sử dụng MATLAB để giải quyết vấn đề này nhưng chức năng tích hợp sẵn có những hạn chế không phù hợp với mục tiêu của tôi. Cùng một giới hạn xảy ra trong NumPy.

Tôi có hai tệp được phân tách bằng tab. Đầu tiên là một tập tin hiển thị amin dư lượng axit, tần số và đếm cho một cơ sở dữ liệu nội bộ của cấu trúc protein, tức là

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

Các tập tin thứ hai bao gồm bộ tứ của các axit amin và số lần chúng xảy ra, tức là

ASTP 1 

Lưu ý, có> 8.000 tứ phân như vậy.

Dựa trên tần suất xuất hiện của từng axit amin và số lượng tứ phân, tôi tính toán hàm mật độ xác suất đa thức cho mỗi tứ phân và sau đó sử dụng nó làm giá trị kỳ vọng trong phép tính khả năng tối đa.

Sự phân bố đa thức như sau:

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

trong đó x là số của mỗi kết quả k trong các thử nghiệm n với xác suất cố định p. n là 4 bốn trong mọi trường hợp trong tính toán của tôi.

Tôi đã tạo bốn hàm để tính toán phân phối này.

# 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 

Vấn đề là cách tốt nhất để cấu trúc dữ liệu để giải quyết việc tính toán một cách hiệu quả nhất, theo cách mà tôi có thể đọc (các bạn viết một số mã :-) khó hiểu) và điều đó sẽ không tạo ra một tràn hoặc lỗi thời gian chạy.

Cho đến nay dữ liệu của tôi được thể hiện dưới dạng danh sách lồng nhau.

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

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

Tôi dự định ban đầu gọi các hàm này trong vòng lặp lồng nhau nhưng điều này dẫn đến lỗi thời gian chạy hoặc lỗi tràn. Tôi biết rằng tôi có thể thiết lập lại giới hạn đệ quy nhưng tôi thà làm điều này một cách tao nhã hơn.

tôi đã có sau:

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

tôi haven'te thực sự nhận được để làm thế nào để kết hợp các chức năng khác được nêu ra. Tôi nghĩ rằng sắp xếp danh sách lồng nhau hiện tại của tôi là phụ tối ưu.

Tôi muốn so sánh từng chữ cái trong chuỗi 'ASTP' với thành phần đầu tiên của mỗi danh sách phụ trong amino_acids. Trong trường hợp một trận đấu tồn tại, tôi muốn chuyển các giá trị số thích hợp cho các hàm bằng cách sử dụng các chỉ mục.

Đây có phải là cách tốt hơn không? Tôi có thể nối thêm các số thích hợp cho mỗi axit amin và tứ phân đoạn vào một cấu trúc dữ liệu tạm thời trong vòng lặp hay không, chuyển nó tới các hàm và xóa nó cho lần lặp tiếp theo?

Cảm ơn, S :-)

Answer 1

Đây có thể là tiếp tuyến với câu hỏi ban đầu của bạn, nhưng tôi khuyên bạn không nên tính giai thừa một cách rõ ràng do tràn.Thay vào đó, hãy sử dụng thực tế là factorial(n) = gamma(n+1), sử dụng lôgarit của hàm gamma và sử dụng các phép cộng thay vì phép nhân, phép trừ thay vì các phần. scipy.special chứa một hàm có tên gammaln, cung cấp cho bạn logarit của hàm 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) 

Nếu bạn không muốn cài đặt scipy chỉ vì lợi ích của gammaln, đây là một sự thay thế bằng Python tinh khiết (tất nhiên nó là chậm hơn và nó không phải là vector hóa như một trong 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) 

Một mẹo nhỏ khác là sử dụng dict cho amino_acids, được lập chỉ mục bởi chính phần dư đó. Với gốc amino_acids cấu trúc của bạn, bạn có thể làm điều này:

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]} 

Bạn có thể sau đó nhìn lên các tần số hoặc số lượng bằng cách cặn dễ dàng hơn:

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

này giúp bạn tiết kiệm một số thời gian trong vòng lặp chính:

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)