  G1=[98  81  38  
  103  84  38  
  103  86  42  
  105  86  42  
  109  88  44  
  123  92  50  
  123  95  46  
  133  99  51  
  133  102  51  
  133  102  51  
  134  100  48  
  136  102  49  
  138  98  51  
  138  99  51  
  141  105  53  
  147  108  57  
  149  107  55  
  153  107  56  
  155  115  63  
  155  117  60  
  158  115  62  
  159  118  63  
  162  124  61  
  177  132  67];  
 G2= [93  74  37  
  94  78  35  
  96  80  35  
  101  84  39  
  102  85  38  
  103  81  37  
  104  83  39  
  106  83  39  
  107  82  38  
  112  89  40  
  113  88  40  
  114  86  40  
  116  90  43  
  117  90  41  
  117  91  41  
  119  93  41  
  120  89  40  
  120  93  44  
  121  95  42  
  125  93  45  
  127  96  45  
  128  95  45  
  131  95  46  
  135  106  47  ];

G=[G1;G2];
[n1,p]=size(G1);
[n2,p]=size(G2);
%
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
%
xbar=(n1*xbar1+n2*xbar2)/(n1+n2)
%
Spl=(S1+S2)/(n1+n2-2)
%
W=(n1-1)*S1+(n2-1)*S2
W/(n1+n2-2)
%
B=n1*(xbar1-xbar)*(xbar1-xbar)'+n2*(xbar2-xbar)*(xbar2-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))
