I think this is indeed not very friendly. I will write here the code and explain why it works.
The equation of a multivariate gaussian is as follows:
In the 2D case,
and are 2D column vectors,
So in the 2D case, the vector
is actually a point [x,y], for which we want to compute function value, given the 2D mean vector , which we can also write as [mX, mY], and the covariance matrix .To make it more friendly to implement, let's compute the result of
:So
[CI[0,0] * [x - mX] + CI[1,0] * [y - mY] , CI[0,1] * [x - mX] + CI[1,1] * [y - mY]], where CI is the inverse of the covariance matrix, shown in the equation as
, which is a 2x2 matrix, like is.Then, the current result, which is a 2D row vector, is multiplied [inner product] by the column vector
CI[0,0][x - mX]^2 + [CI[1,0] + CI[0,1]][y - mY][x - mX] + CI[1,1][y - mY]^2
This is going to be easier to implement this expression using NumPy, in comparison to
, even though they have the same value.>>> m = np.array[[[0.2],[0.6]]] # defining the mean of the Gaussian [mX = 0.2, mY=0.6]
>>> cov = np.array[[[0.7, 0.4], [0.4, 0.25]]] # defining the covariance matrix
>>> cov_inv = np.linalg.inv[cov] # inverse of covariance matrix
>>> cov_det = np.linalg.det[cov] # determinant of covariance matrix
# Plotting
>>> x = np.linspace[-2, 2]
>>> y = np.linspace[-2, 2]
>>> X,Y = np.meshgrid[x,y]
>>> coe = 1.0 / [[2 * np.pi]**2 * cov_det]**0.5
>>> Z = coe * np.e ** [-0.5 * [cov_inv[0,0]*[X-m[0]]**2 + [cov_inv[0,1] + cov_inv[1,0]]*[X-m[0]]*[Y-m[1]] + cov_inv[1,1]*[Y-m[1]]**2]]
>>> plt.contour[X,Y,Z]
>>> plt.show[]
The result:
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
import numpy as np | |
import pdb | |
from matplotlib import pyplot as plt | |
from scipy.stats import multivariate_normal | |
def gauss2d[mu, sigma, to_plot=False]: | |
w, h = 100, 100 | |
std = [np.sqrt[sigma[0, 0]], np.sqrt[sigma[1, 1]]] | |
x = np.linspace[mu[0] - 3 * std[0], mu[0] + 3 * std[0], w] | |
y = np.linspace[mu[1] - 3 * std[1], mu[1] + 3 * std[1], h] | |
x, y = np.meshgrid[x, y] | |
x_ = x.flatten[] | |
y_ = y.flatten[] | |
xy = np.vstack[[x_, y_]].T | |
normal_rv = multivariate_normal[mu, sigma] | |
z = normal_rv.pdf[xy] | |
z = z.reshape[w, h, order='F'] | |
if to_plot: | |
plt.contourf[x, y, z.T] | |
plt.show[] | |
return z | |
MU = [50, 70] | |
SIGMA = [75.0, 90.0] | |
z = gauss2d[MU, SIGMA, True] |