load schizo.dat
grp=schizo(:,1);
Y=schizo(:,2:5);
ii=find(grp==1);
Y=sqrt(Y(ii,:));
N=size(Y,1)
X=ones(N,1);
S=(Y'*Y-Y'*ppo(X)*Y)/(N-1);
ybar=mean(Y)';

plot([0 15 30 45]',ybar,'-')
axis([-10 50 4.6 6])
set(gca,'FontSize',16)
xlabel('Time after Insulin')
ylabel('Root fatty Acid')
set(gca,'LineWidth',2)
hold on
plot([0 15 30 45]',ybar,'*','MarkerSize',10)
hold off
print -deps profile1m.eps

T2=N*(ybar-7)'*inv(S)*(ybar-7)
d=4;
F=T2*(N-d)/(d*(N-1))
crit=finv(.95,d,N-d)
pval = 1-fcdf(F,d,N-d)

% Profile Level
M=ones(4,1)/4
psi=M'*ybar;
T2=N*psi'*inv(M'*S*M)*psi
d=1;
F=T2*(N-d)/(d*(N-1))
crit=finv(.95,d,N-d)
pval = 1-fcdf(F,d,N-d)
se2=diag(M'*S*M/N);
crit=sqrt(d*(N-1)/(N-d)*finv(.95,d,N-d))
[psi sqrt(se2) psi-crit*sqrt(se2) psi+crit*sqrt(se2)]

% Profile Shape
M=[-3 -1 1 3; 1 -1 -1 1; -1 3 -3 1]';
M=M*inv(sqrtm(M'*M))
psi=M'*ybar;
T2=N*psi'*inv(M'*S*M)*psi
d=3;
F=T2*(N-d)/(d*(N-1))
crit=finv(.95,d,N-d)
pval = 1-fcdf(F,d,N-d)
se2=diag(M'*S*M/N);
crit=sqrt(d*(N-1)/(N-d)*finv(.95,d,N-d))
[psi sqrt(se2) psi-crit*sqrt(se2) psi+crit*sqrt(se2)]

% Combined Level/Shape 
M=[eye(4) ones(4,1)/4 M]
d=4;
psi=M'*ybar;
se2=diag(M'*S*M/N);
crit=sqrt(d*(N-1)/(N-d)*finv(.95,d,N-d))
[psi sqrt(se2) psi-crit*sqrt(se2) psi+crit*sqrt(se2)]