Multivariate Guassian Distribution

• I want to learn Multivariate Gussian distribution so I written the following code.*I am implementing following formula http://upload.wikimedia.org/math/a/0/a/a0ad9db46854c4a616ce6959095cf21d.png

Iam expecting this type of plot

http://mlwkdb.files.wordpress.com/2009/09/bivariategauss2.jpeg?w=300

Where I am doing mistake?

clear all
clc

% Taking two guassian random variables
x=randn(1000,1);
y=randn(1000,1);

X=[x y];
X=X';
d=size(X,1);

% find means of x,y
mx=mean(x);
my=mean(y);

mumat=[mx my]';
mumat=repmat(mumat,1,size(X,2));

Dif_mat=X-mumat;
% The above step (Dif_mat) is (X-mu) in the formula

cov_mat=cov(X'); % covariance matrix
det_cv=det(cov_mat); % det of cov matrix
inv_cov=inv(cov_mat); % inverse of cov matrix

% scale term before exp in forumala
scale=((2*pi)^(d/2))*sqrt((abs(det_cv))); 
scale=inv(scale);

% Mahabolis distance in formula 
MB=Dif_mat'*cov_mat*Dif_mat;

% find the final probability
p=scale*exp((-1/2)*MB);

surf(x,y,p)

Answer this question please!

Answer 1

[Andrew Newell]

You're really trying to do two things here. The first is, you have some random data and you want to fit it to a multivariate normal distribution. Your approach to this part works, although it can be streamlined:

n = 1000; d=2;
X = randn(n,2);

Get mean and covariance:

mumat=mean(X);
cov_mat=cov(X);

The second part is plotting the resulting distribution. Here you need a regular grid for your variables, not the random values you generated above:

x = -3:.2:3; y = -3:.2:3;
[X,Y] = meshgrid(x,y);
X = X-mumat(1); Y = Y-mumat(2);

Combine X and Y in a way that each row represents one 2D variable.

Z = [X(:) Y(:)];

Now calculate the probabilities.

scale=((2*pi)^(d/2))*sqrt(abs(det(cov_mat))); 
p = zeros(length(Z),1);
for ii=1:length(Z)
  p(ii) = exp(-Z(ii,:)/cov_mat*Z(ii,:)'/2)/scale;
end

Reshape and plot.

p = reshape(p,length(x),length(y));
surf(x,y,p)
xlabel('x'), ylabel('y')

[Raviteja]

Thank You Andrew! This is what I really trying to do.

Thanks alot, you speed up my learning....

[Raviteja]

Here I have a small doubt is,

you are not at all using any value of X=randn(n,2);

you are just using its mean and variances.

Later you are finding probabilities for Z which is the results of meshgrid values

x = -3:.2:3; y = -3:.2:3;

[X,Y] = meshgrid(x,y);

By this code you are changing all the input data given except mean and cov matrix.

[sarath l]

Why can't we use the random numbers generated to plot the graph (instead of using some interval for X,Y)?

Answer 2

[UJJWAL]

Hi RaviTeja,

Unfortunately in my knowledge there is no known single line statement to manually evaluate the multivariate normal distribution. You will have to create a meshgrid and go through the meshgrid to evaluate the function at every pair of x and y and store it correspondingly and then plot it. The problem with the above code are many and as I have told that there is no single line statement, so it is not really important to debug the above code. For Example look at the part where you have defined X. It does not fit in anyway as in the real relation.

In short the straight single statement implementation is not possible. Use a function mvnpdf(). It will find the pdf for you and actually internally it also implements in the same manner.

For more information mail me back or leave a message. Hope this helps

Happy to help

UJJWAL

[UJJWAL]

If you try to understand the reason there is no single line statement, then think like this :-

Look at the equation for the MVN distribution. Now there is the term X-mu. X corresponds to each pair of combination of x and y. However when we evaluate 2-D or 3-D or multivariate functions in MATLAB using meshgrid we basically manipulate the whole meshgrids (Meshgrid basically saves us from the for loops by providing all combinations of x and y). However in this case the situation is slightly different . You need to calculate X-mu and pre-multiply and post-multiply it to the inverst of the covariance matrix.

Now the problem that occours is that when you type X-mu then it considers the whole meshgrid and tends to subtract mu from it and it results in dimension mismatch error. If you put X=[x y] like in your example, you are not basically taking all possible pairs of x and y.so that method does not work. You might think of getting a meshgrid and repmat the mu vector to the size of meshgrid and then proceed but during the evaluation it wont allow you to look at individual entries in meshgrid matrices. It will manipulate the matrices as a whole .So you dont use that approach.

Basically the problem is that during plotting or evaluation MATLAB does not allow you to handle the individual entries in the vectors or matrices which you really need in this case. So you will have to do it by using loops. Or use the mvnpdf which is a much better and elegant way to do it .

[Raviteja]

Thanks Ujjwal, I need some more info.

The thing is, I want to apply the formula given in the link - by using two Gaussian random processes x and y (like in my code). Please can you give me the code.

[UJJWAL]

Belo Andrew has posted a code. He is also essentially doing the same thing as he has to calculate for all the possible pairs of x and y. You have to use the for loop to do that.

I hope that code will be useful. Else mail back.