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

format compact ;


ieg = 1 ;     %  1,2

n = 50 ;       %  Number of observations
d = 3 ;
xgrid = 1:d ;

if ieg == 1 ; 
  titstr = '1:  Dog Legs' ;
  mu = [1; 0; .5] ;
  msig = [.3 0 .29; 0 .05 0; .29 0 .3] ;
  cmax = 2 ;
  randn('seed',23729877) ;
elseif ieg == 2 ; 
  titstr = '2:  Fans' ;
  mu = [0; .5; 1] ;
  msig = [.01 .01 .01; .01 .25 .25; .01 .25 .5] ;
  cmax = 2.5 ;
  randn('seed',97987373) ;
end ;


%  Generate multivariate normal data
mz = randn(3,n) ;
mdata = sqrtm(msig) * mz + vec2mat(mu,n) ;
          %  Careful about sqrt vs. sqrtm



%  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) ;


%  Make draftsman's plot
subplot(2,2,1) ;
  plot(mdata(1,:),mdata(2,:),'+') ;
    axis([-cmax cmax -cmax cmax]) ;
    title(titstr) ;
    xlabel('x1') ;
    ylabel('x2') ;

subplot(2,2,2) ;
  plot(mdata(3,:),mdata(2,:),'+') ;
    axis([-cmax cmax -cmax cmax]) ;
    title(titstr) ;
    xlabel('x3') ;
    ylabel('x2') ;
                  
subplot(2,2,3) ;
  plot(mdata(1,:),mdata(3,:),'+') ;
    axis([-cmax cmax -cmax cmax]) ;
    title(titstr) ;
    xlabel('x1') ;
    ylabel('x3') ;
       
subplot(2,2,4) ;
  plot((1:3)',mdata,'-') ;
    axis([.5 3.5 -cmax cmax]) ;
    title(titstr) ;
    xlabel('x1') ;
    ylabel('x3') ;

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



disp('Hit any key to make Spectral Decomposition Plots') ;
pause ;

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

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

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

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




