/* Example 10.3.1 */
proc iml;
y1={251 233 250 141 418 229 271 312 194 211 205 191 249 301 201 277
    294 212 230 246 245 179 165 262 212 285 166 179 298 238 191};
y2={262 218 258 143 371 218 289 323 220 232 299 248 217 270 214 242
    313 236 315 205 192 202 142 274 216 292 171 206 280 267 208};
y3={239 230 258 157 363 228 270 318 214 189 278 283 236 282 247 249
    295 235 300 249 215 194 188 245 228 300 166 214 280 269 162};
y4={234 251 286 162 384 244 296 383 256 230 259 268 266 287 274 293
    295 272 305 225 214 239 192 275 221 319 186 189 328 268 218};
y5={248 273 240 169 387 179 346 310 204 231 266 233 235 268 224 306
    271 287 341 236 242 234 200 278 223 277 220 250 318 280 206};
y6={178 254 185 219 205 182 310 191 245 229 245 240 234 210 275 269
    148 181 165 293 195 210 212 243 259 202 184 238 263 144 220 225
    307 313 206 285};
y7={246 260 232 268 232 213 334 204 270 200 293 313 281 252 231 332
    180 194 242 276 190 230 224 271 279 214 192 272 283 226 272 260
    252 300 177 291};
y8={295 278 215 241 265 173 290 227 209 238 261 251 277 275 285 300
    184 212 250 276 205 249 246 304 296 192 205 297 248 261 222 253
    316 313 194 291};
y9={228 245 220 260 242 200 286 228 255 259 297 307 235 235 238 320
    231 217 249 278 217 240 271 273 262 239 253 282 334 227 246 202
    258 317 194 268};
y10={274 340 292 320 230 193 248 196 213 221 231 291 210 237 251 335
     184 205 312 306 238 194 256 318 283 172 217 251 271 283 253 265
     283 397 212 260};
n1=31;n2=36; p=4; q=2;  m=3;f=n1+n2-q;
/* Testing Adequacy */
y=((y2-y1)//(y3-y1)//(y4-y1)//(y5-y1))||((y7-y6)//(y8-y6)//
   (y9-y6)//(y10-y6));
yb=y*J(n1+n2,n1+n2,1/(n1+n2));
A=((J(1,n1,1)||J(1,n2,0))//(J(1,n1,0)||J(1,n2,1)))`;
B={1  1  1, 1  2  4, 1 3.3333 11.1111, 1  4  16}`;
w=y*(I(n1+n2)-A*inv(A`*A)*A`)*y`;
psi=inv(A`*A)*A`*y`*inv(w)*B`*inv(B*inv(w)*B`);
beta=psi*B;
c={-8 21 -27 14};
v1=y*A*inv(A`*A)*A`*y`;
lambda1=abs(c*w*c`)/abs(c*(w+v1)*c`);
chi0=-(f-(p-m-q+1)/2)*log(lambda1);
chia=cinv(1-0.05,q*(p-m));
pv1=1-probchi(chi0,q*(p-m));
print   psi ,lambda1 ci0 cia pv1;
if chi0>chia then print "chi0 > chia: reject";
             else print "chi0 < chia: accept";
/* Comparison of two groups */
L={1 -1};
G=A*inv(A`*A); E=inv(B*inv(W)*B`);
R=G`*(I(n1+n2)+Y`*inv(W)*y)*G-psi*inv(E)*psi`;
Pp=psi`*L`*inv(L*R*L`)*L*psi;
lambda2=det(E)/det(Pp+E);
f1=(n1+n2-p+m-q); t=1; r1=3;
F0=(f1+1-r1)/r1*(1-lambda2)/lambda2;
Fa=finv(1-0.05, r1,f1+1-r1);
pv2=1-probf(F0,r1,f1+1-r1);
print lambda2 F0 Fa pv2 ;
if F0>Fa then print "F0 > Fa: reject";
             else print "F0 < Fa: accept";
/* Testing constant cholestrol values*/
L1={1 0}; L2={ 0 1 }; M={0 0, 1 0, 0 1};
G=A*inv(A`*A); E=inv(B*inv(W)*B`);
R=G`*(I(n1+n2)+Y`*inv(W)*y)*G-psi*inv(E)*psi`;
Q1=M`*E*M;
Pp1=M`*psi`*L1`*inv(L1*R*L1`)*L1*psi*M;
Pp2=M`*psi`*L2`*inv(L2*R*L2`)*L2*psi*M;
lambda21=det(Q1)/det(Pp1+Q1);
lambda22=det(Q1)/det(Pp2+Q1);
r2=2;
F10=(f1+1-r2)/r2*(1-lambda21)/lambda21;
F20=(f1+1-r2)/r2*(1-lambda22)/lambda22;
F1a=finv(1-0.05, r2,f1+1-r2);
pv21=1-probf(F10, r2, f1+1-r2);
pv22=1-probf(F20, r2, f1+1-r2);
print lambda21 F10 F1a pv21;
if F10>F1a then print "F10 > F1a: reject";
             else print "F10 < F1a: accept";
print lambda22 F20 F1a pv22;
if F20>F1a then print "F20 > F1a: reject";
             else print "F20 < F1a: accept";
quit;
        PSI
  0.8173853   4.8024413   0.4319885
  21.173957   1.8070186   0.0039904
    LAMBDA1       CHI0       CHIA         PV1
  0.9572911   2.8371031   5.9914645   0.2420644
                   CHI0 < CHIA: accept
    LAMBDA2          F0          FA         PV2
  0.9060173   2.1437888   2.7529698   0.1037248
                     F0 < Fa: accept
   LAMBDA21         F10         F1A        PV21
  0.8113408   7.3246216   3.1428085   0.0013801
                    F10 > F1A: reject
   LAMBDA22         F20         F1A        PV22
  0.9804974   0.6265513   3.1428085   0.5377285
                    F20 < F1A: accept