I was writing some code that needed to compute the probability density function (PDF) value for a multidimensional Gaussian vector. The regular one-dimensional Gaussian function PDF is the bell-shaped curve. For example if mean = u = 0.0, and standard deviation = sd = 1.0, and x = 2.5 then pdf(x) is the value of the curve at x = 2.5, which is about 0.01753. The pdf() value isn’t a probability, but larger pdf() values mean more likely.

In many scenarios, because pdf() values are often very small, using the log of the pdf() is more convenient mathematically. For the example above, logpdf(x) = log(0.01753) = -4.0438.

The multivariate case extends this idea. For example, suppose dimension = k = 3, mean = u = [1.0, 3.0, 2.0], covariance = the 3×3 identity matrix, and x = [1.5, 1.5, 1.5]. Each value is a vector with 3 elements, and instead of a single standard deviation value, a 3×3 covariance matrix is used. Because x = [1.5, 1.5, 1.5] is fairly close to the mean of u = [1.0, 3.0, 2.0], the pdf() value for x would be smaller than the pdf() value for x2 = [7.0, 9.5, 8.0] because x2 is farther away from u.

To compute logpdf() I used the built-in scipy.stats.multivariate_normal.logpdf() function from the SciPy library. But I wondered how easy it would be to remove the SciPy dependency and compute a logpdf() function using just NumPy — sort of from scratch.

I found the equation for the pdf() for a multivariate Gaussian distribution on Wikipedia (after wading through a lot of noise). The equation uses matrix transpose, inverse, determinant, and multiplication, but NumPy has functions for all of those. I coded up a program-defined logpdf() function using just NumPy functions. Good fun.

There are always tradeoffs between implementing a function from scratch vs. using an off-the-shelf library function.

Three from-scratch origami style paper dresses. More interesting than off-the-shelf clothing but probably less practical.

Demo code:

# multivariate_normal_pdf_demo.py

import numpy as np
import scipy.stats as sps

def my_logpdf(x, u, covar):
  k = len(x)  # dimension
  a = np.transpose(x - u)
  b = np.linalg.inv(covar)
  c = x - u
  d = np.matmul(a, b)
  e = np.matmul(d, c)
  numer = np.exp(-0.5 * e)
  f = (2 * np.pi)**k
  g = np.linalg.det(covar)
  denom = np.sqrt(f * g)
  pdf = numer / denom
  return np.log(pdf)

print("\nBegin multivariate Gaussian log-pdf demo ")

u = np.array([1.0, 3.0, 2.0], dtype=np.float32)
v = np.array([1.0, 1.0, 1.0], dtype=np.float32)
covar = np.diag(v)
x = np.array([1.5, 1.5, 1.5], dtype=np.float32)

print("\nu = ", end=""); print(u)
print("\ncovar = ")
print(covar)
print("\nx = ", end=""); print(x)

logp1 = sps.multivariate_normal.logpdf(x, u, covar)
print("\nlogpdf(x) using scipy logpdf() = %0.4f " % logp1)

logp2 = my_logpdf(x, u, covar)
print("\nlogpdf(x) using my_logpdf()    = %0.4f " % logp2)

print("\nEnd demo ")