# Plot two-dimensional Gaussian density function in MATLAB

I am learning ML. I have encountered a lot probability density functions. I want to use MATLAB to illustrate them. For example,the density function of the Normal distribution function in two dimensions:

$$p(\vec{x};\mu,\Sigma) = \frac{1}{2\pi \sqrt{\det(\Sigma)}}\exp\left(-\frac{1}{2}(\vec{x}-\mu)^T\Sigma^{-1}(\vec{x}-\mu)\right)$$

When we use MATLAB to plot three-dimensional graphs, it plots a lot of points in space, and then uses straight lines to connect adjacent points. The steps are:

1. generating points of the domain in the XY plane: [X,Y]=meshgrid(-3;.5;3)

2. using vectors $X$ and $Y$ to generate $Z$. If the function is $z= \sin(xy)$, it is very simple. Z=sin(X.*Y) is OK. But here we need to do $(x,y)\cdot\Sigma^{-1}\cdot(x,y)^T$ for every element in $[X,Y]$. If I do it directly: X.*$\Sigma^{-1}$.*Y. I will encounter the error that the dimensions do not match, because the matrix $\Sigma$ is $2 \times 2$. and the vector $X$ is about $12 \times 1$.

How can I solve it in MATLAB? Thanks.

## locked by whuber♦Oct 7 '15 at 17:59

This question exists because it has historical significance, but it is not considered a good, on-topic question for this site so please do not use it as evidence that you can ask similar questions here. This question and its answers are frozen and cannot be changed. See the help center for guidance on writing a good question.

Read more about locked posts here.

## 2 Answers

It is not terribly enlightening, I am afraid, but you can permute the parallel dimensions out to the third and fourth, use bsxfun to do the broadcasting, then bring the results back. It involves a few shiftdim operations:

% set the random seed:
randn('state',[2979743726;1541610269]);

% create a random sigma, mu:
Sigma = cov(randn(100,2));
mu = randn(2,1);

% compute, plot the cdf
[X,Y] = meshgrid(-3:.2:3);
XY = cat(3,X,Y);
% subtract mu
XYmmu = bsxfun(@minus,XY,shiftdim(mu(:),-2));

isigXY = squeeze(sum(bsxfun(@times,shiftdim(inv(Sigma),-2),XYmmu),3));
XYisXY = sum(isigXY .* XYmmu,3);

normcdf = (1/(2*pi*sqrt(det(Sigma)))) * exp(-0.5 * XYisXY);
surf(X,Y,normcdf); To do that operation I would try something along the lines of

product = [X(:), Y(:)] * ( SigmaInv * [X(:)'; Y(:)'] )';


where SigmaInv is your 2x2 matrix. The variable product should be a matrix having the values for all the points of the grid.