proc iml;

start main;


y={9,

   13,

   4,

   12,

   11,

   11,

   13,

    9,

   10};


X={1,

   1,

   1,

   1,

   1,

   1,

   1,

   1,

   1};


A={1     0      0     0     0.5  0    0.25   0      0.125,

  0     1      0     0     0.5  0    0.25   0      0.125,

  0     0      1     0     0    0.5  0.5    0.25   0.375,

  0     0      0     1     0    0.5  0      0.75   0.375,

  0.5   0.5    0     0     1    0    0.5    0      0.25,

  0     0      0.5   0.5   0    1    0.25   0.75   0.5,

  0.25  0.25   0.5   0     0.5  0.25 1      0.125  0.5625,

  0     0      0.25  0.75  0    0.75 0.125  1.25   0.6875,

  0.125 0.125  0.375 0.375 0.25 0.5  0.5625 0.6875 1.0625};


AINV=INV(A);



lam=1;


Z={1 0 0 0 0 0 0 0 0,

   0 1 0 0 0 0 0 0 0,

   0 0 1 0 0 0 0 0 0,

   0 0 0 1 0 0 0 0 0,

   0 0 0 0 1 0 0 0 0,

   0 0 0 0 0 1 0 0 0,

   0 0 0 0 0 0 1 0 0,

   0 0 0 0 0 0 0 1 0,

   0 0 0 0 0 0 0 0 1};


LHS=((X`*X)||(X`*Z))//((Z`*X)||(Z`*Z+AINV#LAM));

RHS=(X`*Y)//(Z`*Y);

C=INV(LHS);

BU=C*RHS;


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

SIGMA_A=RMSE#(1/LAM);

print BU RMSE SIGMA_A;

CSUB=C[{2 3 4 5 6 7 8 9 10},{2 3 4 5 6 7 8 9 10}];


PEV=CSUB#RMSE;

PV=A*SIGMA_A+PEV;

PRINT PEV PV;

finish main;

run;

quit;