mu = [0 0]; %Mu vector in the first distribution
sigma = [1 0.9; 0.9 1]; %Sigma matrix in the first distribution
mu2 = [1 1]; %Mu vector in the second distribution
sigma2 = [1 -0.9; -0.9 1]; %Sigma matrix in the second distribution
r = mvnrnd(mu,sigma,100);
r2 = mvnrnd(mu2,sigma2,100);
plot(r(:,1),r(:,2),'b+')
hold on
plot(r2(:,1),r2(:,2),'ro')
title('First Case, sigmax1=1,sigmay1=1,rhoxy1=.9,sigmax2=1,sigmay2=1,rhoxy2=-.9');
clear mu sigma;
disp('press a key to continue');
pause;

mu=[0 0];
sigma=[1 0.5; 0.5 1];
figure;
t = mvnrnd(mu,sigma,100);
subplot(1,2,1),plot(t(:,1),t(:,2),'b*');
title('2-D plot of data points');
x1 = -3:.2:3; x2 = -3:.2:3;
[X1,X2] = meshgrid(x1,x2);
F = mvnpdf([X1(:) X2(:)],mu,sigma);
F = reshape(F,length(x2),length(x1));
subplot(1,2,2), surf(x1,x2,F);
caxis([min(F(:))-.5*range(F(:)),max(F(:))]);
axis([-3 3 -3 3 0 1])   %height 1 for rho=.99 or -.99, height 0.2 for rho=0.5 or 0
xlabel('x1'); ylabel('x2'); zlabel('pdf');
title('2D plot of probability density function');