disp(['Running MATLAB script s322eg12.m']) ;
%
%    EXAMPLE 12 FOR STATISTICS 322
%    POSTED 10/31/96
%    EXPLORES CURVES AS DATA, SHOWING SOME TOY EXAMPLES
%    

format compact ;


ieg = 2 ;     %  1,2


d = 10 ;        %  Dimension of data vectors
xgrid = (.5:1:d)' ;
left = 0 ;
right = d ;

n = 40 ;       %  Number of observations

if ieg <= 2 ;
  titstr = 'E.g. 1' ;
  mdata = (xgrid - 6).^2 ;
    randn('seed',88769874) ;
    eps1 = 4 * randn(1,n) ;
    eps2 = .5 * randn(1,n) ;
    eps3 = 1 * randn(d,n) ;
  mdata = vec2mat(mdata,n) + vec2mat(eps1,d) + ...
                 vec2mat(eps2,d) .* vec2mat(xgrid-d/2,n) + eps3 ;
  bottom1 = -10 ;   top1 = 40 ;
  bottom2 = -12 ;   top2 = 12 ;

  if ieg == 2 ;
    titstr = 'E.g. 2' ;
    mdata = [mdata (30 * sin(pi * xgrid))] ;
    mdata = [mdata (20 * cos(pi * xgrid / 2))] ;
    n = n + 2 ;
    bottom1 = -30 ;   top1 = 40 ;
    bottom2 = -30 ;   top2 = 30 ;
  end ;
end ;


%  Get statistics and do eigenvectors
sst = sum(sum(mdata.^2)) ;
vmean = mean(mdata') ;
vsd = std(mdata') ;
ssm = n * sum(vmean.^2) ;
ssmr = sst - ssm ;
mresid = mdata - vec2mat(vmean',n) ;

mcov = cov(mresid') ;
          %  Get covariance matrix, transpose, since want "coordinates
          %                      as variables"
[meigvec,meigval] = eig(mcov) ;
veigval = (diag(meigval)) ;
          %  Get eigenvalues and eigenvectors
[veigval,vind] = sort(veigval) ;
meigvec = meigvec(:,vind) ;
          %  Sort to avoid possible ordering problems
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) ;

ev2 = meigvec(:,d-1) ;
          %  eigenvector with 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) ;

ev3 = meigvec(:,d-2) ;
          %  eigenvector with largest eigenvalue
pc3 = mresid1' * 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) ;

ev4 = meigvec(:,d-3) ;
          %  eigenvector with largest eigenvalue
pc4 = mresid1' * 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) ;



%  Make plots
subplot(3,4,1) ;
  plot(xgrid,mdata) ;
    title([titstr ', Raw Data']) ;
    axis([0 d bottom1 top1]) ;
      ssstr = ['SS,T = ' num2str(sst)] ;
      textx = 0.5 + .1 * (d - 0) ;
      texty = bottom1 + .9 * (top1 - bottom1) ;
    text(textx,texty,ssstr) ;  

  
subplot(3,4,3) ;
  plot(xgrid,vmean,'-',...
       xgrid,vmean + 2 * vsd,':',...
       xgrid,vmean - 2 * vsd,':') ;
    title([titstr ', Mean +- 2 sd']) ;
    axis([0 d bottom1 top1]) ;
      ssstr = ['% of SS,T = ' num2str(100 * ssm / sst)] ;
      textx = 0.5 + .1 * (d - 0) ;
      texty = bottom1 + .9 * (top1 - bottom1) ;
    text(textx,texty,ssstr) ;  

subplot(3,4,4) ;
  plot(xgrid,mresid) ;
    title([titstr ', Res. from Mean']) ;
    axis([0 d bottom2 top2]) ;
      ssstr = ['% of SS,T = ' num2str(100 * ssmr / sst)] ;
      textx = 0.5 + .1 * (d - 0) ;
      texty = bottom1 + .9 * (top2 - bottom2) ;
    text(textx,texty,ssstr) ;  
  
subplot(3,4,5) ;
  plot(xgrid,mpc1) ;
    title([titstr ', PC1']) ;
    axis([0 d bottom2 top2]) ;
      ssstr = ['% of SS,MR = ' num2str(100 * ss1 / ssmr)] ;
      textx = 0.5 + .1 * (d - 0) ;
      texty = bottom2 + .9 * (top2 - bottom2) ;
    text(textx,texty,ssstr) ;  
    
subplot(3,4,9) ;
  plot(xgrid,mresid1) ;
    title([titstr ', Res. from PC1']) ;
    axis([0 d bottom2 top2]) ;
      ssstr = ['% of SS,MR = ' num2str(100 * ss1r / ssmr)] ;
      textx = 0.5 + .1 * (d - 0) ;
      texty = bottom2 + .9 * (top2 - bottom2) ;
    text(textx,texty,ssstr) ;  

  
subplot(3,4,6) ;
  plot(xgrid,mpc2) ;
    title([titstr ', PC2']) ;
    axis([0 d bottom2 top2]) ;
      ssstr = ['% of SS,1R = ' num2str(100 * ss2 / ss1r)] ;
      textx = 0.5 + .1 * (d - 0) ;
      texty = bottom2 + .9 * (top2 - bottom2) ;
    text(textx,texty,ssstr) ;  
    
subplot(3,4,10) ;
  plot(xgrid,mresid2) ;
    title([titstr ', Res. from PC2']) ;
    axis([0 d bottom2 top2]) ;
      ssstr = ['% of SS,1R = ' num2str(100 * ss2r / ss1r)] ;
      textx = 0.5 + .1 * (d - 0) ;
      texty = bottom2 + .9 * (top2 - bottom2) ;
    text(textx,texty,ssstr) ;  

  
subplot(3,4,7) ;
  plot(xgrid,mpc3) ;
    title([titstr ', PC3']) ;
    axis([0 d bottom2 top2]) ;
      ssstr = ['% of SS,2R = ' num2str(100 * ss3 / ss2r)] ;
      textx = 0.5 + .1 * (d - 0) ;
      texty = bottom2 + .9 * (top2 - bottom2) ;
    text(textx,texty,ssstr) ;  
    
subplot(3,4,11) ;
  plot(xgrid,mresid3) ;
    title([titstr ', Res. from PC3']) ;
    axis([0 d bottom2 top2]) ;
      ssstr = ['% of SS,2R = ' num2str(100 * ss3r / ss2r)] ;
      textx = 0.5 + .1 * (d - 0) ;
      texty = bottom2 + .9 * (top2 - bottom2) ;
    text(textx,texty,ssstr) ;  
  
  
subplot(3,4,8) ;
  plot(xgrid,mpc4) ;
    title([titstr ', PC4']) ;
    axis([0 d bottom2 top2]) ;
      ssstr = ['% of SS,3R = ' num2str(100 * ss4 / ss3r)] ;
      textx = 0.5 + .1 * (d - 0) ;
      texty = bottom2 + .9 * (top2 - bottom2) ;
    text(textx,texty,ssstr) ;  
    
subplot(3,4,12) ;
  plot(xgrid,mresid4) ;
    title([titstr ', Res. from PC4']) ;
    axis([0 d bottom2 top2]) ;
      ssstr = ['% of SS,3R = ' num2str(100 * ss4r / ss3r)] ;
      textx = 0.5 + .1 * (d - 0) ;
      texty = bottom2 + .9 * (top2 - bottom2) ;
    text(textx,texty,ssstr) ;  

  orient landscape ;
eval(['print -dps \matlab\steve\ps\eg12p' num2str(ieg) '.ps']) ;




