G1=[ 69  75
69  70
71  73
78  82
79  81
73  75];
G2=[ 69  70
68  74
75  80
78  85
68  68
63  68
72  74
63  66
71  76
72  78
71  73
70  73
56  59
77  83];
G3=[72  79
64  65
74  74
72  75
82  84
69  68
76  76
68  65
78  79
70  71
60  61 ];
G=[G1;G2;G3];
[n1,p]=size(G1);
[n2,p]=size(G2);
[n3,p]=size(G3);
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
xbar=(n1*xbar1+n2*xbar2+n3*xbar3)/(n1+n2+n3)
W=(n1-1)*S1+(n2-1)*S2+(n3-1)*S3
B=n1*(xbar1-xbar)*(xbar1-xbar)'+n2*(xbar2-xbar)*(xbar2-xbar)'+n3*(xbar3-xbar)*(xbar3-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))

