%let program = CHOLINE.SAS;
options   NOCENTER  PS=55 LINESIZE=100;
title1    "PHD Part II Exam of 1998";

proc iml       workspace=2000   symsize=2000 ;

     ***Enter the mean values from Table 1 ***;
               *-Serum Choline--   ---- Plasma ALT-----*;
               *  1  2  3  4  5      1   2   3   4   5 *;
               *----------------   --------------------*;
     MEAN=     { 12 12 12 12 12    .30 .30 .30 .30 .30           /* F+ A */
               , 12 12 12 12 12    .30 .30 .30 .30 .30           /* F+ B */
               , 12 11 10  9 12    .30 .30 .32 .34 .30           /* F+ C */
               , 12 10  9  8 12    .30 .31 .34 .36 .31           /* F+ D */

               ,  9  9  9  9  9    .35 .35 .35 .35 .35           /* F- A */
               ,  9  9  9  9  9    .35 .35 .35 .35 .35           /* F- B */
               ,  9  8  7  7  9    .35 .35 .50 .60 .40           /* F- C */
               ,  9  7  6  6  9    .35 .38 .55 .65 .45           /* F- D */

               , 10 10 10 10 10    .35 .35 .35 .35 .35           /* M+ A */
               , 10 10 10 10 10    .35 .35 .35 .35 .35           /* M+ B */
               , 10  9  8  7 10    .35 .35 .50 .60 .43           /* M+ C */
               , 10  8  6  4 10    .35 .40 .65 .70 .45           /* M+ D */

               ,  9  9  9  9  9    .35 .35 .35 .35 .35           /* M- A */
               ,  9  9  9  9  9    .35 .35 .35 .35 .35           /* M- B */
               ,  9  8  7  6  9    .35 .35 .55 .65 .45           /* M- C */
               ,  9  7  4  3  9    .35 .48 .60 .75 .48 };        /* M- D */

     ROW_ID=   { 'F+   A'
               , 'F+   B'
               , 'F+   C'
               , 'F+   D'

               , 'F-   A'
               , 'F-   B'
               , 'F-   C'
               , 'F-   D'

               , 'M+   A'
               , 'M+   B'
               , 'M+   C'
               , 'M+   D'

               , 'M-   A'
               , 'M-   B'
               , 'M-   C'
               , 'M-   D' };

     ***Enter the variances and correlations from Table 2 ***;
               *-----Serum Choline-----  ------ Plasma ALT--------*;
               *  1   2    3    4    5    1    2    3    4    5   *;
               *-----------------------  -------------------------*;
     VARIANCE= {  2.  2.   2.   2.   2. .009 .009 .009 .009 .009 };

     CORR=     {  1  .57  .46  .51  .63  .18  .29  .15  .19  .18
               ,  .   1   .77  .30  .34  .17  .21  .38  .23  .14
               ,  .   .    1   .29  .34  .18  .24  .26  .31  .13
               ,  .   .    .    1   .25  .13  .19  .22  .17  .32
               ,  .   .    .    .    1   .08  .15  .17  .12  .17
               ,  .   .    .    .    .    1   .80  .56  .54  .59
               ,  .   .    .    .    .    .    1   .56  .54  .56
               ,  .   .    .    .    .    .    .    1   .74  .31
               ,  .   .    .    .    .    .    .    .    1   .30
               ,  .   .    .    .    .    .    .    .    .    1  };
     N=10;
     do i=1 to N-1;  ***to make CORR symmetric***;
        do j=i+1 to N;
           CORR[j,i] = CORR[i,j];
        end;
     end;
     ***Compute Var-Cov matrix "SIGMA" and verify positive-definite ***;
     VARIANCE= { 2  2  2  2  2   .009 .009 .009 .009 .009 };
     SIGMA=    diag(sqrt(VARIANCE)) * CORR * diag(sqrt(VARIANCE)) ;
     call      eigen( EIGENVAL, EIGENVEC, SIGMA );
     print     CORR, VARIANCE, SIGMA, EIGENVAL;

     *** draw a simulated sample (N x p) from a Gaussian distribution ***;
     start sample( MU, SIGMA, SEED );
       HALF=   root( SIGMA );          *** Cholesky decomposition   ***;
       N=      nrow( MU );             *** ( SIGMA = HALF` * HALF ) ***;
       P=      ncol( MU );
       Z=      normal( J(N,P,SEED) );  *** pseudo-random number generator ***;
       Y=      MU +  Z*HALF;
       return(Y);
     finish sample;

     CELLSIZE= 13;
     MU=       J( CELLSIZE, 1, 1 ) @ MEAN;
     SEED=     9491; * STUDENT: CHANGE THIS NUMBER!;
     Y=        sample( MU, SIGMA, SEED );
     ID=       (1:nrow(MU))`;
     INDEX=    J( CELLSIZE, 1, 1 ) @ (1:nrow(ROW_ID))`;
     ROWNAMES= ROW_ID[ INDEX, *];
     ID_Y=     (  ID   ||  Y   );
     VARNAMES= ({"ID"} || {"C1" "C2" "C3" "C4" "C5" "A1" "A2" "A3" "A4" "A5"});
     create    work.FROM_IML
               from ID_Y  [ rowname=ROWNAMES colname=VARNAMES ];
     append    from ID_Y  [ rowname=ROWNAMES colname=VARNAMES ];
     close     work.FROM_IML;
     quit;

data work.FROM_IML;
     set       work.FROM_IML;
     drop      rownames;
     length    stratum $ 2   diet $ 2;
     stratum=  substr(rownames,1,2); * first two characters *;
     diet=     substr(rownames,6,1); * last character *;
     run;

data work.CHOLINE; * transpose rows and columns *;
     set       work.FROM_IML;
     drop      C1-C5 A1-A5;
     period=1; choline=C1; alt=A1; output;
     period=2; choline=C2; alt=A2; output;
     period=3; choline=C3; alt=A3; output;
     period=4; choline=C4; alt=A4; output;
     period=5; choline=C5; alt=A5; output;
     run;

data work.CHOLINE;
     set       work.CHOLINE;
     random1=  rannor( 54574 ); * STUDENT: CHANGE THIS NUMBER!;
     random2=  rannor( 63750 ); * STUDENT: CHANGE THIS NUMBER!;
     run;

proc sort      data=work.CHOLINE;
     by        random1;
     run;
data work.CHOLINE;
     set       work.CHOLINE;
     by        random1;
     if (_N_<=20) then choline=.;
     run;

proc sort      data=work.CHOLINE;
     by        random2;
     run;
data work.CHOLINE;
     set       work.CHOLINE;
     by        random2;
     if (_N_<=20) then alt=.;
     run;

* endsas *;