/* Example 10.3.2 */
proc iml;
x1={ 0,0,0,7,6,0,0,0,1,0,0,1,9,10,0,0,0,0,9,0,0,0,1,4,0,0,0,7,10};
x2={0,0,0,9,10,0,1,1,8,3,1,1,10,10,0,0,0,7,10,1,4,5,5,10,0,2,1,
    10,10};
x3={0,1,3,9,10,3,6,3,10,5,2,3,10,10,0,0,0,8,10,2,4,7,7,10,5, 6,7,
    10,10};
x4={1,3,5,9,10,7,7,7,10,5,2,3,10,10,0,0,3,8,10,4,6,8,7,10,6,7,10,
    10,10};
x5={1,3,5,9,10,7,8,9,10,5,3,3,10,10,0,0,4,8,10,7,8,10,7,10,6,7,10,
    10,10};
d={10,21,30,100,150,40,75,148,256, 27,21,30,90,150,16,43,49,140,
   290,21, 27, 30, 30, 100, 20, 21, 48, 80, 145};
n=29; p=5; q=2; f=n-q; m=2;
y=arsin(sqrt((x1||x2||x3||x4||x5)/10))*180/3.14158;
yb=y`*J(n,n,1/n);
A=log10(d)||J(n,1,1);
B=J(1,5,1)//{1 2 3 4 5};
c={16.8 -6.6 -18 -11.4 19.2, 102.8 -31 -114.8 -88.6 131.6,
  564 -129 -642 -585 792};
w=y`*(I(n)-A*inv(A`*A)*A`)*Y;
v1=y`*A*inv(A`*A)*A`*Y;
u=det(c*w*c`)/det(c*(w+v1)*c`);
chi0=-(f-(p-m-q+1)/2)*log(u);
chia=cinv(1-0.05,q*(p-m));
print   u chi0 chia;
if chi0>chia then print "chi0 > chia: reject";
             else print "chi0 < chia: accept";
/*Try quadratic model*/
m2=3;
B2={1 1 1 1 1, 1 2 3 4 5, 1 4 9 16 25};
c2={-1.2 2.4 0 -2.4 1.2, -14.06 27.43 2.06 -30.17 14.74};
u2=det(c2*w*c2`)/det(c2*(w+v1)*c2`);
chi20=-(f-(p-m2-q+1)/2)*log(u2);
chi2a=cinv(1-0.05,q*(p-m2));
print   u2 chi20 chi2a;
if chi20>chi2a then print "chi20 > chi2a: reject";
             else print "chi20 < chi2a: accept";
psi=inv(B2*inv(w)*B2`)*B2*inv(w)*y`*A*inv(A`*A);
print psi;
quit;
                  U        CHI0        CHIA
          0.3460619   27.589581   12.591587
                  chi0 > chia: reject
                  U2       CHI20       CHI2A
           0.8762065   3.5020685    9.487729
                  chi20 < chi2a: accept
                    PSI
              41.879822    -73.9567
              8.4783284    10.05687
              -1.488394   0.1388218