CODEWITHSUNDEEP
numpy
user51966
user51966

Reputation: 1057

Numpy: Loglikelihood of Multivariate Normal Distribution

I would like to calculate the loglikelihood of multivariate normal distribution.

Data:

data = np.random.multivariate_normal(mean=[2,5], cov=[[1, 0], [0, 10]], size=1000)

Likelihood (I followed Wikipedia):

def likelihood(mean, cov): 
    # mean = np.array([y1, y2]), cov = np.array([[c1, 0], [0, c2]])

    loglikelihood = -0.5*(  np.log(np.linalg.det(cov))  + (data - mean).transpose() * np.linalg.inv(cov) * (data - mean) + 2 * np.log(2 * np.pi)  )

    loglikelihoodsum = loglikelihood.sum()

    return loglikelihoodsum

It returns following error:

> likelihood(mean, cov)
---------------------------------------------------------------------------
ValueError                         Traceback (most recent call last)
<ipython-input-54-3b20c2eefea4> in <module>()
----> 1 likelihood(mean, cov)

<ipython-input-53-8a2a7219131c> in likelihood(mean, cov)
      2     # param = mean[y1, y2], cov = [[c1, 0], [0, c2]]
      3 
----> 4     loglikelihood = -0.5*(  np.log(np.linalg.det(cov))  + (data - mean).transpose() * np.linalg.inv(cov) * (data - mean) + 2 * np.log(2 * np.pi)  )
      5 
      6     loglikelihoodsum = loglikelihood.sum()

ValueError: operands could not be broadcast together with shapes (2,1000) (2,2)

How can I fix it?

Upvotes: 0

Views: 13061

Answers (3)

np_king
np_king

Reputation: 33

You can use scipy.stats.multivariate_normal.logpdf

Upvotes: 1

Till Hoffmann
Till Hoffmann

Reputation: 9887

Multiplication with the * operator in numpy refers to elementwise multiplication. You want to compute the inner product instead using np.einsum:

mean = np.random.normal(0, 1, 2)
cov = np.random.normal(0, 1, (2, 2))
data = np.random.normal(0, 1, (1000, 2))

residuals = data - mean
loglikelihood = -0.5 * (
    np.log(np.linalg.det(cov)) 
    + np.einsum('...j,jk,...k', residuals, np.linalg.inv(cov), residuals) 
    + len(mean) * np.log(2 * np.pi)
)
np.sum(loglikelihood)

Upvotes: 0

user51966
user51966

Reputation: 1057

I find an answer:

def likelihood(mean, cov): # Wikipedia
    def calc_loglikelihood(residuals):
        return -0.5 * (np.log(np.linalg.det(cov)) + residuals.T.dot(np.linalg.inv(cov)).dot(residuals) + 2 * np.log(2 * np.pi))

    # mean = np.array([y1, y2]), cov = np.array([[c1, 0], [0, c2]])
    residuals = (data - mean)

    loglikelihood = np.apply_along_axis(calc_loglikelihood, 1, residuals)
    loglikelihoodsum = loglikelihood.sum()

    return loglikelihoodsum

Upvotes: 2

Related Questions