G1=[ 5.4  6.0  6.3  6.7 
5.2  6.5  6.0  5.8 
6.1  5.9  6.0  7.0 
4.8  5.0  4.9  5.0 
5.0  5.7  5.0  6.5 
5.7  6.1  6.0  6.6 
6.0  6.0  5.8  6.0 
4.0  5.0  4.0  5.0 
5.7  5.4  4.9  5.0 
5.6  5.2  5.4  5.8 
5.8  6.1  5.2  6.4 
5.3  5.9  5.8  6.0 
]
G2=[5.0  5.3  5.3  6.5 
4.8  4.9  4.2  5.6 
3.9  4.0  4.4  5.0 
4.0  5.1  4.8  5.8 
5.6  5.4  5.1  6.2 
6.0  5.5  5.7  6.0 
5.2  4.8  5.4  6.0 
5.3  5.1  5.8  6.4 
5.9  6.1  5.7  6.0 
6.1  6.0  6.1  6.2 
6.2  5.7  5.9  6.0 
5.1  4.9  5.3  4.8 
]
G3=[4.8  5.0  6.5  7.0 
5.4  5.0  6.0  6.4 
4.9  5.1  5.9  6.5 
5.7  5.2  6.4  6.4 
4.2  4.6  5.3  6.3 
6.0  5.3  5.8  6.4 
5.1  5.2  6.2  6.5 
4.8  4.6  5.7  5.7 
5.3  5.4  6.8  6.6 
4.6  4.4  5.7  5.6 
4.5  4.0  5.0  5.9 
4.4  4.2  5.6  5.5 
]
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))
