options nodate nonumber ps=55 ls=80;

title 'RPD -- Various Ch. 2 Examples';

proc iml;



/* Example 2.1 */

T = {1 2, 4 5, 3 0};         /* define a matrix */

W = {-1, 3};                 /* define a vector */

TW = T*W;                    /* matrix multiplication */

print T W TW;                /* displaying matrices */

TTRANS = T`;                 /* transpose */

print TTRANS;



/* Example 2.2: Singular Matrix */

A = {2 4 6, 1 2 3, 5 7 9};

DET_A = det(A);              /* determinant */

print A DET_A;



/* Example 2.3: Matrix Inverse */

B = {1 3 2, 4 5 6, 8 7 9};

DET_B = det(B);

INV_B = inv(B);              /* inverse */

print B DET_B INV_B;



/* Example 2.4: Matrix Rank */

A = {1 2 3, 2 4 6, 3 3 3};

call qr(Q,R,PIV,LINDEP,A);   /* a numerical decomposition */

RANK_A = 3 - LINDEP;         /* LINDEP = smallest dimension - rank */

print A RANK_A;

Y = {6, 10, 9};

A_Y = A||Y;                  /* combine columns into one matrix */

call qr(Q,R,PIV,LINDEP,A_Y);

RANK_A_Y = 3 - LINDEP;

print A_Y RANK_A_Y;

   /* Note:  There is an IML function named "rank",
      but it does something else. */



/* Example 2.10: Projections */

A = (1/6) # {5 2 -1, 2 2 2, -1 2 5};
                             /* scalar multiplication */

AA = A*A;

TRACE_A = trace(A);          /* trace */

print A AA TRACE_A;

I_A = i(3) - A;              /* i(n) = n-by-n identity */

I_AI_A = I_A*I_A;

TRACE_I_A = trace(I_A);

print I_A I_AI_A TRACE_I_A;

AI_A = A*I_A;                /* zero, except for numerical error */

print AI_A;