%
%  Computation of optimal contrast coefficients for Profile Analysis
%
%  See Rencher 2002 page 201 for the data set
%
load pig.dat;
grp=pig(:,1);
Y=pig(:,2:end);
Y=(Y/500).^5;

d=size(Y,2);
x=dummyvar(grp);
n=length(grp);
X=[ones(n,1) x];
Px=ppo(X);
P1=ppo(ones(n,1));
In=speye(n);
B=pinv(X'*X)*X'*Y;
E_save=Y'*(In-Px)*Y;

Ybar=([1 1 0 0
      1 0 1 0
      1 0 0 1]*B)'
week=[1 3 4 5 6 7]';
plot(week,Ybar(:,1),'-',week,Ybar(:,2),'--',week,Ybar(:,3),'-.')
set(gca,'FontSize',16)
xlabel('Week')
ylabel('Weight')
set(gca,'LineWidth',2)
legend('Zero Vit E','Low Vit E', 'High Vit E','Location','NorthWest')
title('Growth of Guinea Pigs')
hold on
plot(week,Ybar(:,1),'*',week,Ybar(:,2),'*',week,Ybar(:,3),'*','MarkerSize',10)
hold off
print -deps profile3m.eps


%
% Optimal Coefficient For Group Effect
%
disp(['Group Coefficient'])
M=ones(d,1);
E=M'*E_save*M;
L=[0 1 -1 0
   0 1 0 -1];
c_grp=L'*inv(L*pinv(X'*X)*L')*L*B*M;
c_grp=c_grp/sqrt(c_grp'*c_grp);
disp(c_grp)

%
% Optimal Coefficient For Time Effect
%
disp(['Time Coefficient'])
M=[ 1  1  1  1  1
   -1  0  0  0  0
    0 -1  0  0  0
    0  0 -1  0  0
    0  0  0 -1  0
    0  0  0  0 -1];
E=M'*E_save*M;
L=[3 1 1 1]/3;
c_time=M*inv(E)*M'*B'*L';
c_time=c_time/sqrt(c_time'*c_time);
disp(c_time);
%
% Optimal Coefficients For Time x Group Effect
%
disp(['Group x Time Coefficients'])

L=[0 1 -1 0; 0 1 0 -1];
H=(L*B*M)'*inv(L*pinv(X'*X)*L')*L*B*M;

[Evectm,Lam]=eig(inv(E)*H);
lam=diag(Lam);
[tmp,jj]=sort(-abs(lam));
lam=abs(lam(jj));
Evectm=Evectm(:,jj);
q=length(E);
f=size(L,1);
s=min(q,f);
lam=lam(1:s);
disp(['Roots'])
disp([lam cumsum(lam)/sum(lam)])
c_grp_gxt = L'*inv(L*pinv(X'*X)*L')*L*B*M*Evectm(:,1:s);
c_time_gxt = M*Evectm(:,1:s);
c_time_gxt = c_time_gxt*inv(sqrtm(c_time_gxt'*c_time_gxt));
c_grp_gxt = c_grp_gxt*inv(sqrtm(c_grp_gxt'*c_grp_gxt));
disp(['Coefficients'])
disp(c_time_gxt)
disp(c_grp_gxt)