% Group 1 data
data1=[7 5 9 8
8 8 5 6
16 18 11 9
8 3 7 9
6 3 13 9
11 8 10 10
12 7 9 8
8 11 9 3
14 12 11 4
13 13 13 6
13 9 9 9
13 10 15 7
14 11 12 8
15 11 11 10
13 10 15 9
10 5 8 6
10 3 7 7
17 13 13 7
10 6 10 7
10 10 15 8
14 7 11 5
16 11 12 11
10 7 14 6
10 10 9 6
10 7 10 10
7 6 5 9
15 12 10 6
17 15 15 8
16 13 16 9
13 10 17 8
13 10 17 10
19 12 16 10
19 15 17 11
13 10 7 8
15 11 12 8
16 9 11 11
14 13 14 9];

% Group 2 data
data2=[9 5 10 8
10 0 6 2
8 9 11 1
13 7 14 9
4 0 4 0
4 0 6 0
11 9 9 8
5 3 3 6
9 7 8 6
7 2 6 4
12 10 14 3
13 12 11 10];

[n1, p]=size(data1);
[n2, p]=size(data2);

S1=cov(data1)*(n1-1)/n1;  % sample variance of data1
S2=cov(data2)*(n2-1)/n2;  % sample variance of data2
mean1=mean(data1)';
mean2=mean(data2)';



alpha=0.05; % significance level

S=(n1*S1+n2*S2)/(n1+n2);

%% Parallel Profiles
C=zeros(p-1,p);
for idx_p=1:p-1
    C(idx_p,idx_p)=1; C(idx_p,idx_p+1)=-1;
end

n1*n2*(n1+n2-p)/(n1+n2)^2/(p-1)*(C*(mean1-mean2))'/(C*S*C')*C*(mean1-mean2)
finv(1-alpha,p-1,n1+n2-p)

%% Coincident Profiles
n1*n2*(n1+n2-2)/(n1+n2)^2*(ones(1,p)*(mean1-mean2))^2/(ones(1,p)*S*ones(p,1))
finv(1-alpha,1,n1+n2-2)

%% Level Profiles
mean12=(n1*mean1+n2*mean2)/(n1+n2);
(n1+n2-p)/(p-1)*(C*(mean12))'/(C*S*C')*C*(mean12)
finv(1-alpha,p-1,n1+n2-p)