/*Example 9.8.1----The comparison between PC and LSE  */
proc iml;
xx={10.0     8.9     13.9     14.0     18.8,
    7.4     6.0     10.4     11.1     13.1,
    8.9     5.8     12.5     13.4     12.4,
    6.6     8.4      9.1      9.8     17.2,
    5.7     8.3      8.1      8.7     16.9,
    4.7     9.1      6.7      7.1     18.9,
    9.1     9.5     12.8     13.7     20.0,
    8.7     6.3     12.0     13.0     13.3,
    8.7     6.5     12.1     13.0     14.3,
    6.0     6.8      8.4      9.0     14.7,
    5.8     5.9      8.0      8.6     12.7,
    7.0     5.3      9.9     10.5     11.0,
    7.9     5.8     10.9     11.8     12.7,
    7.9     8.3     10.9     11.7     17.5,
    8.9     5.2     12.4     13.4     10.9,
    9.7     9.9     13.5     14.0     20.0,
   10.0     5.9     14.0     14.0     12.7,
    7.5     8.3     10.4     11.2     17.8,
    7.7     7.1     10.9     11.6     15.1,
    5.0     5.3      6.9      7.5     11.6,
    8.5     4.6     11.8     12.7     10.2,
    4.9     7.4      6.8      7.3     15.4,
    7.3     9.4     10.1     11.0     20.0,
    3.8     8.4      5.3      5.7     17.4,
    5.3     4.9      7.3      7.9     10.2,
    7.5     9.0     10.6     11.3     18.7,
    6.2     8.7      8.7      9.2     18.4,
    8.0     7.6     11.3     12.0     16.1,
    9.9     5.5     13.7     14.0     12.3,
    6.7     6.6      9.4     10.1     13.6};
yy={17.9, 16.2, 16.8, 13.3, 15.8, 12.6, 17.6, 16.8, 16.6, 15.1,
    16.1, 17.0, 17.3, 16.3, 15.7, 17.2, 16.6, 17.1, 16.7, 16.2,
    16.7, 15.4, 15.9, 16.3, 14.9, 17.5, 16.4, 16.1, 17.9, 16.0};
k=5;n=30;
xxb=J(n,n,1/n)*xx;
Z1=xx-xxb;
S1=Z1`*Z1;
D1=sqrt(diag(S1));
U1=Z1*inv(D1);
call eigen(Dlambda,Ga,U1`*U1);
V1=U1*Ga;
print Dlambda  ;
/*So we choose the first two principal components */
/* Following is the comparison between PC and LSE */
sump=0;suml=0;
do i=1 to 30;
if i=1 then  do; x=xx[1:29,]; y=yy[1:29,];end;
       else if i=30 then do;x=xx[2:30,]; y=yy[2:30,];end;
                    else do;x=xx[1:i-1,]//xx[i+1:n,];
                            y=yy[1:i-1,]//yy[i+1:n,];end;
/*PC*/
xb=J(n-1,n-1,1/(n-1))*x;
Z=x-xb;
S=Z`*Z;
D=sqrt(diag(S));
U=Z*inv(D);
call eigen(Dlambdai,Gai,U`*U);
V=U*Gai;
delta=inv(V[,1:2]`*V[,1:2])*V[,1:2]`*y;
yip=J(1,n,1/n)*yy+V1[i,1:2]*delta;
sump=sump+(yy[i,]-yip)**2;
/*LSE*/
eta=inv(U`*U)*U`*y;
beta=inv(D)*eta;
yil=J(1,n,1/n)*yy+(xx[i,]-J(1,n,1/n)*xx)*beta;
suml=suml+(yy[i,]-yil)**2;
end;
DsqPC=sump/30;
DsqLSE=suml/30;
print  DsqPC DsqLSE;
quit;
                     DLAMBDA
                   3.0031544
                     1.98487
                   0.0067562
                   0.0044753
                   0.0007441

                    GA2
              0.5614929  0.133771
              -0.175168 0.6841331
              0.5604725 0.1405369
              0.5601185 0.1383529
              -0.161801 0.6893294

                    DELTA
                  2.0315957
                  0.2647653

                 DSQPC    DSQLSE
              0.970287 1.2052036