G1=[1.11 2.569 3.58 0.760
 1.19 2.928 3.75 0.821
 1.09 2.865 3.93 0.928
 1.25 3.844 3.94 1.009
 1.11 3.027 3.60 0.766
 1.08 2.336 3.51 0.726
 1.11 3.211 3.98 1.209
 1.16 3.037 3.62 0.750]
G2=[ 1.05 2.074 4.09 1.036
 1.17 2.885 4.06 1.094
 1.11 3.378 4.87 1.635
 1.25 3.906 4.98 1.517
 1.17 2.782 4.38 1.197
 1.15 3.018 4.65 1.244
 1.17 3.383 4.69 1.495
 1.19 3.447 4.40 1.026]
G3=[ 1.07 2.505 3.76 0.912
  0.99 2.315 4.44 1.398
 1.06 2.667 4.38 1.197
 1.02 2.390 4.67 1.613
 1.15 3.021 4.48 1.476
 1.20 3.085 4.78 1.571
 1.20 3.308 4.57 1.506
 1.17 3.231 4.56 1.458]
G4=[ 1.22 2.838 3.89 0.944
 1.03 2.351 4.05 1.241
 1.14 3.001 4.05 1.023
 1.01 2.439 3.92 1.067
 0.99 2.199 3.27 0.693
 1.11 3.318 3.95 1.085
 1.20 3.601 4.27 1.242
 1.08 3.291 3.85 1.017]
G5=[ 0.91 1.532 4.04 1.084
 1.15 2.552 4.16 1.151
 1.14 3.083 4.79 1.381
 1.05 2.330 4.42 1.242
 0.99 2.079 3.47 0.673
 1.22 3.366 4.41 1.137
 1.05 2.416 4.64 1.455
 1.13 3.100 4.57 1.325]
G6=[ 1.11 2.813 3.76 0.800
 0.75 0.840 3.14 0.606
 1.05 2.199 3.75 0.790
 1.02 2.132 3.99 0.853
 1.05 1.949 3.34 0.610
 1.07 2.251 3.21 0.562
 1.13 3.064 3.63 0.707
 1.11 2.469 3.95 0.952]
G=[G1;G2;G3;G4;G5;G6];
[n1,p]=size(G1);
[n2,p]=size(G2);
[n3,p]=size(G3);
[n4,p]=size(G4);
[n5,p]=size(G5);
[n6,p]=size(G6);
%
I1=eye(n1);
j1=ones(n1,1);
xbar1=(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);
xbar2=(1/n2)*G2'*j2
J2=ones(size(I2));
S2=(1/(n2-1))*G2'*(I2-(1/n2)*J2)*G2
%
I3=eye(n3);
j3=ones(n3,1);
xbar3=(1/n3)*G3'*j3
J3=ones(size(I3));
S3=(1/(n3-1))*G3'*(I3-(1/n3)*J3)*G3
%
I4=eye(n4);
j4=ones(n4,1);
xbar4=(1/n4)*G4'*j4
J4=ones(size(I4));
S4=(1/(n4-1))*G4'*(I4-(1/n4)*J4)*G4
%
I5=eye(n5);
j5=ones(n5,1);
xbar5=(1/n5)*G5'*j5
J5=ones(size(I5));
S5=(1/(n5-1))*G5'*(I5-(1/n5)*J5)*G5
%
I6=eye(n6);
j6=ones(n6,1);
xbar6=(1/n6)*G6'*j6
J6=ones(size(I6));
S6=(1/(n6-1))*G6'*(I6-(1/n6)*J6)*G6
%
xbar=(n1*xbar1+n2*xbar2+n3*xbar3+n4*xbar4+n5*xbar5+n6*xbar6)/(n1+n2+n3+n4+n5+n6)
%
W=(n1-1)*S1+(n2-1)*S2+(n3-1)*S3+(n4-1)*S4+(n5-1)*S5+(n6-1)*S6
%
B=n1*(xbar1-xbar)*(xbar1-xbar)'+n2*(xbar2-xbar)*(xbar2-xbar)'+n3*(xbar3-xbar)*(xbar3-xbar)'...
    +n4*(xbar4-xbar)*(xbar4-xbar)'+n5*(xbar5-xbar)*(xbar5-xbar)'+n6*(xbar6-xbar)*(xbar6-xbar)';
%
Wilks = det(W)/det(W+B)
H_L_trace=trace(inv(W)*B)
P_trace=trace(B*inv(B+W))
Roy=max(eig(inv(W)*B))
