bone=[...
    1  47.8 48.8 49.0 49.7
2  46.4 47.3 47.7 48.4
3  46.3 46.8 47.8 48.5
4  45.1 45.3 46.1 47.2
5  47.6 48.5 48.9 49.3
6  52.5 53.2 53.3 53.7
7  51.2 53.0 54.3 54.4
8  49.8 50.0 50.3 52.7
9  48.1 50.8 52.3 54.4
10 45.0 47.0 47.3 48.3
11 51.2 51.4 51.6 51.9
12 48.5 49.2 53.0 55.5
13 52.1 52.8 53.7 55.0
14 48.2 48.9 49.3 49.8
15 49.6 50.4 51.2 51.8
16 50.7 51.7 52.7 53.3
17 47.2 47.7 48.4 49.5
18 53.3 54.6 55.1 55.3
19 46.2 47.5 48.1 48.4
20 46.3 47.6 51.3 51.8];
bone=bone(:,2:end);
[n1,p]=size(bone);
I1=eye(n1);
j1=ones(n1,1);
ybar1=(1/n1)*bone'*j1
J1=ones(size(I1));
S1=(1/(n1-1))*bone'*(I1-(1/n1) *J1)*bone

MLS=zeros(p,p);
for i=1:n1
    MLS=(bone(i,:)'-ybar1)*(bone(i,:)'-ybar1)'+MLS;
end
MLS=(1/n1).*MLS

g=zeros(n1,n1);

for i=1:n1
    for j=1:n1
    g(i,j)=[(bone(i,:)'-ybar1)'*inv(MLS)*(bone(j,:)'-ybar1)];
    end
end

sum1=sum(g.^3,1);
sum2=sum(sum1);
b1p=(1/n1^2)*sum2
b2p=sum(diag(g).^2)*(1/n1)
z1=(((p+1)*(n1+1)*(n1+3))/(6*((n1+1)*(p+1)-6)))*b1p
z2=(b2p-p*(p+2))/sqrt(8*p*(p+2)/n1)
z3=((b2p-p*(p+2)*(n1+p+1)/n1))/sqrt(8*p*(p+2)/(n1-1))
