Numpy: Get Rid Of A Loop By Broadcasting

Question

I am trying to implement a Maximize Expectations Algorithm for a Gaussian mixture model in python.

I have the following line to calculate the Gaussian probability p of my X data, given the average mu and the covariance sigma of the Gaussian distribution:

for i in range(len(X[0])):  
   p[i] = scipy.stats.multivariate_normal.pdf(X[:,i],mu,sigma)

I wanted to know if I could somehow get rid of the for loop to get something like

p[:] = scipy.stats.multivariate_normal.pdf(X[:,:]??)

I did some research on broadcasting and thought about using a function numpy.einsum, but I can’t figure out how this will work in this case.

Answer 1

Flatten, use the call pdfand change the shape back -

from scipy import stats

out = stats.multivariate_normal.pdf(X.ravel(),mu,sigma).reshape(-1,len(X[0])).T