% Iraqcasualties.m
% Americans killed in in Iraq, post May 1
% From the NYTimes 6/11/03, p A10.
% DATA(:,1) Days since June 30, DATA(1,1)=1 is June 1
% DATA(:,2) US Deaths per day.
% calls poiss_pdf.m, Stixbox: pchisq.m
% revised by E Gallagher 5/6/05

DATA=[1 1;2 0;3 1;4 1;5 0;6 0;7 0;8 4;9 2;10 1;11 0;12 2;
    13 2;14 1;15 0;16 1;17 0;18 2;19 6;20 0;21 1;22 0;
    23 0;24 0;25 1;26 4;27 2;28 2;29 0;30 3;31 1;32 0;
    33 0;34 0;35 1;36 1;37 1];
DEATHS=zeros(1,7); % To tally the number of deaths
for i=0:6
        k=(DATA(:,2)==i); % returns a 1 if equality met
        DEATHS(i+1)=sum(k); % sum those 1's
        % The above statements could be compressed to:
        DEATHS(i+1)=sum(DATA(:,2)==i);
end
[r,c]=size(DATA);
lambda=mean(DATA(:,2))
p=poiss_pdf(lambda,0:100);
Expdeathperday=p*r;  % this assumes no 'missing' days in the data matrix.

EXP=fliplr(Expdeathperday);
i=min(find(cumsum(EXP)>1));

EXP=fliplr([sum(EXP(1:i)) EXP(i+1:end)]);  % fliplr after concantenated to parts of vector
EXP=EXP';                   % change to a column vector
maxdeath=max(DATA(:,2));
if exist('histc')==3
    OBS = histc(DATA(:,2),0:max(DATA(:,2)));
else
    OBS=zeros(maxdeath+1,1);
    for i=0:maxdeath
        OBS(i+1)=sum(DATA(:,2)==i);
    end
end

le=length(EXP);
OBS=[OBS(1:le-1);sum(OBS(le:end))];


 disp('U.S.soldier Iraq Deaths per Day May 1 to June 7');
 disp('Deaths/day   Expected  Observed') 
 disp([[0:le-1]' EXP OBS])
 % That final category isn't deaths but 3+ deaths, using string
 % commands as shown in LMcs020401, but I won't bother here
chisquare=sum((OBS-EXP).^2 ./EXP);
df=length(EXP)-1-1; % One d.f. is lost from estimating lambda, once lambda is estimated another df is 
                    % lost because the final category can be estimated by simple arithmetic.

% p=1-gammainc(chisquare/2,df/2); % Chi square p value.
p = 1- chi2cdf(chisquare,df);
fprintf(...
'P(Chi square statistic >= %4.2f | Poisson model & %1.0f df)= %4.2f, 1-tailed\n',...
chisquare,df,p)
if exist('bar')~=0
    % Create a side-by-side bar chart. This is quite tricky in Matlab.
    deaths=[0:3]';
    h=bar(deaths,[OBS EXP],'grouped');
    lims=axis; % get axis limits, need to plot Y to 20
    axis([lims(1:2) lims(3) 20]);
    ax1=gca;% save the handle of the graph, allows customizing size & colors
    ylabel('Number of days','Color','blue','FontSize',14)
    set(ax1,'Xtick',deaths,'XColor','k')
    set(ax1,'XtickLabel',{'0 ','1 ','2','3+'},'FontSize',9)
    set(ax1,'Ytick',[0:5:20],'YColor','b'); % Manually set ticks & color
    xlabel('Deaths/day','FontSize',14)
    legend(ax1,'Observed','Expected')
    title('Iraqi Casualties May 1 to June 7, 2003','FontSize',14)
    text(2.96,17.7,sprintf('lambda = %4.2f',lambda)) % find the coordinates using [x,y]=ginput(1);
    % [x,y]=ginput(1)
    figure(gcf)   % display the current figure gcf stands for 'get current figure'
    pause     % pause until a key is hit, use this time to copy figure
    print -depsc -tiff -r300 Iraq  % Use this to save as an eps file      
    % clf   % Clear the figure when done
end
% Pool the final category too
OBS=[OBS(1:2);sum(OBS(3:end))];
EXP=[Expdeathperday(1:2) sum(Expdeathperday(3:end))]';
le=length(EXP);
 disp('Deaths/day   Expected  Observed') 
 disp([[0:le-1]' EXP OBS])
 % That final category isn't deaths but 3+ deaths, using string
 % commands as shown in LMcs020401, but I won't bother here
chisquare=sum((OBS-EXP).^2 ./EXP);
df=length(EXP)-1-1;  

% p=1-gammainc(chisqstat/2,df/2);
p = 1- chi2cdf(chisquare,df);
fprintf(...
'P(Chi square statistic >= %4.2f | Poisson model & %1.0f df)= %4.2f, 1-tailed\n',...
chisquare,df,p)
    
if exist('bar')~=0
    % Create a side-by-side bar chart. This is quite tricky in Matlab.
    deaths=[0:2]';
    h=bar(deaths,[OBS EXP],'grouped');
    lims=axis; % get axis limits, need to plot Y to 20
    axis([lims(1:2) lims(3) 20]);
    ax1=gca;% save the handle of the graph, allows customizing size & colors
    ylabel('Number of days','Color','blue','FontSize',14)
    set(ax1,'Xtick',deaths,'XColor','k')
    set(ax1,'XtickLabel',{'0 ','1 ','2+'},'FontSize',9)
    set(ax1,'Ytick',[0:5:20],'YColor','b'); % Manually set ticks & color
    xlabel('Deaths/day','FontSize',14)
    legend(ax1,'Observed','Expected')
    title('Iraqi Casualties May 1 - June 7, 2003','FontSize',14)
    text(2.1,17.7,sprintf('lambda = %4.2f',lambda)) % find the coordinates using 
    % [x1,y1]=ginput(1)
    figure(gcf)   % display the current figure gcf stands for 'get current figure'
    pause     % pause until a key is hit, use this time to copy figure
    print -depsc -tiff -r300 Iraq2  % Use this to save as an eps file
    clf   % Clear the figure when done
end