%Example 6.9.5
G1=[1 20   21   21   32   42   37   32   32   32;
    1 67   48   29   43   56   48   39   40   41;
    1 37   31   25   27   28   30   31   33   34;
    1 42   40   38   37   36   28   19   27   35;
    1 57   45   32   27   21   25   30   29   29;
    1 39   39   38   46   54   43   31   29   28;
    1 43   32   20   33   46   44   42   37   31;
    1 35   34   34   39   43   39   35   39   42;
    1 41   32   23   37   51   39   27   28   30;
    1 39   32   24   30   35   31   26   29   32];

G2=[2 47   36   25   31   36   29   21   24   27;
    2 53   43   32   40   48   47   46   50   54;
    2 38   35   33   38   42   45   48   48   49;
    2 60   51   41   54   67   60   53   52   50;
    2 37   36   35   40   45   40   34   40   46;
    2 59   48   37   45   52   44   36   44   52;
    2 67   50   33   47   61   46   31   41   50;
    2 43   35   27   32   36   35   33   33   32;
    2 64   59   53   58   62   51   40   42   43;
    2 41   38   34   41   47   42   37   41   46];
A=[2 2 2 -1 -1 -1 -1 -1 -1;0 0 0 1 1 1 -1 -1 -1];
B=[2 -1 -1 2 -1 -1 2 -1 -1 ;0 1 -1 0 1 -1 0 1 -1];
G=[4 -2 -2 -2  1 1 -2  1 1;
    0  2 -2  0 -1 1  0 -1 1;
    0 0 0 2 -1 -1 -2 1 1;
    0 0 0 0 1 -1 0 -1 1];
G1=G1(:,2:size(G1,2))
G2=G2(:,2:size(G2,2))
[n1,p]=size(G1);
[n2,p]=size(G2);
I1=eye(n1);
j1=ones(n1,1);
gbar1=(1/n1)*G1'*j1
J1=ones(size(I1));
S1=(1/(n1-1))*G1'*(I1-(1/n1)*J1)*G1
I2=eye(n2);
j2=ones(n2,1);
gbar2=(1/n2)*G2'*j2
J2=ones(size(I2));
S2=(1/(n2-1))*G2'*(I2-(1/n2)*J2)*G2
gbar=(n1*gbar1+n2*gbar2)/(n1+n2)
E=((n1-1)*S1+(n2-1)*S2)
Spl=E/(n1+n2-2)
H=n1*(gbar1-gbar)*(gbar1-gbar)'+n2*(gbar2-gbar)*(gbar2-gbar)'
% Testing Factor A
T_squareA=(n1+n2)*(A*gbar)'*inv(A*Spl*A')*(A*gbar)
%Testing factor B
T_squareB=(n1+n2)*(B*gbar)'*inv(B*Spl*B')*(B*gbar)
%Testing factor B
T_squareAB=(n1+n2)*(G*gbar)'*inv(G*Spl*G')*(G*gbar)
%
WilksAC = det(A*E*A')/det(A*(E+H)*A')
WilksBC = det(B*E*B')/det(B*(E+H)*B')
WilksABC = det(G*E*G')/det(G*(E+H)*G')
