PROC IML;

sigma_a2 =1;                        
sigma_e2 =2;                        

y={ 2.3, 2.2,  1.7, 1.8,  1.9, 1.9,
               1.8, 1.7,  1.9, 1.7};
                                    

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

                                   
Z={ 1 0 0 0 0,
    1 0 0 0 0,
    0 1 0 0 0,
    0 1 0 0 0,
    0 0 1 0 0,
    0 0 1 0 0,
    0 0 0 1 0,
    0 0 0 1 0,
    0 0 0 0 1,
    0 0 0 0 1};
                                   

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

G = A     * sigma_a2;
R = I(10) * sigma_e2;

PRINT y,X,Z;
PRINT G, R;

V  = Z*G*Z`  + R;                    
IV = INV(V);
b  = INV (X`*IV*X) *  (X`*IV*y);     

a  = G*Z`*IV*(y-X*b);                
e  = R   *IV*(y-X*b);                

PRINT b, a, e;

IR   = INV(R);
XRX  = X`*IR*X;
XRZ  = X`*IR*Z;
ZRX  = XRZ`;
ZRZG = Z`*IR*Z + INV(G);
                            

C1   = XRX || XRZ;          
C2   = ZRX || ZRZG;         
C    = C1  // C2;           
IC   = INV(C);              
                            

XRy  = X`*IR*y;             
ZRy  = Z`*IR*y;             
                            

RHS  = XRY//ZRy;            
                            

SOL  = IC * RHS;            
                            

b = SOL(|1:2 ,|);
a = SOL(|3:7 ,|);

e = y - X*b - Z*a;

PRINT  XRX, XRZ, ZRX, ZRZG,     XRy, ZRy;
PRINT  RHS,  SOL, a,  b;