dm 'log;clear;output;clear;';
options ps=50 ls=70 pageno=1;
goptions reset=global border ftext=swiss gunit=cm htext=0.4 htitle=0.5;
goptions display noprompt;

**********************************************************************;
**                                                                  **;
** AUTHOR: Chris Bilder                                             **;
** COURSE: STAT 873                                                 **;
** DATE: 1-26-01                                                    **;
** UPDATE: 8-24-03                                                  **;
** PURPOSE: Bivariate normal distribution                           **;
**                                                                  **;
** NOTES:                                                           **;
**                                                                  **;
**********************************************************************;


title1 'Chris Bilder, STAT 873'; 

proc iml;

  mu={15,20};
  sigma={1 0.5,
         0.5 1.25};

  tr_sigma = trace(sigma);
  print 'The trace of sigma' tr_sigma;

  det_sigma = det(sigma);
  print 'The determinant of sigma' det_sigma;

  eig = eigval(sigma);
  print 'The eigenvalues of sigma' eig;

  eigenvec = eigvec(sigma);
  print 'The eigenvectors of sigma' eigenvec;

  *verify eigenvectors;
  check1_left = sigma*eigenvec[,1];
  check1_right = eig[1,1]*eigenvec[,1];

  check2_left = sigma*eigenvec[,2];
  check2_right = eig[2,1]*eigenvec[,2];

  print 'lambda1' check1_left check1_right; 
  print 'lambda2' check2_left check2_right;

  det = det(sigma);
  inv_sigma = inv(sigma);
  pi = 3.141592654;
  p=2;
  save = repeat(0, 30*30*10*10, 3); *initialize save matrix;
  i=1;

  do x1 = 10 to 25 by 0.25;
    do x2 = 10 to 25 by 0.25;
	  x = x1 // x2;
      f_x = (1/(2*pi)**(p/2)) * det**(1/2) * exp(-1/2 * t(x-mu)*inv_sigma*(x-mu));
	  save[i,] = x1 || x2 || f_x; 
	  i=i+1;
    end;
  end;


  col={ "x1" "x2" "f_x"};
  create set1 from save [colname=col];
  append from save;

quit;


data set1;
  set set1;
  if x1=0 and x2=0 then delete;
run;

*Create 3D plot of the surface;
proc g3d data=set1;
  title2 '3D surface ';
  plot x1*x2=f_x / grid zmin=0 zmax=0.3 zticknum=5 rotate=180;
  * rotate=180; *use the rotate option on the above plot statement if you need to; 
run;


*Create contour plot;
proc gcontour data=set1;
  plot x1*x2=f_x  / grid autolabel=(reveal) 
     haxis=axis1 vaxis=axis2 legend=legend1 
	 levels=0.01 0.03 0.05 0.07 0.09 0.11 0.13 0.15;
  title2 'Contour plot';                                
  axis1 label = ('x2')
        length=10
        order = (10 to 30 by 5);
  axis2 label=(a=90 'x1')
        length=8.71
        order = (10 to 30 by 5);
  legend1 label=('f(x)')
          down=2;
run;


*Elipisoid of concentration - try to duplicate book;
proc gcontour data=set1;
  plot x1*x2=f_x  / grid autolabel=(reveal) 
     haxis=axis1 vaxis=axis2 legend=legend1  
	 levels=0.01 0.03 0.05 0.07 0.09 0.11 0.13 0.15;
  title2 'Contour plot';                                
  axis1 label = ('x2')
        length=10
        order = (10 to 30 by 5);
  axis2 label=(a=90 'x1')
        length=8.71
        order = (10 to 30 by 5);
  legend1 label=('f(x)')
          down=2;
run;




quit;