How to find Sigma Matrix for Bivariate Distribution Given Data

Hi,

I'm trying to construct a Bivariate Normal Distribution with Given Data the Location (Lattiude, Longitude).

So far, I have the following code:

xc = 32.7419; % Center of Latitude Coordinate
yc = -117.0904; % Center of Longtitude Coordinate
pop = 72994; % Population at Specified Location

% Creating a 2 dimensional space
x1 = linspace(xc-0.05,xc+0.05,31);
x2 = linspace(yc-0.05,yc+0.05,31);
[X1,X2] = meshgrid(x1,x2);
X = [X1(:) X2(:)];
% Estabilishing Mean Vector,mu & Covariance Matrix,Sigma
mu = [xc yc]
[muHat, Sigma] = normfit(X);
S1 = Sigma(1);
S2 = Sigma(2);
Sigma = [S1 0; 0 S2]
% Sigma = [0.25 0.3; 0.3 1];
% Sigma = [0.85 0.1; 0.1 0.85];

% Evaluating pdf of the normal distribution:
Y = mvnpdf(X,mu,Sigma)*pop;
Y = reshape(Y,length(x2),length(x1));

% Plotting:
figure(1);
surf(x1,x2,Y)
xlabel('x1')
ylabel('x2')
zlabel('Probability Density')
title("Example")

The code produces a distribution, but I would like to create a distribution with minimal standard deviation so that it resembles a "typical" Bivariate Normal Distribution. But I'm having a hard finding a Sigma Matrix that does such.

Answer 1

[Jeff Miller]

You can't really use normfit to estimate sigma in the manner you are doing it, especially because x1 & x2 have such different numerical values (i.e., the standard deviation of your x = [X1(:) X2(:)] vector is far more than the sd of either vector on its own).

I suggest getting rid of normfit and doing something like this:

x1range = 0.1;
x2range = 0.1;
x1 = linspace(xc-x1range/2,xc+x1range/2,31);
x2 = linspace(yc-x2range/2,yc+x2range/2,31);
[X1,X2] = meshgrid(x1,x2);
X = [X1(:) X2(:)];
mu = [xc yc]
% 2 notes about the following line:
%  (a) the constant 8 reflects the fact that the range of a normal
%  distribution is about 8 sd's, from Z=-4 to Z=+4.
%  Admittedly, 6 might be a better choice Z=-3 to Z=+3.
%  (b) mvnpdf wants variances, not sd's
Sigma = [(x1range/8)^2 0; 0 (x2range/8)^2];
Y = mvnpdf(X,mu,Sigma)*pop;

[Thomas Rodriguez]

Thank you so much!

[Thomas Rodriguez]

I just one more question.

Why is that when I plot the same distribution but adjust the x1 and x2 values it shifts the distribution from its original center to a new one?

Is there a way when I plot the distribution remains at the center [xc,yc]?

xmin = 32.5781;
xmax = 33.0175;
ymin = -117.2423;
ymax = -117.0904;
x = linspace(xmin,xmax,31);
y = linspace(ymin,ymax,31);
figure(2)
surf(x,y,Y)
xlabel('x1')
ylabel('x2')
zlabel('Probability Density')
title("Example")

This code shifts the center of the distriubtion from [xc = 32.7419, yc = -117.0904] to [x = 32.7978, -117.1282].

[Jeff Miller]

I am not sure that I understand this question.

Are you trying to plot the original bivariate normal with the original mu and Sigma but just look at a section of that distribution plotted over a slightly different spatial region?

If so, then I think you have to compute a new meshgrid for that new region and compute new Y values corresponding to the meshgrid of that new region (but giving mvnpdf the same mu and Sigma used before, since you want the same center point and spread of the distribution).

[Thomas Rodriguez]

Thank you so much. It worked!

xmin = min(Lat);
xmax = max(Lat);
ymin = min(Lon);
ymax = max(Lon);
x = linspace(xmin,xmax,31);
y = linspace(ymin,ymax,31);
[X1,X2] = meshgrid(x,y);
X = [X1(:) X2(:)];
% Mu:
mu = [xc yc];
% Sigma:
Sigma = [(x1range/8)^2 0; 0 (x2range/8)^2];
% Recompute PDF:
Y = mvnpdf(X,mu,Sigma)*pop;
Y = reshape(Y,length(x),length(y));
figure(2)
surf(x,y,Y)
xlabel('x1')
ylabel('x2')
zlabel('Probability Density')
title("Example")