2015-09-11 3 views
2

У меня есть следующая функция в python, которую я не могу понять, как выразить в векторизованной форме. Для меня cov - это массивный массив формы (2,2), mu - средний вектор с формой (2) и xtp имеет форму (~ 50000,2). Я знаю, что scipy предоставляет scipy.stats.multivariate_normal.pdf, но я пытаюсь научиться писать эффективный векторизованный код. ПожалуйстаКак записать этот код в векторном формате?

def mvnpdf(xtp, mu, cov): 
    temp = np.zeros(xtp.shape[0]) 
    i = 0 
    length = xtp.shape[0] 
    const = 1/(((2* np.pi)**(len(mu)/2)) * (np.linalg.det(cov)**(1/2))) 
    inv = np.linalg.inv(cov) 
    while i < length: 
     x = xtp[i]-mu 
     exponent = (-1/2) * (x.dot(inv).dot(x)) 
     temp[i] = (const * np.exp(exponent)) 
     i+=1 
    return temp 

ответ

4

Единственная сложная часть векторизации является то, что двойной .dot. Оставим это:

x = xtp - mu # move this out of the loop 
ddot = [i.dot(inv).dot(i) for i in x] 
temp = const * np.exp(-0.5 * ddot) 

Положите это в свой код и посмотрите, создает ли оно то же самое.

Существует несколько способов «векторизации» dot. Тот, который я хотел бы попробовать первым, это einsum. В моих тестах это эквивалентно:

ddot = np.einsum('ij,jk,ik->i',x,inv,x) 

Я предлагаю попробовать, чтобы увидеть, работает ли оно и ускоряет работу. И играйте с этими вычислениями с меньшими массивами (а не ~ 50000) в интерактивной оболочке.

Я проверяю вещи с

In [225]: x 
Out[225]: 
array([[ 0., 2.], 
     [ 1., 3.], 
     ... 
     [ 7., 9.], 
     [ 8., 10.], 
     [ 9., 11.]]) 
In [226]: inv 
Out[226]: 
array([[ 1., 0.], 
     [ 0., 1.]]) 

Поскольку это обучение упражнение, которое я оставлю детали к вам.

С (2,2), расчеты один cov можно быстрее сделать это явно, а не с функциями det и inv. Но это length итерация, которая является потребителем времени.

+0

Большое спасибо. einsum выглядит потрясающе. Полученный код намного быстрее. Есть ли у вас хорошие ресурсы для изучения использования enisum? – enitihas

+1

@enitihas: относительно немного ресурсов, разбросанных по сети, но вот вопрос SO с некоторыми ответами и ссылками: http://stackoverflow.com/questions/26089893/understanding-einsum-numpy –

Смежные вопросы