proc iml;
y1={21,21,20.5,23.5,21.5,20,21.5,23,20,16.5,24.5};
y2={20,21.5,24,24.5,23,21,22.5,23,21,19,25};
y3={21.5,24,24.5,25,22.5,21,23,23.5,22,19,28};
y4={23,25.5,26,26.5,23.5,22.5,25,24,21.5,19.5,28};
M=11; p=4;  f=M-1; m=2;
y=y1||y2||y3||y4;
yb=y`*J(11,1,1/11);
s= (y`-yb*J(1,11,1))*(y`-yb*J(1,11,1))`/(11-1);
print s;
X={ 1   1   1   1,
   -3  -1   1   3};
G={ 1  -1  -1   1,
  -1   3  -3   1};
b_egls=inv(X*inv(s)*X`)*X*inv(s)*Yb;
print b_egls;
e=inv(X*inv(s)*X`);
print e;
Tsq=11*yb`*G`*inv(G*s*G`)*G*yb;
print Tsq;
U={0 1}; k=1; d=U*b_egls;
print x;
F0=(f-k-p+m+1)*11/(f*k)*d`*inv(U*inv(X*inv(s)*X`)*U`)*d/(1+Tsq/f);
print F0;
Fa=Finv(1-0.05,k,f-k-p+m+1);
print Fa;