proc iml;
start main;

y={400,
   380,
   410,
   350,
   420,
   360,
   390,
   390,
   430,
   370};

  X={1 0 0 1,
     1 0 0 0,
     1 0 0 1,
     0 1 0 0,
     0 1 0 1,
     0 1 0 0,
     0 0 1 1,
     0 0 1 0,
	 0 0 1 1,
	 0 0 1 0};

*Ainv={
      2.5  .5  1   0    0  0  -1  -1   0   0  -1  0  0  0  0  0,
      .5  2.5  0   1    0  0  -1   0   0  -1  0   0 -1  0  0  0,
      1    0   2   0    0  0   0  -1   0   0  -1  0  0  0  0  0,
      0    1   0   3    .5 0   0  .5  -1  -1  0  -1 -1  0  0  0,
      0    0   0  .5   1.5 0   0   0  -1   0  0   0  0  0  0  0,
      0    0   0   0   0  2.5  0  .5   0  .5  0  .5  0 -1  -1 -1,
      -1  -1   0   0   0   0   2   0   0   0  0   0  0  0  0  0,
      -1   0  -1  .5   0  .5   0   3   0   0  0  -1  0 -1  0  0,
	  0    0   0  -1   -1  0   0   0   2   0  0   0  0  0  0  0,
	  0   -1   0  -1   0  .5   0   0   0  2.5  0  0  0  0  -1 0,
	  -1   0  -1   0   0   0   0   0   0   0  2   0  0  0  0  0,
	  0    0   0  -1   0  .5   0  -1   0   0  0  2.5 0  0  0  -1,
	  0   -1   0  -1   0   0   0   0   0   0  0   0  2  0  0  0,
	  0    0   0   0   0  -1   0  -1   0   0  0   0  0  2  0  0,
      0    0   0   0   0  -1   0   0   0   -1 0   0  0  0  2  0,
      0    0   0   0   0  -1   0   0   0   0  0   -1 0  0  0  2};

A={1  0  0  0  0  0  0.5  0.5  0 0 0.5 0.25 0 0.25 0 0.125,
      0  1  0  0  0  0  0.5  0  0 0.5 0 0 0.5 0 0.25 0,
      0  0  1  0  0  0  0  0.5  0 0 0.5 0.25 0 0.25 0 0.125,
      0  0  0  1  0  0  0  0  0.5 0.5 0 0.5 0.5 0 0.25 0.25,
      0  0  0  0  1  0  0  0  0.5 0 0 0 0 0 0 0,
      0  0  0  0  0  1  0  0  0 0 0 0 0 0.5 0.5 0.5,
    0.5  0.5  0  0  0  0  1  0.25  0 0.25 0.25 0.125 0.25 0.125 0.125 
0.0625,
    0.5  0  0.5  0  0  0  0.25  1  0 0 0.5 0.5 0 0.5 0 0.25,
      0  0  0  0.5  0.5  0  0  0  1 0.25 0 0.25 0.25 0 0.125 0.125,
      0  0.5  0  0.5  0  0  0.25  0  0.25 1 0 0.25 0.5 0 0.5 0.125,
    0.5  0  0.5  0  0  0  0.25  0.5  0 0 1 0.25 0 0.25 0 0.125,
   0.25  0  0.25  0.5  0  0  0.125  0.5  0.25 0.25 0.25 1 0.25 0.25 
0.125 0.5,
      0  0.5  0  0.5  0  0  0.25  0  0.25 0.5 0 0.25 1 0 0.25 0.125,
   0.25  0  0.25  0  0  0.5  0.125  0.5  0 0 0.25 0.25 0 1 0.25 0.375,
      0  0.25  0  0.25  0  0.5  0.125  0  0.125 0.5 0 0.125 0.25 0.25 1 
0.3125,
   0.125  0  0.125  0.25  0  0.5  0.0625  0.25  0.125 0.125 0.125 0.5 
0.125 0.375 0.3125 1};

AINV=INV(A);
PRINT AINV;

Z1={  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0,
      0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0,
	  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0,
	  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0,
	  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0,
	  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0,
	  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0,
	  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0,
      0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0,
      0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1};

Z2={  1  0  0  0  0  0,
      0  1  0  0  0  0,
	  0  0  1  0  0  0,
	  1  0  0  0  0  0,
	  0  1  0  0  0  0,
	  0  0  0  1  0  0,
	  1  0  0  0  0  0,
	  0  0  0  1  0  0,
      0  0  0  0  1  0,
      0  0  0  0  0  1};

I={1  0  0  0  0  0,
   0  1  0  0  0  0,
   0  0  1  0  0  0,
   0  0  0  1  0  0,
   0  0  0  0  1  0,
   0  0  0  0  0  1};
 
K11=3.25;
K12=0;
K21=0;
K22=13;

LHS=((X`*X)||(X`*Z1)||(X`*Z2))
//((Z1`*X)||(Z1`*Z1+AINV#K11)||(Z1`*Z2))
//((Z2`*X)||(Z2`*Z1)||(Z2`*Z2+I#K22));

RHS=(X`*Y)//(Z1`*Y)//(Z2`*Y);

C=INV(LHS);

BU=C*RHS;

RMSE=(Y`*Y-BU`*RHS)#(1/6);
print BU ;

finish main;
run;
quit;