Multivariate Guassian Distribution

7 views (last 30 days)
Raviteja
Raviteja on 26 Sep 2011
Commented: sarath l on 9 Sep 2018
Iam expecting this type of plot
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!

Accepted Answer

Andrew Newell
Andrew Newell on 26 Sep 2011
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')
  3 Comments
Raviteja
Raviteja on 27 Sep 2011
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
sarath l on 9 Sep 2018
Why can't we use the random numbers generated to plot the graph (instead of using some interval for X,Y)?

Sign in to comment.

More Answers (1)

UJJWAL
UJJWAL on 26 Sep 2011
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
  3 Comments
Raviteja
Raviteja on 26 Sep 2011
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
UJJWAL on 27 Sep 2011
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.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!