disp('Running MATLAB script file lmtzcf54.m') ;
%
%    Makes Figure 5.4, for the Locantore, Marron, Tripoli,
%    Zhang, Cohen paper.
%
%    Copied from gpcd2t.m and gpcd2.m  (itest == 4)


%  Example from classes\322\s322eg12.m  (with outlier added)
d = 10 ;
n = 50 ;
xgrid = (.5:1:d)' ;
data = (xgrid - 6).^2 ;
  randn('seed',88769874) ;
  eps1 = 4 * randn(1,n) ;
  eps2 = .5 * randn(1,n) ;
  eps3 = 1 * randn(d,n) ;
data = vec2mat(data,n) + vec2mat(eps1,d) + ...
               vec2mat(eps2,d) .* vec2mat(xgrid-d/2,n) + eps3 ;
data = [data, 15 * (sin(pi * xgrid) + 1)] ;



%  set preliminary stuff
%
rand('seed',20983743) ;



%  Get statistics and do eigenvectors
%
d = size(data,1) ;
         %  dimension of each data curve
n = size(data,2) ;
         %  number of data curves
nev = min([d,n,3]) ;
         %  number of eigenvectors to display


%  get mean and subtract
%
vmean = mean(data')' ;
ssm = n * sum(vmean.^2) ;
mresid = data - vec2mat(vmean,n) ;
ssmr = sum(sum(mresid .^ 2)) ;

vrad = sqrt(sum(mresid .^2)) ;
        %  Transpose, since want "coordinates as variables"
        %  vector of radii of each curve
sphereresid = mresid' ./ vec2mat(vrad',d) ;
        %  make each curve have "length" one
mcov = cov(sphereresid) ;
        %  Get covariance matrix of "sphered" data


[meigvec,meigval] = eig(mcov) ;
veigval = (diag(meigval)) ;
          %  Get eigenvalues and eigenvectors
[veigval,vind] = sort(veigval) ;
meigvec = meigvec(:,vind) ;
          %  Sort to avoid possible ordering problems
vev = flipud(veigval) ;
          %  put largest one first

vev = vev / sum(vev) ;
          %  put on scale of R^2
cev = cumsum(vev) ;
          %  get cumulative sum


ev1 = meigvec(:,d) ;
          %  eigenvector with largest eigenvalue
pc1 = mresid' * ev1 ;
          %  vector of projections of data onto ev1
mpc1 = ev1 * pc1' ;
          %  matrix for display
ss1 = sum(sum(mpc1.^2)) ;
ss1r = ssmr - ss1 ;
mresid1 = (mresid - mpc1) ;
vr2 = ss1 / ssmr ;
vflast = ss1r / ssmr ;

ev2 = meigvec(:,d-1) ;
          %  eigenvector with 2nd largest eigenvalue
pc2 = mresid1' * ev2 ;
          %  vector of projections of data onto ev2
mpc2 = ev2 * pc2' ;
          %  matrix for display
ss2 = sum(sum(mpc2.^2)) ;
ss2r = ss1r - ss2 ;
mresid2 = (mresid1 - mpc2) ;
vr2 = [vr2, ss2 / ssmr] ;
vflast = [vflast, ss2r / ss1r] ;


if nev > 2 ;

  ev3 = meigvec(:,d-2) ;
          %  eigenvector with 3rd largest eigenvalue
  pc3 = mresid2' * ev3 ;
          %  vector of projections of data onto ev3
  mpc3 = ev3 * pc3' ;
          %  matrix for display
  ss3 = sum(sum(mpc3.^2)) ;
  ss3r = ss2r - ss3 ;
  mresid3 = (mresid2 - mpc3) ;
  vr2 = [vr2, ss3 / ssmr] ;
  vflast = [vflast, ss3r / ss2r] ;

  if nev > 3 ;

    ev4 = meigvec(:,d-3) ;
          %  eigenvector with 3rd largest eigenvalue
    pc4 = mresid3' * ev4 ;
          %  vector of projections of data onto ev4
    mpc4 = ev4 * pc4' ;
          %  matrix for display
    ss4 = sum(sum(mpc4.^2)) ;
    ss4r = ss3r - ss4 ;
    mresid4 = (mresid3 - mpc4) ;
    vr2 = [vr2, ss4 / ssmr] ;
    vflast = [vflast, ss4r / ss3r] ;

  end ;

end ;

cr2 = cumsum(vr2) ;


%  make output plot
%
npx = 3 ;
npy = 4 ;
xgrid = (1:d)' ;
min1 = min(min(data)) ;
max1 = max(max(data)) ;
range1 = max1 - min1 ;
min1 = min1 - 0.05 * range1 ;
max1 = max1 + 0.05 * range1 ;
min2 = min(min(mresid)) ;
max2 = max(max(mresid)) ;
range2 = max2 - min2 ;
min2 = min([min2; -max2]) ;
max2 = max([max2; -min2]) ;
min2 = min2 - 0.05 * range2 ;
max2 = max2 + 0.05 * range2 ;


subplot(npy,npx,1) ;
  plot(xgrid,data) ;
    title('Raw Data') ;
    axis([1,d,min1,max1]) ;
    set(gca,'XTickLabel',[]) ;
    set(gca,'FontSize',8) ;
    set(get(gca,'Title'),'FontSize',8) ;


subplot(npy,npx,2) ;
  plot(xgrid,vmean) ;
    title('Mean') ;
    axis([1,d,min1,max1]) ;
    set(gca,'XTickLabel',[]) ;
    set(gca,'FontSize',8) ;
    set(get(gca,'Title'),'FontSize',8) ;


subplot(npy,npx,3) ;
  plot(xgrid,mresid) ;
    title('Mean Resid.') ;
    axis([1,d,min2,max2]) ;
    set(gca,'XTickLabel',[]) ;
    set(gca,'FontSize',8) ;
    set(get(gca,'Title'),'FontSize',8) ;



%  loop through principal components
%
for iev = 1:nev ;

  istr = num2str(iev) ;

  subplot(npy,npx,3 * iev + 1) ;
  eval(['plot(xgrid,mpc' istr ')']) ;
    title(['PC' num2str(iev)]) ;
    axis([1,d,min2,max2]) ;
    set(gca,'XTickLabel',[]) ;
      set(gca,'FontSize',8) ;
    set(get(gca,'Title'),'FontSize',8) ;


  %  do residual plot
  subplot(npy,npx,3 * iev + 3) ;
  eval(['plot(xgrid,mresid' istr ')']) ;
    title(['resid, PC' num2str(iev)]) ;
    axis([1,d,min2,max2]) ;
    set(gca,'XTickLabel',[]) ;
    set(gca,'FontSize',8) ;
    set(get(gca,'Title'),'FontSize',8) ;


  %  do mean +- plot
  subplot(npy,npx,3 * iev + 2) ;
    plot(xgrid,vmean) ;
      title(['Mean +- PC' istr]) ;
      axis([1,d,min1,max1]) ;
      hold on ;
        eval(['plot(xgrid,vmean + max(pc' istr ') * ev' istr ',''--'') ;']) ;
        eval(['plot(xgrid,vmean + min(pc' istr ') * ev' istr ','':'') ;']) ;
      hold off ;
      set(gca,'XTickLabel',[]) ;
      set(gca,'FontSize',8) ;
      set(get(gca,'Title'),'FontSize',8) ;


end ;





  orient tall ;
  set(gcf,'PaperPosition',[1.35, 3.85, 5.8, 5.8]) ;
        %  1.35 inches on left and right, 5.8 interior, 8.5 total
        %  3.85 from bottom, 5.8 interior, 11 total
print -deps \matlab\steve\ps\lmtzcf54.eps ;


