Bivariate Normal Distribution different for expression and mvnpdf()

Hi,

I am plotting 3D graph of bivariate gaussian distribution using 2 methods for the same range of x and y but still i am getting different results.

Method 1:

    figure6b = figure('Name','Figure 1  ','NumberTitle','off');
    mu = [0 ,0];
    sigma =  [0.5 0.8; 0.8 2.0];
    x1 = -3:0.2:3;
    x2 = -3:0.2:3;
    x = [x1; x2];
    [X1,X2] = meshgrid(x1,x2);
    X = [X1(:) X2(:)];
    y = mvnpdf(X,mu,sigma);
    y = reshape(y,length(x1),length(x2));
    surf(X1,X2,y)
    xlabel('x')
    ylabel('y')
    zlabel('Bivariate Gaussian Distributions')
    
    
    

Output:

Method 2:

    figure6b = figure('Name','Figure 2 ','NumberTitle','off');
    mu = [0 ;0];
    covariance =  [0.5 0.8; 0.8 2.0];
    x1 = -3:0.2:3;
    x2 = -3:0.2:3;
    x = [x1; x2];
    [X1,X2] = meshgrid(x1,x2);
    Z = 1/(2*pi*(det(covariance))^(-0.5))* exp(-(1/2).*(x-mu)'*pinv(covariance)*(x-mu));
    Z = reshape(Z,size(X1));
    surf(X1,X2,Z)
    xlabel('x')
    ylabel('y')
    zlabel('Bivariate Gaussian Distributions')

Output:

Why are the outputs from the two methods different? How can I get the correct output for a bivariate gaussian distribution by the 2nd method? Please help

Answer 1

[David Goodmanson]

Hello Aishwarya,

The problem here is that the array x is 2x31 and does not have every possible combination of an element of X1 and an element of X2, which is 2x961. The method 2 code below fixes that by implementing the same idea as in method 1,

   X = [X1(:) X2(:)];

For method 2 to be correct you need the factor in front to be to the +1/2 power, since you have already put the factor into the denominator.

    mu = [0,0];
    sigma =  [0.5 0.8; 0.8 2.0];
    x1 = -3:0.2:3;
    x2 = -3:0.2:3;
    X = [x1; x2]
    [X1,X2] = meshgrid(x1,x2);
    X = [X1(:) X2(:)];
    y = mvnpdf(X,mu,sigma);
    y = reshape(y,length(x1),length(x2));
    figure(1)
    surf(X1,X2,y)
    xlabel('x')
    ylabel('y')
    zlabel('Bivariate Gaussian Distributions')
   
    % method 2       
    mu = [0;0];
    covar =  [0.5 0.8; 0.8 2.0];
    x1 = -3:0.2:3;
    x2 = -3:0.2:3;
    [X1,X2] = meshgrid(x1,x2);
    X = [X1(:), X2(:)]';
    arg = sum((X-mu).*(pinv(covar)*(X-mu)));
    y2 = 1/(2*pi*(det(covar))^(1/2)) * exp(-(1/2)*arg);
    y2 = reshape(y2,size(X1));
    figure(2)
    surf(X1,X2,y2)
    xlabel('x')
    ylabel('y')
    zlabel('Bivariate Gaussian Distributions')
   
    max(max(abs(y-y2)))
   
    ans = 1.8041e-16

[Aishwarya Radhakrishnan]

Thank you ! That cleared out my doubts conceptually as well!