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-8-01                                                     **;
** UPDATE: 2-6-01 - corrected error in finding f(x), 8-22-03        **;
** PURPOSE: Graph multivariate normal distribution                  **;
**                                                                  **;
** NOTES:                                                           **;
**                                                                  **;
**********************************************************************;


title1 'Chris Bilder, STAT 873'; 

proc iml;

  pi = 3.141592654;
  mu={15,20};
  sigma={1 0.5,
         0.5 1.25};
  *sigma={1 0.5,
          0.5 3};
  *sigma={1 0,
          0 3};
  *sigma={1 0,
          0 1};
  *sigma={1 1,
         1 0.5};
 
  det = det(sigma);          
  inv_sigma = inv(sigma);
  p=2;

  *initialize save matrix, note there are 61*61 x1 and x2 combinations in the do loop below;
  save = repeat(0, 3721, 3); 

  *initialize counter for do loop;
  i=1;  

  *evaluate f(x) for many x1 and x2 values;
  do x1 = 10 to 25 by 0.25;
    do x2 = 10 to 25 by 0.25;

	  *Creates 2x1 vector;
	  x = x1 // x2; 

      *evaluate f(x);
      f_x = (1/(2*pi)**(p/2)) * det**(-1/2) * exp(-1/2 * t(x-mu)*inv_sigma*(x-mu)); 

	  *Create 1x3 vector for the results and then save it into the ith row of the save matrix;
	  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;


*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=50;
  * 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; 
  *The reveal statement puts the level labels on the plot;
  *  If there are too many to fit, SAS will not plot them;
  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;