disp('Running MATLAB script file lmtzcf31.m') ;
%
%    Makes Figure 3.1, for the Locantore, Marron, Tripoli,
%    Zhang, Cohen paper.
%
%    Copied from cs1.m  (part 3)


seed = 29030987 ;
rand('seed',seed) ;
randn('seed',seed) ;

n = 20 ;
xdat = randn(n,1) ;
%ydat = xdat + 0.7 * randn(n,1) ;
ydat = xdat + 0.4 * randn(n,1) ;
 
plot(xdat,ydat,'o') ;
  axis([-2,2,-2,2]) ;
  axis('square') ;
  set(gca,'XTickLabels',[]) ;
  set(gca,'YTickLabels',[]) ;


vmean = mean([xdat, ydat]) ;
mcov = cov([xdat, ydat]) ;
[meigv,eigval] = eig(mcov) ; 
eigval = diag(eigval) ;
[temp,visort] = sort(eigval) ;
visort = flipud(visort) ;
meigv = meigv(:,visort) ;
eigval = eigval(visort) ;


line1x = [-2; 2] ; 
line1y = vmean(2) + (meigv(2,1) / meigv(1,1)) * (line1x - vmean(1)) ;

mdat = [xdat, ydat] ;
datp1 = (mdat - vec2mat(vmean,n)) * ...
             meigv(:,1) / sqrt(meigv(:,1)' * meigv(:,1)) ;



line2y = [-2; 2] ; 
line2x = vmean(1) + (meigv(1,2) / meigv(2,2)) * (line2y - vmean(2)) ;

datp2 = (mdat - vec2mat(vmean,n)) * vec2mat(meigv(:,2),n) ./ ...
                   sqrt(meigv(:,2)' * meigv(:,2)) ;


%  np = [.2; -1.2] ;
%  npp1 = (np - vmean')' * meigv(:,1) / sqrt(meigv(:,1)' * meigv(:,1)) ;



hold on ;
  plot(line1x,line1y,'-') ;
    vhand = get(gca,'Children') ;
    set(vhand(1),'LineWidth',1.5) ;
  plot(line2x,line2y,'-') ;
    vhand = get(gca,'Children') ;
%    set(vhand(1),'LineWidth',1) ;

  plot((vmean(1) + datp1 * meigv(1,1)), ...
           (vmean(2) + datp1 * meigv(2,1)),'r+') ;
    vhand = get(gca,'Children') ;
    set(vhand(1),'MarkerSize',9) ;
    set(vhand(1),'LineWidth',1.5) ;
  plot((vmean(1) + datp2 * meigv(1,2)), ...
           (vmean(2) + datp2 * meigv(2,2)),'r+') ;
    vhand = get(gca,'Children') ;
%    set(vhand(1),'LineWidth',1) ;

hold off ;



  orient tall ;
  set(gcf,'PaperPosition',[2.5, 4.5, 4.5, 4.5]) ;
        %  2.0 inches on left and right, 4.0 interior, 8.5 total
        %  4.5 from bottom, 4.5 interior, 11 total
print -deps \matlab\steve\ps\lmtzcf31.eps ;


