disp(['Running MATLAB script s322hw9.m']) ;
%
%    SOLUTION TO HOMEWORK 9 FOR STATISTICS 322
%    POSTED 11/7/96
%    EXPLORES ROBUSTNESS IN PRINCIPAL COMPONENTS
%    

format compact ;

iquest = 3 ;      %  1,2,3

iview = 2 ;       %  1,2 - two different movies for answer to Quest. 3

%  GENERATE DATA SET
if iquest <= 2 ;
  titstr = 'Graph for Question 1' ;
    n = 100 ;   
    d = 2 ;
    mu = [-6; -4] ;
    msig = [2 -1.4;-1.4 1] ;
    randn('seed',82038650) ;
  mz = randn(d,n) ;
  mdata = sqrtm(msig) * mz + vec2mat(mu,n) ;
          %  CAREFUL ABOUT MATLAB'S sqrt VS. sqrtm
  if iquest == 2 ; 
    titstr = 'Graph for Question 2' ;
    mdata = [mdata [19; 14]] ;
    n = n + 1 ;
  end ;
elseif iquest >= 3 ;
    n = 60 ;   
    d = 3 ;
    mu1 = [-5; -2; 0] ;
    mu2 = [0; 4; 0] ;
    mu3 = [5; -3; 0] ;
    msig = [.5 0 0; 0 .5 0; 0 0 3] ;
          %  COVARIANCE OF A DISTRIBUTION WHICH IS "STRETCHED 
          %          IN Z DIRECTION"
    randn('seed',12234602) ;
  mz = randn(d,n) ;
  mdata = sqrtm(msig) * mz ;
          %  PSEUDO-DATA CENTERED AT 0 WITH THIS COVARIANCE ;
    rand('seed',16796234) ;
  vind = ceil(3 * rand(n,1)) ;
          %  VECTOR OF RANDOM 1,2,3'S 
  mmu = vec2mat(mu1,n) ;
          %  START WITH ALL mu1
    vind2 = find(vind == 2)' ;
  mmu(:,vind2) = vec2mat(mu2,size(vind2,2)) ;
          %  REPLACE APPROPRIATE COLS WITH mu2
    vind3 = find(vind == 3)' ;
  mmu(:,vind3) = vec2mat(mu3,size(vind3,2)) ;
          %  REPLACE APPROPRIATE COLS WITH mu3
  mdata = mdata + mmu ;
end ;


%  DO EIGENVECTOR-VALUE DECOMPOSITION
vmean = mean(mdata')' ;
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
ev2 = meigvec(:,d-1) ;
          %  EIGENVECTOR WITH SECOND LARGEST EIGENVALUE
if iquest >= 3 ;
  ev3 = meigvec(:,1) ;
          %  EIGENVECTOR WITH THIRD LARGEST EIGENVALUE

  %   AND ALSO GET SUBPOPULATION EIGENVECTORS
  [meigvecs,meigvals] = eig(msig) ;
  veigvals = (diag(meigvals)) ;
  [veigvals,vinds] = sort(veigvals) ;
  meigvecs = meigvecs(:,vinds) ;
          %  SORT TO AVOID POSSIBLE ORDERING PROBLEMS
  ev1s = meigvecs(:,3) ;
  ev2s = meigvecs(:,2) ;
  ev3s = meigvecs(:,1) ;
end ;


%  DO GRAPHICS
if iquest <= 2 ;
  plot(mdata(1,:),mdata(2,:),'x') ;
    title(titstr) ;
    hold on ;
        sev1 = 3 * sqrt(veigval(d)) * ev1 ;
          %  SCALED VERSION OF ev1 TO USE FOR AXIS, ADJUSTED SO LENGTH
          %  SHOWS "3 STANDARD DEVIATIONS"
        mlev1 = [(vmean - sev1) (vmean + sev1)] ;
      plot(mlev1(1,:),mlev1(2,:),'-r') ;
          %  PUT FIRST EIGENVECTOR IN SOLID RED
        sev2 = 3 * sqrt(veigval(d-1)) * ev2 ;
        mlev2 = [(vmean - sev2) (vmean + sev2)] ;
      plot(mlev2(1,:),mlev2(2,:),':b') ;
          %  PUT SECOND EIGENVECTOR IN DOTTED BLUE
    hold off ;
    axis('equal') ;
          %  MAKE ASPECT RATIO 1, SO PERPENDICULAR VECTORS LOOK THAT WAY
    legend('Data','EV1','EV2') ;
    if iquest == 2 ;
      text(-12,13,'PC directions mostly driven by outlier,') ;
      text(-12,11,'i.e. they "ignore the bulk of the data"') ;
    end ;
    orient landscape ;

    eval(['print -dps \matlab\steve\ps\hw9p' num2str(iquest) '.ps']) ;

elseif iquest == 3 ;

  if iview == 1 ;
    hfig = figure(1) ;
    subplot(2,2,1) ;
      plot3(mdata(1,:),mdata(2,:),mdata(3,:),'.') ;
        title('Raw Data') ;
        axis([-6 6 -6 6 -6 6]) ;
      hobj1 = gca ;
    subplot(2,2,2) ;
          sev1 = 3 * sqrt(veigvals(3)) * ev1 ;
          mlev1 = [(vmean - sev1) (vmean + sev1)] ;
          sev2 = 3 * sqrt(veigvals(2)) * ev2 ;
          mlev2 = [(vmean - sev2) (vmean + sev2)] ;
          sev3 = 3 * sqrt(veigvals(1)) * ev3 ;
          mlev3 = [(vmean - sev3) (vmean + sev3)] ;
      plot3(mlev1(1,:),mlev1(2,:),mlev1(3,:),'-y',...
                  mlev2(1,:),mlev2(2,:),mlev2(3,:),'-r',...
                  mlev3(1,:),mlev3(2,:),mlev3(3,:),'-b') ;
        title('P. C. Directions') ;
        axis([-6 6 -6 6 -6 6]) ;
      hobj2 = gca ;
    subplot(2,2,3) ;
          sev1s = 3 * sqrt(veigvals(3)) * ev1s ;
          mlev1s1 = [(mu1 - sev1s) (mu1 + sev1s)] ;
          mlev1s2 = [(mu2 - sev1s) (mu2 + sev1s)] ;
          mlev1s3 = [(mu3 - sev1s) (mu3 + sev1s)] ;
          sev2s = 3 * sqrt(veigvals(2)) * ev2s ;
          mlev2s1 = [(mu1 - sev2s) (mu1 + sev2s)] ;
          mlev2s2 = [(mu2 - sev2s) (mu2 + sev2s)] ;
          mlev2s3 = [(mu3 - sev2s) (mu3 + sev2s)] ;
          sev3s = 3 * sqrt(veigvals(1)) * ev3s ;
          mlev3s1 = [(mu1 - sev3s) (mu1 + sev3s)] ;
          mlev3s2 = [(mu2 - sev3s) (mu2 + sev3s)] ;
          mlev3s3 = [(mu3 - sev3s) (mu3 + sev3s)] ;
      plot3(mlev1s1(1,:),mlev1s1(2,:),mlev1s1(3,:),'-y',...
                  mlev2s1(1,:),mlev2s1(2,:),mlev2s1(3,:),'-r',...
                  mlev3s1(1,:),mlev3s1(2,:),mlev3s1(3,:),'-b',...
                  mlev1s2(1,:),mlev1s2(2,:),mlev1s2(3,:),'-y',...
                  mlev2s2(1,:),mlev2s2(2,:),mlev2s2(3,:),'-r',...
                  mlev3s2(1,:),mlev3s2(2,:),mlev3s2(3,:),'-b',...
                  mlev1s3(1,:),mlev1s3(2,:),mlev1s3(3,:),'-y',...
                  mlev2s3(1,:),mlev2s3(2,:),mlev2s3(3,:),'-r',...
                  mlev3s3(1,:),mlev3s3(2,:),mlev3s3(3,:),'-b') ;
        title('Ind. P. C. Dir''s') ;
        axis([-6 6 -6 6 -6 6]) ;
      hobj3 = gca ;


    nrotate = 10 ;
    nframe = 50 ;

    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

      subplot(2,2,1) ;
        chobj1 = get(hobj1,'Children') ;
          %  OBJECT HANDLE FOR PONT CLOUD, WHICH IS THE CHILD
          %  OF THE AXES OBJECT
        rotate(chobj1,azel,alpha,origin) ;
          %  ROTATES POINT CLOUD
        plot3(get(chobj1,'XData'),get(chobj1,'YData'),...
                  get(chobj1,'ZData'),'.') ;
        title('Raw Data') ;
        axis([-6 6 -6 6 -6 6]) ;
        hobj1 = gca ;

      subplot(2,2,2) ;
        vchobj2 = get(hobj2,'Children') ;
          %  VECTOR OF FILE HANDLE FOR LINES IN PLOT
        for ic = 1:length(vchobj2) ;
          rotate(vchobj2(ic),azel,alpha,origin) ;
        end ;
        plot3(get(vchobj2(3),'XData'),get(vchobj2(3),'YData'),...
                  get(vchobj2(3),'ZData'),'-y',...
              get(vchobj2(2),'XData'),get(vchobj2(2),'YData'),...
                  get(vchobj2(2),'ZData'),'-r',...
              get(vchobj2(1),'XData'),get(vchobj2(1),'YData'),...
                  get(vchobj2(1),'ZData'),'-b') ;
          %  THIS LONG AND CLUMSY FORM SEEMED TO BE NEEDED, BECAUSE
          %  THE USUAL hold on WOULD NOT INTERACT PROPERLY WITH THE
          %  VECTOR OF OBEJCT HANDLES, vchobj2.
        title('P. C. Directions') ;
        axis([-6 6 -6 6 -6 6]) ;
        hobj2 = gca ;

      subplot(2,2,3) ;
        vchobj3 = get(hobj3,'Children') ;
        for ic = 1:length(vchobj3) ;
          rotate(vchobj3(ic),azel,alpha,origin) ;
        end ;
        plot3(get(vchobj3(9),'XData'),get(vchobj3(9),'YData'),...
                  get(vchobj3(9),'ZData'),'-y',...
              get(vchobj3(8),'XData'),get(vchobj3(8),'YData'),...
                  get(vchobj3(8),'ZData'),'-r',...
              get(vchobj3(7),'XData'),get(vchobj3(7),'YData'),...
                  get(vchobj3(7),'ZData'),'-b',...
              get(vchobj3(6),'XData'),get(vchobj3(6),'YData'),...
                  get(vchobj3(6),'ZData'),'-y',...
              get(vchobj3(5),'XData'),get(vchobj3(5),'YData'),...
                  get(vchobj3(5),'ZData'),'-r',...
              get(vchobj3(4),'XData'),get(vchobj3(4),'YData'),...
                  get(vchobj3(4),'ZData'),'-b',...
              get(vchobj3(3),'XData'),get(vchobj3(3),'YData'),...
                  get(vchobj3(3),'ZData'),'-y',...
              get(vchobj3(2),'XData'),get(vchobj3(2),'YData'),...
                  get(vchobj3(2),'ZData'),'-r',...
              get(vchobj3(1),'XData'),get(vchobj3(1),'YData'),...
                  get(vchobj3(1),'ZData'),'-b') ;
          %  THIS IS EVEN CLUMSIER, BUT I DON'T KNOW A BETTER WAY TO
          %  GET ALL THESE ON THE PLOT
        title('Ind. P. C. Dir''s') ;
        axis([-6 6 -6 6 -6 6]) ;
        hobj3 = gca ;

      M(:,iframe) = getframe(hfig) ;
    end ;

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

    movie(hfig,M,nrotate) ;

  elseif iview == 2 ;
    hfig = figure(1) ;
    subplot(1,1,1) ;
      plot3(mdata(1,:),mdata(2,:),mdata(3,:),'.') ;
        title('Solution to Part 3') ;
        axis([-6 6 -6 6 -6 6]) ;
        hold on ;
            sev1 = 3 * sqrt(veigvals(3)) * ev1 ;
            mlev1 = [(vmean - sev1) (vmean + sev1)] ;
            sev2 = 3 * sqrt(veigvals(2)) * ev2 ;
            mlev2 = [(vmean - sev2) (vmean + sev2)] ;
            sev3 = 3 * sqrt(veigvals(1)) * ev3 ;
            mlev3 = [(vmean - sev3) (vmean + sev3)] ;
          plot3(mlev1(1,:),mlev1(2,:),mlev1(3,:),'-g',...
                mlev2(1,:),mlev2(2,:),mlev2(3,:),'-g',...
                mlev3(1,:),mlev3(2,:),mlev3(3,:),'-g') ;

            sev1s = 3 * sqrt(veigvals(3)) * ev1s ;
            mlev1s1 = [(mu1 - sev1s) (mu1 + sev1s)] ;
            mlev1s2 = [(mu2 - sev1s) (mu2 + sev1s)] ;
            mlev1s3 = [(mu3 - sev1s) (mu3 + sev1s)] ;
            sev2s = 3 * sqrt(veigvals(2)) * ev2s ;
            mlev2s1 = [(mu1 - sev2s) (mu1 + sev2s)] ;
            mlev2s2 = [(mu2 - sev2s) (mu2 + sev2s)] ;
            mlev2s3 = [(mu3 - sev2s) (mu3 + sev2s)] ;
            sev3s = 3 * sqrt(veigvals(1)) * ev3s ;
            mlev3s1 = [(mu1 - sev3s) (mu1 + sev3s)] ;
            mlev3s2 = [(mu2 - sev3s) (mu2 + sev3s)] ;
            mlev3s3 = [(mu3 - sev3s) (mu3 + sev3s)] ;
          plot3(mlev1s1(1,:),mlev1s1(2,:),mlev1s1(3,:),'-r',...
                mlev2s1(1,:),mlev2s1(2,:),mlev2s1(3,:),'-r',...
                mlev1s2(1,:),mlev1s2(2,:),mlev1s2(3,:),'-r',...
                mlev2s2(1,:),mlev2s2(2,:),mlev2s2(3,:),'-r',...
                mlev1s3(1,:),mlev1s3(2,:),mlev1s3(3,:),'-r',...
                mlev2s3(1,:),mlev2s3(2,:),mlev2s3(3,:),'-r') ;
          hold off ;


    nrotate = 10 ;
    nframe = 50 ;

    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

        hobj = gca ;

        vchobj = get(hobj,'Children') ;
           %  VECTOR OF OBJECT HANDLES FOR "CHILDREN" OF AXES,
           %  WHICH ARE THE PIECES OF THE PLOT.


        for ic = 1:length(vchobj) ;
          rotate(vchobj(ic),azel,alpha,origin) ;
        end ;      

        plot3(get(vchobj(10),'XData'),get(vchobj(10),'YData'),...
                          get(vchobj(10),'ZData'),'.',...
              get(vchobj(9),'XData'),get(vchobj(9),'YData'),...
                          get(vchobj(9),'ZData'),'-g',...
              get(vchobj(8),'XData'),get(vchobj(8),'YData'),...
                          get(vchobj(8),'ZData'),'-g',...
              get(vchobj(7),'XData'),get(vchobj(7),'YData'),...
                          get(vchobj(7),'ZData'),'-g',...
              get(vchobj(6),'XData'),get(vchobj(6),'YData'),...
                          get(vchobj(6),'ZData'),'-r',...
              get(vchobj(5),'XData'),get(vchobj(5),'YData'),...
                          get(vchobj(5),'ZData'),'-r',...
              get(vchobj(4),'XData'),get(vchobj(4),'YData'),...
                          get(vchobj(4),'ZData'),'-r',...
              get(vchobj(3),'XData'),get(vchobj(3),'YData'),...
                          get(vchobj(3),'ZData'),'-r',...
              get(vchobj(2),'XData'),get(vchobj(2),'YData'),...
                          get(vchobj(2),'ZData'),'-r',...
              get(vchobj(1),'XData'),get(vchobj(1),'YData'),...
                          get(vchobj(1),'ZData'),'-r') ;
          %  HERE I EVEN RAN INTO A PROBLEM WITH TOO MANY ARGUMENTS
          %  OF THE COMMAND plot3, WHICH IS WHY I WENT TO ONLY
          %  SHOWING THE FIRST TWO PC'S FOR THE INDIVIDUAL POINT
          %  CLOUDS.
          title('Solution to Part 3') ;
          axis([-6 6 -6 6 -6 6]) ;

      M(:,iframe) = getframe(hfig) ;
    end ;

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

    movie(hfig,M,nrotate) ;

  end ;

end ;

