%% computes and displays bivariate normal densities
% Also relies on the function MVN in the separate file, MVN.m
% 

corr = -0.99; % correlation coefficient between X1 and X2. 
G = ones(2,2);   % set Var(X1) = Var(X2) = 1
G(1,2) = corr;
G(2,1) = corr;

Ginv = inv(G);
detG = det(G);
D = size(G);
N = D(1);
factor = (2*pi)^(N/2) *sqrt(det(G))  % miltiplier in front of the exp(...) in the MVN density


K = 20;
D = 2*K+1;
fGrid = ones(D,D);
for xi=-K:K,
 for yj = -K:K,  
  x = [xi*2/K, yj*2/K]';    
  [f] = MVN(x,Ginv,factor);
  fGrid(xi+1+K,yj+1+K) = f;  
 end
end
xvec = (-K:K)*2/K;
surf(xvec,xvec,fGrid)