G1=[...   
    1  112.5  100.5  102.5  102.5  107.5  107.5
    1   92.5  102.5  105.0  100.0  110.0  117.5
    1  132.5  125.0  115.0  112.5  110.0  110.0
    1  102.5  107.5  107.5  102.5  90.0  112.5
    1  110.0  130.0  115.0  105.0  112.5  110.0
    1   97.5   97.5   80.0   82.5   82.5  102.5
    1   90.0   70.0   85.0   85.0   92.5   97.5];
G2=[...
    2  115.0  115.0  107.57 107.5  112.5  107.5
    2  125.0  125.0  120.0  120.0  117.5  125.0
    2   95.0   90.0   95.0   90.0  100.0  107.5
    2   87.5   65.5   85.0   90.0  105.0   90.0
    2   90.0   87.5   97.5   95.0  100.0   95.0
    2   97.5   92.5   57.5   55.0   90.0   97.5
    2  107.5  107.5  145.0  110.0  105.0  112.5
    2  102.5  130.0   85.0   80.0  127.5   97.5 ];
G3=[...
    3  107.5  107.5  102.5  102.5  102.5   97.5
    3   97.5  108.5   94.5  102.5  102.5  107.5
    3  100.0  105.0  105.0  105.0  110.0  110.0
    3   95.0   95.0   90.0  100.0  100.0  100.0
    3   85.0   92.5   92.5   92.5   90.0  110.0
    3   82.5   77.5   75.0   65.5   65.0   72.5
    3   62.5   75.0  115.0  110.0  100.0  100.0
    ];
G4=[...
    4   70.0   67.5   67.5   77.5   77.5   77.5
    4   45.0   37.5   45.0   45.0   47.5   45.0
    4   52.5   22.5   90.0   65.0   60.0   65.5
    4  100.0  100.0  100.0  100.0   97.5   92.5
    4  115.0  110.0  100.0  110.0  105.0  105.0
    4   97.5   97.5   97.5  105.0   95.0   92.5
    4   95.0  125.0  130.0  125.0  115.0  117.5
    4   72.5   87.5   65.0   57.5   92.5   82.5
    4  105.0  105.0  105.0  105.0  102.5  100.0
    ];
G1=G1(:,2:7);
G2=G2(:,2:7);
G3=G3(:,2:7);
G4=G4(:,2:7);
%
[n1,p]=size(G1);
[n2,p]=size(G2);
[n3,p]=size(G3);
[n4,p]=size(G4);
%
I1=eye(n1);
j1=ones(n1,1);
gbar1=(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);
gbar2=(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);
gbar3=(1/n3)*G3'*j3;
J3=ones(size(I3));
S3=(1/(n3-1))*G3'*(I3-(1/n3)*J3)*G3;
%
I4=eye(n4);
j4=ones(n4,1);
gbar4=(1/n4)*G4'*j4;
J4=ones(size(I4));
S4=(1/(n4-1))*G4'*(I4-(1/n4)*J4)*G4;
%
gbar=(n1*gbar1+n2*gbar2+n3*gbar3+n4*gbar4)/(n1+n2+n3+n4);
W=(n1-1)*S1+(n2-1)*S2+(n3-1)*S3+(n4-1)*S4/(n1+n2+n3+n4)-4;
%
SX=[G1;G2;G3;G4];
IS=eye(size(SX,1));
JS=ones(size(IS));
S2=(1/(size(SX,1)-1))*SX'*(IS-(1/size(SX,1))*JS)*SX;
%
tM30=([1 5 10 15 30 60]'-(ones(6,1))*30)
A=[ones(6,1) tM30 tM30.^2 tM30.^3 tM30.^4]
betahat=inv(A'*inv(S2)*A)*A'*inv(S2)*gbar
T21=28*(gbar-A*betahat)'*inv(S2)*(gbar-A*betahat)
T22=28*(gbar'*inv(S2)*gbar-gbar'*inv(S2)*A*betahat)
%
minutes=[1 5 10 15 30 60];
axis([0 60 80 110])
hold on
plot(minutes,gbar1,':')
text('Position',[minutes(6:6),gbar1(6:6)],'String','Group 1')
plot(minutes,gbar2,'-.')
text('Position',[minutes(6:6),gbar2(6:6)],'String','Group 2')
plot(minutes,gbar3,'-.')
text('Position',[minutes(6:6),gbar3(6:6)],'String','Group 3')
plot(minutes,gbar4,'-.')
text('Position',[minutes(6:6),gbar4(6:6)],'String','Group 4')
plot(minutes,gbar,'-')
text('Position',[minutes(6:6),gbar(6:6)],'String','All groups')
title ('Profiles and mean profile for four rat groups')
ylabel('Weight (in grams)')
xlabel('Minutes')
hold off


