/*Example 9.4.3.----- */
proc iml;
y1={12.7,  7.2, 14.8,  5.6, 13.9, 10.4, 16.6, 13.9, 11.9, 10.7, 10.1,
    10.8, 14.7,  9.7,  5.5, 13.1,  8.8,  8.3, 16.2,  7.7, 14.5, 11.6,
     5.3, 11.2, 13.5, 11.6, 10.0, 12.3, 13.6, 12.8, 10.9,  8.3};
y2={10.4,  8.4, 13.9, 10.2, 12.1, 11.2, 17.8, 14.3,  7.0, 11.3,  6.9,
     5.2, 14.4, 12.1,  7.1,  6.5, 17.7,  3.2, 11.9,  4.9, 14.0,  2.5,
     6.1, 11.0,  9.7,  2.4, 14.8, 16.2, 14.9, 13.2, 14.3, 10.5};
y3={ 5.1,  8.5,  8.4,  7.8,  8.8, 12.5, 14.0,  5.9,  5.9,  4.4,  7.8,
     1.3, 11.6,  5.2,  7.8,  1.3, 11.5,  5.2, 10.2, 11.7, 16.9,  5.5,
     1.5,  5.7, 12.3,  9.7,  9.1,  6.6,  9.3, 11.6, 10.8,  7.5};
y4={ 8.6,  6.8,  7.3,  6.1,  8.8,  7.0,  6.8,  7.7,  6.1,  3.9,  6.4,
     1.3,  7.4,  9.1,  3.1,  0.9,  6.6,  8.9, 15.6, 12.7,  8.4,  4.5,
     4.1,  8.1, 13.4, 14.0,  9.6,  9.2, 10.2, 11.5,  9.6, 10.6};
y5={ 7.1,  6.3,  8.0,  6.4,  8.1,  6.9,  5.9,  9.2,  0.8,  4.7,  9.5,
     2.1,  7.8,  9.7,  1.5,  0.8,  5.4,  4.3, 15.3, 13.2, 13.1,  5.8,
     6.2, 10.0, 14.0, 10.8,  8.2,  8.3, 11.5, 11.1, 13.2,  8.5};
y6={ 9.7,  4.5, 10.4, 16.5,  7.8,  6.9,  5.3,  5.7,  5.1,  5.3,  7.9,
     6.3, 14.8,  5.2,  8.4,  0.5, 12.0,  4.4, 13.9, 10.7,  9.8,  8.6,
     2.1,  8.1,  6.1, 10.3,  9.3, 12.6, 15.8, 10.5, 10.0,  6.2};
y7={ 3.0,  7.9,  6.0, 16.7,  4.8,  6.0,  3.7,  7.3,  4.2,  6.6,  5.7,
     6.6,  9.7,  8.8, 12.6,  0.9,  2.3,  8.1,  5.6,  9.0,  6.6,  6.0,
    11.7, 11.2,  3.3,  5.7,  7.2,  8.1,  2.7,  6.2,  6.0, 13.1};
y8={ 4.9,  7.0,  7.4, 12.8,  6.5,  7.3,  2.7,  6.1,  9.7,  8.1,  9.8,
     8.9,  9.9,  6.5, 13.9, 10.7,  0.6, 11.8,  2.0,  7.5,  6.0, 16.4,
    14.7, 12.1,  1.2, 15.5,  4.1,  4.8,  5.6,  6.4,  2.2,  4.1};
y9={ 3.8,  7.7,  7.2, 11.6,  4.1,  1.2,  1.0,  9.2, 10.1,  9.9,  4.7,
    12.9,  8.9,  7.3,  4.2, 13.4,  1.0,  7.3,  0.0,  0.7,  9.7,  9.9,
    13.4,  9.2,  1.9,  3.8,  7.0,  7.9,  2.7,  2.5,  2.0,  5.1};
y10={5.5,  8.7, 10.5, 15.9,  7.1,  9.0,  7.7,  9.7,  7.7,  8.8,  5.9,
    15.0, 15.5,  5.7, 11.1,  7.9,  7.0,  7.0,  0.8,  1.4,  9.8,  8.8,
     7.0,  6.6,  1.3,  1.8,  3.4,  8.2,  4.3,  4.1,  7.0,  6.2};
y11={6.3,  7.1, 12.7, 10.6,  7.0, 10.8,  7.4,  7.3, 15.3, 11.4,  4.5,
    13.0, 11.0,  6.0, 16.0,  9.8,  9.8, 10.3,  3.6,  0.4,  5.1,  9.4,
     4.2,  8.9,  2.8,  6.3,  6.5,  9.9,  3.0,  4.9,  4.1,  8.3};
y12={5.3, 11.4,  8.9,  4.5,  7.6,  8.1, 11.5,  9.0, 10.9, 10.1,  7.6,
     8.6,  5.0, 12.2,  8.9, 12.9,  4.8,  9.5,  8.9,  0.3,  6.0,  8.4,
    13.3, 10.4,  5.8,  1.5,  1.3,  5.2,  1.3,  4.0,  4.2,  6.0};
x1={    0.25,     0.25,     0.25,     0.25,     0.25,     0.25,
        0.25,     0.25,     0.65,     0.65,     0.65,     0.65,
        0.65,     0.65,     0.65,     0.65,     1.30,     1.30,
        1.30,     1.30,     1.30,     1.30,     1.30,     1.30,
        1.71,     1.71,     1.71,     1.71,     1.71,     1.71,
        1.71,     1.71};
x2={  0.0625,   0.0625,   0.0625,   0.0625,   0.0625,   0.0625,
      0.0625,   0.0625,   0.4225,   0.4225,   0.4225,   0.4225,
      0.4225,   0.4225,   0.4225,   0.4225,   1.6900,   1.6900,
      1.6900,   1.6900,   1.6900,   1.6900,   1.6900,   1.6900,
      2.9241,   2.9241,   2.9241,   2.9241,   2.9241,   2.9241,
      2.9241,   2.9241};
x3={0.015625, 0.015625, 0.015625, 0.015625, 0.015625, 0.015625,
    0.015625, 0.015625, 0.274625, 0.274625, 0.274625, 0.274625,
    0.274625, 0.274625, 0.274625, 0.274625, 2.197000, 2.197000,
    2.197000, 2.197000, 2.197000, 2.197000, 2.197000, 2.197000,
    5.000211, 5.000211, 5.000211, 5.000211, 5.000211, 5.000211,
    5.000211, 5.000211};
x4={     254,      262,      301,      311,      290,      300,
         306,      286,      275,      282,      256,      276,
         337,      296,      309,      296,      275,      192,
         338,      248,      315,      295,      312,      286,
         275,      270,      290,      260,      302,      284,
         280,      329};
y=y1||y2||y3||y4||y5||y6||y7||y8||y9||y10||y11||y12;
x=J(32,1,1)||x1||x2||x3||x4;
beta=inv(x`*x)*x`*y;
H=x*inv(x`*x)*x`;
c={0 0 0 1 0}; n=32; p=12; q=5; f=n-q;
v=y`*(I(n)-H)*y;
W=beta`*c`*inv(c*inv(x`*x)*c`)*c*beta;
u=det(v)/det(v+w);
F0=(f+1-p)/p*(1-u)/u;
Fa=finv(1-0.05, p, f+1-p);
print  F0 Fa;
if F0>Fa then print "F0 > Fa: reject";
         else print "F0 < Fa: accept";
quit;
             F0            FA
      0.5766299       2.42466
            F0 < Fa: accept