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

title1 'Chris Bilder, STAT 873';

*******************************************************************;
**                                                               **;
** NAME:  Chris Bilder                                           **;
** COURSE: STAT 873                                              **;
** DATE:  1-07-01                                                **;
** UPDATE: 8-24-03                                               **;
** Purpose: Generate data from a bivariate normal                **; 
**          distribution                                         **;
**                                                               **;
** NOTES:  1) See http://ftp.sas.com/techsup/download/stat/mvn.html **;
**            for an additional SAS program to do this.          **;
**         2) This program can be used to for multivariate normal**;
**            distributions by changing the mu vector and sigma  **;
**            matrix.                                            **;
**                                                               **;
*******************************************************************;


proc iml;
  
  mu={15,20};
  sigma={1 0.5,
      0.5 1.25};
  *Note that the pop. corr. is 0.5/sqrt(1*1.25) = 0.45;
  
  *CALL VNORMAL( series, mu, sigma, n, seed);
  call vnormal(save, mu, sigma, 1000, 155504946);

  *print save;

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

quit;


title2 'Random sample from bivariate normal distribution';
proc print data=set1;
run;


title2 'Examine the estimated correlation';
proc corr data=set1;
  var x1 x2;
run;


title2 'Scatter plot of the data';
proc gplot data=set1;
  plot x1*x2 / grid vaxis=axis1 haxis=axis2 frame;
  title2 'x1 vs. x2';
  symbol1 v=dot cv=blue h=0.1;
  axis1 label = (a=90 'x1')
        length = 10
        order = (10 to 30 by 5);
  axis2 label = ('x2')
        length = 8.71
        order = (10 to 30 by 5);
run;


*Finds a bivariate estimate of the density and puts it into the out data set;
*Note the following from the help for PROC KDE:                             ; 
*  PROC KDE uses a Gaussian density as the kernel, and its assumed variance determines;
*  the smoothness of the resulting estimate.;
*Since a Gaussian (normal) density is used, the resulting bivariate density estimate;
*  will probably be closer to the bivariate normal density that what it actually is;
proc kde data=set1 out=set2;
  var x1 x2;
run;


*Create 3D plot of the surface;
proc g3d data=set2;
  title2 '3D surface ';
  plot x1*x2=density / grid zmin=0 zmax=0.3 zticknum=5;
run;


*Create contour plot;
proc gcontour data=set2;
  plot x1*x2=density  / 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;


*To help put the above 3D surface plot on the same scale;
*  add this data set to end of set2.  MAKE SURE that these;
*  min x1 and x2 value in set is >10 and max x1 and x2 value;
*  is < 25!;
data extra;
  input x1 x2 density;
  datalines;
  10 10 0
  10 25 0
  25 10 0
  25 25 0
  ;
run;


data set3;
  set set2 extra;
run;


*Create 3D plot of the surface that is same scale as chapter 1 graph;
proc g3d data=set3;
  title2 '3D surface (same scale as chapter 1 graph)';
  plot x1*x2=density / grid zmin=0 zmax=0.3 zticknum=5;
run;

quit;