What is my problem with the implementation of multivariate_gauss pdf?

user_12721568 0
Viewed 47

I use python to calculate the multivariate_gauss distribution, but I don't know what's wrong. The code is here

# calculate multi-d gaussian pdf
def mul_gauss(x, mu, sigma) -> float:
    d = len(x[0])
    front = 1 / math.sqrt(((2 * math.pi) ** d) * np.linalg.det(sigma))
    tmp = (np.array(x) - np.array(mu))
    tmp_T = np.transpose(tmp)
    back = -0.5 * (np.matmul(np.matmul(tmp, np.linalg.inv(sigma)), tmp_T))[0][0]
    return front * math.exp(back)

I compared the result with scipy.stats.multivariate_normal(x,mu,sigma)

x = [[2,2]]
mu = [[4,4]]
sigma = [[3,0],[0,3]]
ret_1 = mul_gauss(x, mu, sigma)
ret_2 = scipy.stats.multivariate_normal(x[0], mu[0], sigma).pdf(x[0])
print('ret_1=',ret1)
print('ret_2=',ret2)

and output is ret_1=0.013984262505331654 ret_2=0.03978873577297383

Could anyone help me?

gaussian math numpy python
1 Answers

In line 5 of your main you call the .pdf() on the object instead as a method. Here is a fix:

# calculate multi-d gaussian pdf
import math

import numpy as np
from scipy import stats


def mul_gauss(x, mu, sigma) -> float:
    d = x[0].shape[0]
    coeff = 1/np.sqrt((2 * math.pi) ** d * np.linalg.det(sigma))
    tmp = x - mu
    exponent = -0.5 * (np.matmul(np.matmul(tmp, np.linalg.inv(sigma)), tmp.T))[0][0]
    return coeff * math.exp(exponent)


x = np.array([[2,2]])
mu = np.array([[4,4]])
sigma = np.array([[3,0],[0,3]])
ret_1 = mul_gauss(x, mu, sigma)
ret_2 = stats.multivariate_normal.pdf(x[0], mu[0], sigma)
print('ret_1=',ret_1)
print('ret_2=',ret_2)

Output:

ret_1= 0.013984262505331654
ret_2= 0.013984262505331658

Cheers.

user_45960780
Related