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


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

n = 20 ;
xdat = -2 + 0.2 * randn(n,1) ;
ydat = 1.2 * randn(n,1) ;
xdat = [xdat; 3] ;
ydat = [ydat; 0] ;

subplot(2,2,1) ;
plot(xdat,ydat,'o') ;
  axis([-4,4,-4,4]) ;
  axis('square') ;
  xlabel('x1') ;
  ylabel('x2') ;
  set(gca,'XTickLabel',[]) ;
  set(gca,'YTickLabel',[]) ;
  set(get(gca,'XLabel'),'FontSize',8) ;
  set(get(gca,'YLabel'),'FontSize',8) ;


vmad = gpmad([xdat, ydat]) ;

rxdat = xdat / vmad(1) ;
rydat = ydat / vmad(2) ;

subplot(2,2,2) ;
plot(rxdat,rydat,'o') ;
  axis([-4/vmad(1),4/vmad(1),-4/vmad(2),4/vmad(2)]) ;
  axis('equal') ;
  xlabel('x1/MAD') ;
  ylabel('x2/MAD') ;
  set(gca,'XTickLabel',[]) ;
  set(gca,'YTickLabel',[]) ;
  set(get(gca,'XLabel'),'FontSize',8) ;
  set(get(gca,'YLabel'),'FontSize',8) ;



vmean = rmean([rxdat, rydat]) ;
r = 4 ;
theta = 2 * pi * (0:(1/100):1) ;
cx = vmean(1) + r * cos(theta) ;
cy = vmean(2) + r * sin(theta) ;

cxdat = rxdat - vmean(1) ;
cydat = rydat - vmean(2) ;
vlength = sqrt(cxdat.^2 + cydat.^2) ;
sxdat = vmean(1) + r * cxdat ./ vlength ;
sydat = vmean(2) + r * cydat ./ vlength ;


subplot(2,2,4) ;
plot(sxdat,sydat,'+',cx,cy,'--') ;
hold off ;
  axis([-4/vmad(1),4/vmad(1),-4/vmad(2),4/vmad(2)]) ;
  axis('equal') ;
  xlabel('x1/MAD') ;
  ylabel('x2/MAD') ;
  set(gca,'XTickLabel',[]) ;
  set(gca,'YTickLabel',[]) ;
  set(get(gca,'XLabel'),'FontSize',8) ;
  set(get(gca,'YLabel'),'FontSize',8) ;




usxdat = vmad(1) * sxdat ;
usydat = vmad(2) * sydat ;

ucx = vmad(1) * cx ;
ucy = vmad(2) * cy ;

subplot(2,2,3) ;
plot(usxdat,usydat,'+',ucx,ucy,'--') ;
hold off ;
  axis([-4,4,-4,4]) ;
  axis('square') ;
  xlabel('x1') ;
  ylabel('x2') ;
  set(gca,'XTickLabel',[]) ;
  set(gca,'YTickLabel',[]) ;
  set(get(gca,'XLabel'),'FontSize',8) ;
  set(get(gca,'YLabel'),'FontSize',8) ;


mcov = cov([usxdat, usydat]) ;
[meigv,eigval] = eig(mcov) ; 
eigval = diag(eigval) ;
[temp,visort] = sort(eigval) ;
visort = flipud(visort) ;
meigv = meigv(:,visort) ;
eigval = eigval(visort) ;


uvmean = vmad .* vmean ;

hold on ;
  plot([uvmean(1); (uvmean(1) + 2 * eigval(1) * meigv(1,1))], ...
       [uvmean(2); (uvmean(2) + 2 * eigval(1) * meigv(2,1))],'-') ;
    vhand = get(gca,'Children') ;
    set(vhand(1),'LineWidth',1.5) ;
  plot([uvmean(1); (uvmean(1) + 2 * eigval(2) * meigv(1,2))], ...
       [uvmean(2); (uvmean(2) + 2 * eigval(2) * meigv(2,2))],'-') ;
    vhand = get(gca,'Children') ;
    set(vhand(1),'LineWidth',0.5) ;

hold off ;


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




