function fig10_5

%  this needs data file crime_data_2009.txt
%  it also needs subaxis.m and parseArgs.m

%  finds best line fit, then plots test data for variable for column IY
%  as function of variable from column IX
%  Columns: 1=Pop; 2=Murder; 3=Rape; 4=Robbery; 5=Assault; 6=Burglary;
%          7=Larceny; 8=Motor vehicle theft

% pick prediction variable
IX=2

load crime_data_2009.txt

AA=crime_data_2009(:,:);
[m n]=size(AA)
means=mean(AA)

%  split using every other row A=training data   B = testing data
i=0;
for ii=1:2:m-1
    i=i+1;
    A(i,:)=AA(ii,:);
    B(i,:)=AA(ii+1,:);
end
[m n]=size(A)
means=mean(A)

% randomize rows then split D into training set A and testing set B
% rows=randperm(m);
% fac=2/3;
% mm=round(fac*m);
% A=AA(rows(1:mm),:);
% B=AA(rows(mm+1:m),:);
% [m n]=size(A)
% means=mean(A)

% normalization of A
D=zeros(m,n);
for j=1:n
    D(:,j)=A(:,j)-means(j)*ones(m,1);
end
DD=zeros(m,n);
for j=1:n
    %DM(j)=max(abs(D(:,j)));
    DM(j)=norm(D(:,j),2)/sqrt(n);
    DD(:,j)=D(:,j)/DM(j);
end

%  DD=U*S*V'
[U,S,V] = svd(DD);
ss=size(S)

for i=1:ss(2)
    id(i)=i;
    sigmas(i)=S(i,i);
    %    	fprintf(' %i   %19.14e \n',i,S(i,i));
end
sigmas

V(:,1)

%  normalize testing data (using same normalization as used for A)
[M N]=size(B)
DB=zeros(M,N);
for j=1:n
    DB(:,j)=B(:,j)-means(j)*ones(M,1);
end
DDB=zeros(M,N);
for j=1:n
    DDB(:,j)=DB(:,j)/DM(j);
end


clf
% get(gcf)
set(gcf,'Position', [4 715 700 630]);
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)

% fontsize
fs=14;

%  IY choices used in different plots
subaxis(2,2,2,2,'MT',0.008,'MB',0.06,'MR',0.003,'ML',0.042,'P',0.03)
IY=1;
% calculate best line
slope=V(IY,1)/V(IX,1);
xL=min(DDB(:,IX));
yL=slope*xL;
xR=max(DDB(:,IX));
yR=slope*xR;
plot([xL xR],[yL yR],'r','LineWidth',1.6)
hold on
plot(DDB(:,IX),DDB(:,IY),'ob','MarkerSize',8,'LineWidth',2)
LX=min(DDB(:,IX)); RX=max(DDB(:,IX)); LY=min(DDB(:,IY)); RY=max(DDB(:,IY));
% LX=-0.2; RX=0.5; LY=-0.2; RY=0.6;
axis([LX RX LY RY])
axis square
xlabel('Murder')
ylabel('Population')
grid on
box on
set(gca,'FontSize',fs,'FontWeight','bold')

% format shortG
% mean squared prediction error
MSPE=norm(slope*DDB(:,IX)-DDB(:,IY),2)^2/length(DDB(:,IY));
% mean absolute deviation
MAD=norm(slope*DDB(:,IX)-DDB(:,IY),1)/length(DDB(:,IY));
% coefficient of determination R^2
R2=(dot(slope*DDB(:,IX),DDB(:,IY))/(norm(DDB(:,IX),2)*norm(DDB(:,IY),2)))^2;
errors_pop=[MSPE  MAD  R2];
fprintf('\n Population  &  %3.2e  &  %3.2e  &  %10.3f',MSPE,MAD,R2)

subaxis(2,2,1,2)
IY=4;
% calculate best line
slope=V(IY,1)/V(IX,1);
xL=min(DDB(:,IX));
yL=slope*xL;
xR=max(DDB(:,IX));
yR=slope*xR;
plot([xL xR],[yL yR],'r','LineWidth',1.6)
hold on
plot(DDB(:,IX),DDB(:,IY),'ob','MarkerSize',8,'LineWidth',2)
LX=min(DDB(:,IX)); RX=max(DDB(:,IX)); LY=min(DDB(:,IY)); RY=max(DDB(:,IY));
% LX=-0.2; RX=0.5; LY=-0.2; RY=0.6;
axis([LX RX LY RY])
axis square
xlabel('Murder')
ylabel('Robbery')
grid on
box on
set(gca,'FontSize',fs,'FontWeight','bold');
% mean squared prediction error
MSPE=norm(slope*DDB(:,IX)-DDB(:,IY),2)^2/length(DDB(:,IY));
% mean absolute deviation
MAD=norm(slope*DDB(:,IX)-DDB(:,IY),1)/length(DDB(:,IY));
% coefficient of determination R^2
R2=(dot(slope*DDB(:,IX),DDB(:,IY))/(norm(DDB(:,IX),2)*norm(DDB(:,IY),2)))^2;
errors_ass=[MSPE  MAD  R2];
fprintf('\n Robbery  &  %3.2e  &  %3.2e  &  %10.3f',MSPE,MAD,R2)

subaxis(2,2,2,1)
IY=5;
% calculate best line
slope=V(IY,1)/V(IX,1);
xL=min(DDB(:,IX));
yL=slope*xL;
xR=max(DDB(:,IX));
yR=slope*xR;
plot([xL xR],[yL yR],'r','LineWidth',1.6)
hold on
plot(DDB(:,IX),DDB(:,IY),'ob','MarkerSize',8,'LineWidth',2)
LX=min(DDB(:,IX)); RX=max(DDB(:,IX)); LY=min(DDB(:,IY)); RY=max(DDB(:,IY));
% LX=-0.2; RX=0.5; LY=-0.2; RY=0.6;
axis([LX RX LY RY])
axis square
xlabel('Murder')
ylabel('Assault')
grid on
box on
set(gca,'FontSize',fs,'FontWeight','bold');
% mean squared prediction error
MSPE=norm(slope*DDB(:,IX)-DDB(:,IY),2)^2/length(DDB(:,IY));
% mean absolute deviation
MAD=norm(slope*DDB(:,IX)-DDB(:,IY),1)/length(DDB(:,IY));
% coefficient of determination R^2
R2=(dot(slope*DDB(:,IX),DDB(:,IY))/(norm(DDB(:,IX),2)*norm(DDB(:,IY),2)))^2;
errors_ass=[MSPE  MAD  R2];
fprintf('\n Assault  &  %3.2e  &  %3.2e  &  %10.3f',MSPE,MAD,R2)

subaxis(2,2,1,1)
IY=8;
% calculate best line
% calculate best line
slope=V(IY,1)/V(IX,1);
xL=min(DDB(:,IX));
yL=slope*xL;
xR=max(DDB(:,IX));
yR=slope*xR;
plot([xL xR],[yL yR],'r','LineWidth',1.6)
hold on
plot(DDB(:,IX),DDB(:,IY),'ob','MarkerSize',8,'LineWidth',2)
LX=min(DDB(:,IX)); RX=max(DDB(:,IX)); LY=min(DDB(:,IY)); RY=max(DDB(:,IY));
% LX=-0.2; RX=0.5; LY=-0.2; RY=0.6;
axis([LX RX LY RY])
axis square
xlabel('Murder')
ylabel('Vehicle')
grid on
box on
set(gca,'FontSize',fs,'FontWeight','bold')
% mean squared prediction error
MSPE=norm(slope*DDB(:,IX)-DDB(:,IY),2)^2/length(DDB(:,IY));
% mean absolute deviation
MAD=norm(slope*DDB(:,IX)-DDB(:,IY),1)/length(DDB(:,IY));
% coefficient of determination R^2
R2=(dot(slope*DDB(:,IX),DDB(:,IY))/(norm(DDB(:,IX),2)*norm(DDB(:,IY),2)))^2;
errors_ass=[MSPE  MAD  R2];
fprintf('\n Vehicle  &  %3.1e  &  %3.1e  &  %7.3f \n\n',MSPE,MAD,R2)

%exportgraphics(gcf,'/Users/mark/Desktop/murder.eps')





























