/* EXAMPLE 7.3.2 */

FILENAME PSYCH 'T5_1_PSYCH.dat';
TITLE 'EXAMPLE 7.3.2';
DATA SCORES;
  INFILE PSYCH;
  INPUT SEX TEST1 TEST2 TEST3 TEST4;
PROC IML;
  USE SCORES;
  READ ALL VAR {TEST1 TEST2 TEST3 TEST4} INTO X;
  N = NROW(X);
  P = NCOL(X);
  X1 = X[1:32,];
  X2 = X[33:64,];
  N1 = NROW(X1);
  X1BAR = 1/N1*X1`*J(N1,1);
  S1 = 1/(N1-1)*X1`*(I(N1)-1/N1*J(N1))*X1;
  N2 = NROW(X2);
  X2BAR = 1/N2*X2`*J(N2,1);
  S2 = 1/(N2-1)*X2`*(I(N2)-1/N2*J(N2))*X2;
  Spl = 1/(N1+N2-2)*((N1-1)*S1+(N2-1)*S2);
  v1 = N1-1;
  v2 = N2-1;
  v = v1+v2;
  detS1 = det(S1);
  detS2 = det(S2);
  detSpl = det(Spl);
  lnM =0.5*(v1*LOG(DET(S1))+v2*LOG(DET(S2)) - v*LOG(DET(Spl)));
  K = 2;
  c1 = (1/v1+1/v2-1/v)*(2*P**2+3*P-1)/(6*(P+1)*(K-1));
  u = -2*(1-c1)*lnM;
  c2 = (1/v1**2+1/v2**2-1/v**2)*(P-1)*(P+2)/(6*(K-1));
  a1 = .5*(K-1)*P*(P+1);
  a2 = (a1+2)/ABS(c2-c1**2);
  b1 = (1-c1-a1/a2)/a1;
  b2 = (1-c1-2/a2)/a2;
  PRINT v1,v2,detS1,detS2,detSpl,lnM,c1,u,c2,a1,a2,b1,b2;
  START FVALUE;
    IF c2>c1**2 THEN DO;
      F = -2*b1*lnM;
      PRINT F;
      END;
    ELSE DO;
      F = -a2*b2*lnM/(a1*(1+2*b2*lnM));
      PRINT F;
      END;
   FINISH;
  RUN FVALUE;
ENDSAS;