/* BIVARIATE NORMAL PROBABILITY DENSITY FUNCTION */

nobs = 31;
xaxis = seqa(-3,0.2,nobs);
yaxis = seqa(-3,0.2,nobs);
sigma = {1 .75, .75 1};
binorm = zeros(nobs,nobs);
counter1 = 1;

do while counter1 <= nobs;
  counter2 = 1;
  do while counter2 <= nobs;
	x = xaxis[counter1,1]|yaxis[counter2,1];	
	Q = x'*inv(sigma)*x;
	binorm[counter1,counter2] = ((2*pi)^(-1))*(det(sigma)^(-0.5))*exp(-0.5*Q);
	counter2 = counter2 +1;
  endo;
  counter1 = counter1 + 1;
endo;

library pgraph;
graphset;
title("Sample Bivariate Normal(0,Sigma) Probability Density Function \L Sigma={1 0.75, 0.75 1}"); 
xlabel("X2");
ylabel("X1");
surface(xaxis',yaxis,binorm);