/* IML commands for problem 1.6 in the book */

 proc iml ;
 x = { 2  4  3,
       2  3  4,
       3  5  4,
       1  2  6,
       2  6  8 } ;
 n = nrow(x) ;
 p = ncol(x) ;
 *mean = ( x[+,])/n ;  *mean = mean` ;
 mean = (j(1,5)*x)/n;  mean = mean` ;
 sumsq = x`*x -n*mean*mean` ;
 s = sumsq/(n - 1) ;
 sv = sqrt(diag(s)) ;
 r = inv(sv)*s*inv(sv) ;
 za = j(5,1)*mean` ;
 zb = (x -j(5,1)*mean`) ;
 z = (x -j(5,1)*mean`)*inv(sv) ;
 pr1 = x[3,2] ;  pr2 = x[3,]` ;  pr3 = x[3,] ;
  pr4 = z[3,2] ;  pr5 = z[3,]` ;
 print "Results from exercise 1.6 from Proc IML" ;
 print p n pr1 pr2 pr3 ;
 print mean  ;
 print sumsq s r ;
 print za zb z ;
 print pr4 pr5 ;
* These lines follow the example in Chapter 4 ;
 sigma = { 6  2,
           2  3 } ;
 call  eigen(eigsig, evecsig, sigma) ;
 print eigsig evecsig ;
 isig = inv(sigma) ;
 print isig ;
* Computing a distance ;
 distx2 = sqrt( (x[2,] -mean`)*inv(s)*(x[2,]` -mean) ) ;
 print distx2 ;
 quit ;