disp(['Running MATLAB script s322eg11.m']) ;
%
%    EXAMPLE 11 FOR STATISTICS 322
%    POSTED 10/29/96
%    EXPLORES PRINCIPAL COMPONENT IDEAS,
%              FOR 3 DIMENSIONAL TOY EXAMPLES
%    

format compact ;

ieg = 1 ;     %  1 - SPHERICAL NORMAL
              %  2 - PARALLEL TO AXIS
              %  3 - ODD ANGLE

iview = 1 ;   %  1 - PRINTED OUTPUT OF SUMMARY STATISTICS
              %  2 - DRAFTSMAN'S PLOT
              %  3 - FIRST ATTEMPT AT SPINNING POINT CLOUD
              %  4 - BETTER VERSION OF SPINNING POINT CLOUD

autop = 1 ;     %  0 FOR NO AUTOPRINT


n = 50 ;       %  Number of observations
d = 3 ;

%  SET EXAMPLE PARAMETERS
if ieg == 1 ; 
  titstr = 'Spherical Normal' ;
  mu = zeros(3,1) ;
  msig = eye(3) ;
  randn('seed',44392887) ;
  cmax = 2.5 ;
elseif ieg == 2 ; 
  titstr = 'Parallel to axes' ;
  mu = [0; 1; 2] ;
  msig = [16 0 0; 0 4 0; 0 0 1] ;
  cmax = 10 ;
  randn('seed',98834531) ;
elseif ieg == 3 ; 
  titstr = 'Odd Angle' ;
  mu = [3; 2; 1] ;
  msig = [1 .1 .9; .1 1 .5; .9 .5 1] ;
  cmax = 5 ;
  randn('seed',11643198) ;
end ;
vax2 = [-cmax cmax -cmax cmax] ;
          %  FOR USE AS 2-d AXIS
vax3 = [vax2 -cmax cmax] ;
          %  FOR USE AS 3-d AXIS


%  GENERATE MULTIVARIATE NORMAL DATA
mz = randn(3,n) ;
mdata = sqrtm(msig) * mz + vec2mat(mu,n) ;
          %  CAREFUL ABOUT MATLAB'S sqrt VS. sqrtm


%  GET SUMMARY STATISTICS, AND DO EIGENVECTOR-VALUE DECOMPOSITION
sst = sum(sum(mdata.^2)) ;
          %  TOTAL SUM OF SQUARES
vmean = mean(mdata') ;
vsd = std(mdata') ;
ssm = n * sum(vmean.^2) ;
          %  SUM OF SQUARES EXPLAINED BY THE MEAN
ssmr = sst - ssm ;
mresid = mdata - vec2mat(vmean',n) ;
          %  RESIDUALS AFTER SUBTRACTION OF THE MEAN FROM THE DATA

mcov = cov(mresid') ;
          %  GET COVARIANCE MATRIX
          %  USE TRANSPOSE, SINCE DATA VECTORS ARE COLUMN VECTORS
          %  WHILE cov ASSUMES THEY ARE IN ROWS
[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
ss1 = sum(pc1.^2) ;
          %  SUM OF SQUARES EXPLAINED BY THE FIRST PRINCIPAL COMPONENT
ss1r = ssmr - ss1 ;

ev2 = meigvec(:,d-1) ;
          %  EIGENVECTOR WITH SECOND LARGEST EIGENVALUE
pc2 = mresid' * ev2 ;
          %  VECTOR OF PROJECTIONS OF DATA ONTO EV1
ss2 = sum(pc2.^2) ;
          %  SUM OF SQUARES EXPLAINED BY THE 2ND PRINCIPAL COMPONENT
ss2r = ss1r - ss2 ;

ev3 = meigvec(:,1) ;
          %  EIGENVECTOR WITH SMALLEST EIGENVALUE
pc3 = mresid' * ev3 ;
          %  VECTOR OF PROJECTIONS OF DATA ONTO EV1
ss3 = sum(pc3.^2) ;
          %  SUM OF SQUARES EXPLAINED BY THE 2ND PRINCIPAL COMPONENT
ss3r = ss2r - ss3 ;


if iview == 1 ;   %  THEN PRINT SUMMARY STATISTICS

  if ieg == 1 ;
    fid = fopen('s322eg11.out','wt') ;
      titlstr1 = 'Output from the MATLAB Script s322eg11.m,   ' ;
      titlstr1 = [titlstr1,date] ;
    cntbytes = fprintf(fid,'%1s\n',titlstr1) ;
  elseif ieg >= 2 ;
    fid = fopen('s322eg11.out','at') ;
  end ;

    titlstr2 = ['For the ' titstr ' Data'] ;
  cntbytes = fprintf(fid,'\n\n\n%1s\n\n',titlstr2) ;

    forstr = 'With Mean Vector:      %5.2f   %5.2f   %5.2f \n\n' ;
  cntbytes = fprintf(fid,forstr,mu) ;
    str = 'And Covar Matrix:' ;
  cntbytes = fprintf(fid,'%1s\n',str) ;
    forstr = '                       %5.2f   %5.2f   %5.2f\n' ;
  cntbytes = fprintf(fid,forstr,msig') ;

    forstr = '\nTotal Sum of Squares is:        %10.4f' ;
  cntbytes = fprintf(fid,forstr,sst) ;

    str = 'For:      S of S       SS Rem      % of Tot.  % of T - M' ;
    str = [str '  % of last'] ;
  cntbytes = fprintf(fid,'\n\n%1s\n',str) ;
    forstr = 'Mean   %10.4f   %10.4f      %6.2f      %6.2f     %6.2f\n' ;
    vout = [ssm ssmr (100*ssm/sst) 100 (100*ssm/sst)] ;
  cntbytes = fprintf(fid,forstr,vout) ;
    forstr = ' pc1   %10.4f   %10.4f      %6.2f      %6.2f     %6.2f\n' ;
    vout = [ss1 ss1r (100*ss1/sst) (100*ss1/ssmr) (100*ss1/ssmr)] ;
  cntbytes = fprintf(fid,forstr,vout) ;
    forstr = ' pc2   %10.4f   %10.4f      %6.2f      %6.2f     %6.2f\n' ;
    vout = [ss2 ss2r (100*ss2/sst) (100*ss2/ssmr) (100*ss2/ss1r)] ;
  cntbytes = fprintf(fid,forstr,vout) ;
    forstr = ' pc3   %10.4f   %10.4f      %6.2f      %6.2f     %6.2f\n' ;
    vout = [ss3 ss3r (100*ss3/sst) (100*ss3/ssmr) (100*ss3/ss2r)] ;
  cntbytes = fprintf(fid,forstr,vout) ;

  fclose('all') ;


elseif iview == 2 ;   %  THEN CONSTRUCT DRAFTSMAN'S PLOTS
  subplot(2,2,1) ;  
    plot(mdata(1,:),mdata(2,:),'+') ;
      axis(vax2) ;
      title(titstr) ;
      xlabel('x1') ;
      ylabel('x2') ;

  subplot(2,2,2) ;
    plot(mdata(3,:),mdata(2,:),'+') ;
      axis(vax2) ;
      title(titstr) ;
      xlabel('x3') ;
      ylabel('x2') ;
                  
  subplot(2,2,3) ;
    plot(mdata(1,:),mdata(3,:),'+') ;
      axis(vax2) ;
      title(titstr) ;
      xlabel('x1') ;
      ylabel('x3') ;
       
  if autop == 1 ;     %  then write plot to postscript file
    eval(['print -dps \matlab\steve\ps\eg11' num2str(ieg) 'a.ps']) ;
           %  CAN DO SOMETHING LIKE THIS TO DIRECT FILES TO
           %  ANOTHER DIRECTORY
  end ;


elseif iview == 3 ;   %  FIRST ATTEMPT AT SPINNING 3-d PLOT
  subplot(1,1,1) ;
          %  THIS GETS RID OF THE 2x2 ARRAY OF SUBPLOTS THAT WERE
          %  SET BEFORE, SO THIS NOW OCCUPIES THE WHOLE SCREEN.

  hfig = figure(1) ;
          %  THIS IS A "HANDLE" TO THE LOW LEVEL GRAPHICS OBJECT,
          %  WHICH IS THE "FIGURE" (ESSENTIALLY THE WINDOW ON THE 
          %  SCREEN.
  nrotate = 10 ;
  nframe = 50 ;
          %  SET THIS TO SOMETHING SMALLER IF THE MOVIE GOES TOO SLOW

  disp('Storing Movie') ;
  M = moviein(nframe,hfig) ;
          %  THIS CREATES A BIG MATRIX IN RAM WHICH WILL BE FILLED
          %  WITH THE DATA FOR SHOWING THE MOVIE.
  for iframe = 1:nframe ;
          %  LOOP THROUGH MOVIE FRAMES
    plot3(mdata(1,:),mdata(2,:),mdata(3,:),'.') ;
          %  MAKE A "3-D SCATTERPLOT"
      title(titstr) ;
      axis([-5 5 -5 5 -5 5]) ;
          %  THIS IS IMPORTANT, BECAUSE OTHERWISE THE VIEW CHANGES
          %  TOO MUCH EACH TIME (SINCE THE AXES ARE RECHOSEN FOR
          %  EACH FRAME).  SEE THIS BY COMMENTING OUT THIS LINE.
      axis('off') ;
          %  THIS IS USEFUL, BECAUSE OTHERWISE THE STANDARD
          %  AXES CHANGE IN A DISTRACTING WAY
      view(360*iframe/nframe,30) ;
          %  AZIMUTHAL ANGLE, ELEVATION ANGLE (BOTH IN DEGRESS)
%  THIS MOVIE IS NOT VERY GOOD, BECAUSE IT "MOVES AROUND IN AN
%  UNPLEASANT WAY", SINCE FOR EACH FRAME THE BOX CHOSEN FILLS
%  THE WINDOW AS MUCH AS IT CAN.  THIS MEANS THE DATA LOOK CLOSER
%  WHEN THE BOX IS "END ON", AND FURTHER WHEN WE ARE LOOKING AT 
%  THE BOX "FROM A DIAGONAL DIRECTION".
%  HERE WAS AN ATTEMPT TO FIX THIS PROBLEM, WHICH DIDN'T WORK EITHER
%        view(viewmtx(360*iframe/nframe,30,25,[0,0,0])) ;
    M(:,iframe) = getframe(hfig) ;
          %  THIS CREATES ONE MOVIE SCREEN FROM WHAT IS ON THE SCREEN
          %  IN THE WINDOW LABELLED "Figure 1" (INDEXED BY THE HANDLE
          %  figh) RIGHT NOW.
  end ;

  disp('Finished storing movie, hit any key to run') ;
  pause ;
          %  TO SEE THIS MESSAGE, AND PROPERLY RESPOND, YOU NEED TO 
          %  HAVE BOTH THE COMMAND WINDOW, AND THE FIGURE WINDOW
          %  CLEARLY VISIBLE

  movie(hfig,M,nrotate) ;
          %  THIS PLAYS THE MOVIE, RUNNING THROUGH nrotate TIMES


elseif iview == 4 ;   %  SECOND ATTEMPT AT SPINNING 3-d PLOT
%  TO SOLVE THE ABOVE PROBLEM, I LOOKED AT SOME OF THE DEMOS, WHICH
%  SHOWED SPINNING OBJECTS, AND THEN LOOKED AT THE SOURCE CODE FOR
%  THESE. THE ONE CALLED "crulspin.m" WAS ESPECIALLY USEFUL IN
%  SHOWING HOW TO MAKE THINGS ROTATE WITHOUT THE "BOX-CORNER" 
%  EFFECT WHICH WAS A PROBLEM ABOVE.
  subplot(1,1,1) ;

  hfig = figure(1) ;
  hobj = plot3(mdata(1,:),mdata(2,:),mdata(3,:),'.') ;
          %  THIS GIVES A "HANDLE TO THE GRAPHICS OBJECT, WHICH IS THE 
          %  POINT CLOUD".  THEN WILL "SPIN THE OBJECT".
          %  THE FUNDAMENTAL DIFFERENCE BETWEEN THIS AND iview = 2
          %  IS THAT NOW WE WILL "SPIN THE OBJECT WITHIN THE AXES"
          %  AS OPPOSED TO "SPINNING THE AXES" (WHICH THEN HAD TO
          %  BE REFIT TO THE AXES EACH TIME).
    title(titstr) ;
    axis(vax3) ;

  nrotate = 10 ;
  nframe = 50 ;
          %  SET THIS TO SOMETHING SMALLER IF THE MOVIE GOES TOO SLOW

  M = moviein(nframe,hfig) ;
  for iframe = 1:nframe ;
      azel = [0,90] ;
          %  AZIMUTH AND ELEVATION, OF AXIS TO ROTATE AROUND
      alpha = 360 / nframe ;
          %  ANGLE TO ROTATE FOR EACH FRAME
      origin = [0 0 0] ;
          %  CENTERPOINT FOR ROTATIONS
    rotate(hobj,azel,alpha,origin);
          %  THIS TAKES THE EXISTING OBJECT, AND REPLACES IT BY
          %  A VERSION THAT HAS BEEN "ROTATED BY ANGLE ALPHA,
          %  AROUND THE GIVEN AXIS, CENTERED AT THE GIVEN ORIGIN"
    hobj=plot3(get(hobj,'XData'),get(hobj,'YData'),...
                                         get(hobj,'ZData'),'.');
          %  get COMMANDS EXTRACT THE NEEDED PIECES TO PLOT.
      title(titstr) ;
      axis(vax3) ;
          %  THIS KEEPS THE AXES FROM CHANGING IN AN UNPLEASANT WAY
    M(:,iframe) = getframe(hfig) ;
  end ;

  disp('Finished storing movie, hit any key to run') ;
  pause ;

  movie(hfig,M,nrotate) ;


%  IT WOULD HAVE BEEN NICE TO ADD THE PRINCIPLE COMPONENT AXES AS 
%  LINES (E.G. USING hold), BUT I COULD NOT MAKE IT WORK IN THE 
%  TIME ALLOTTED.  IF I HAD MORE TIME, I WOULD INVESTIGATE THE 
%  "PARENT CHILD" STURCTURE OF THE GRAPHICS MORE CAREFULLY.

end ;

