% Shark attack
% Written by E. Gallagher for ECOS601 final
% Revised: 12/18/01
% Data from Florida museum of Natural History Shark attack page
% As linked to Journal of Statistical Education Volume 8
clf
Year=[1905:10:1995]';
Shark=[1 3 3.5 4 1.5 12 30 27.5 57 125.5]'*200/135;
Total_attacks=sum(Shark) % The table is for 392 shark attacks, this is close
% 392.5926 is the sum, close to the 392 in the original shark attack
Pop=[4 6 7.5 11.5 15 21 37 52 75 120]'*18/135*1e6;
% Fit of the Poisson model to the data.
% Poisson parameter is total shark attacks divided by total population
SharkP=Total_attacks/sum(Pop);  % This gives the probability of attack per person
Expected=SharkP*Pop;
C=sum((Shark-Expected).^2./Expected)
% C=46.8, df=8
% The degrees of freedom is the number of categories (decades)- the number
% of estimated parameters for the fitted model minus 1 (actually
% a better definition is the number of categories minus the number
% of parameters estimated from the data).  For this problem, one
% must estimate the binomial probability, SharkP, from the data and one
% estimates the sample total, hence 2 degrees of freedom are
% lost.  Since the total must be estimated in all goodness of fit
% tests the books recommendation is correct: df=k-r-1, where k
% is the number of categories, r is the number of parameters.  The
% other 1 is because the sample total must be estimated too. 
v=length(Expected)-2
% Using Matlab's built-in incomplete gamma function, which returns the
% probability of the chi square distribution.  The first input argument is
% the chi-square test statistic divided by 2 and the 2nd input parameter
% is the degrees of freedom divided by 2.  Larsen & Marx discuss this
% in the comment to Definition 7.3.1 on the top of page 411.
Chipf=1-gammainc(C/2,v/2) % p(Chisquared <= 47,8 df)=1.7e-7;
Indyear=[0:length(Year)-1]';
X=[Indyear Pop Shark]; % save the data in a format to be read into SPSS
save sharkdata  X -ASCII -DOUBLE -TABS
%Figures
% Create an overlaid bar plot and line plot as in the original
clf
h=bar(Year,Shark,.4);
ax1=gca; % save the handle for the axes of the first graph
ylabel('Number of Shark Attacks','Color','blue','FontSize',14)
set(ax1,'Xtick',Year,'XColor','k')
set(ax1,'XtickLabel',{'1900-10','1910-20','1920-30','1930-40',...
        '1940-50','1950-60', '1960-70',...
        '1970-80','1980-90','1990-2000'},'FontSize',9)
set(ax1,'Ytick',[0:40:200],'YColor','b')
xlabel('Year','FontSize',14)
% Plot the 2nd graph at the same position, using the handle ax1
ax2=axes('Position',get(ax1,'Position'),...
    'YaxisLocation','right','Color','none','XTickLabel',[],...
    'YColor','r');
h2=line(Year,Pop,'Color','r','Parent',ax2,'Linewidth',3)
set(h2,'marker','square','markersize',5,...
                 'markeredgecolor','r','markerfacecolor','r',...
                 'linestyle','-','color','r','linewidth',3)
set(ax2,'YLim',[0 16e6],'YTick',[0:16e6/5:16e6],'YColor','r')
ylabel('Florida population','Color','red','FontSize',14)
title('Florida Shark Attacks and Population','FontSize',14)
figure(gcf)
pause
print -depsc -tiff -r300 sharkf1
clf
% Create a semilog population & shark graph
plot(Year,Shark,'-bs','LineWidth',3,'MarkerEdgeColor','b',...
    'MarkerFaceColor','b','MarkerSize',3)
ax1=gca;
set(ax1,'Yscale','log')
ylabel('Number of Shark Attacks','Color','blue','FontSize',14)
set(ax1,'Xtick',Year,'XColor','k','YColor','b')
set(ax1,'XtickLabel',{'1900-10','1910-20','1920-30','1930-40',...
        '1940-50','1950-60', '1960-70',...
        '1970-80','1980-90','1990-2000'},'FontSize',9)
xlabel('Year','FontSize',14)
% Plot the 2nd graph at the same position, using the handle ax1
ax2=axes('Position',get(ax1,'Position'),...
    'YaxisLocation','right','Color','none','XTickLabel',[],...
    'YColor','r');
h2=line(Year,Pop,'Color','r','Parent',ax2,'Linewidth',3)
set(ax2,'Yscale','log')
set(ax2,'XLim',get(ax1,'XLim'),'Layer','top')
set(h2,'marker','square','markersize',5,...
                 'markeredgecolor','r','markerfacecolor','r',...
                 'linestyle','-','color','r','linewidth',3)
set(ax2,'YColor','r')
ylabel('Florida population','Color','red','FontSize',14)
title('Florida Shark Attacks and Population','FontSize',14)
figure(gcf)
pause
print -depsc -tiff -r300 sharkf2
clf


% Create an overlaid bar plot and line plot with observed and expected.
bar(Year ,[Shark Expected],'grouped')
ax1=gca;% save the handle of the first graph
set(ax1,'Xtick',Year,'XColor','k','Ytick',[0:40:200],'YColor','b')
ylabel('Number of Shark Attacks','Color','blue','FontSize',14)
set(ax1,'XtickLabel',{'1900-10','1910-20','1920-30','1930-40',...
        '1940-50','1950-60', '1960-70',...
        '1970-80','1980-90','1990-2000'},'FontSize',9)
xlabel('Year','FontSize',14)
% Plot the 2nd graph at the same position, using the handle ax1
ax2=axes('Position',get(ax1,'Position'),...
    'YaxisLocation','right','Color','none','XTickLabel',[],...
    'YColor','r');
h2=line(Year,Pop,'Color','r','Parent',ax2,'Linewidth',3)
set(h2,'marker','square','markersize',5,...
                 'markeredgecolor','r','markerfacecolor','r',...
                 'linestyle','-','color','r','linewidth',3)
set(ax2,'YLim',[0 16e6],'YTick',[0:16e6/5:16e6],'YColor','r')
ylabel('Florida population','Color','red','FontSize',14)
title('Florida Shark Attacks and Population','FontSize',14)
legend(ax1,'Observed shark attacks','Expected shark attacks',2)
figure(gcf)
pause
print -depsc -tiff -r300 sharkf3
clf

% Now plot the probability of shark attack.
h2=line(Year,Shark./Pop,'Color','b','Linewidth',3);
set(h2,'marker','o','markersize',6,...
    'MarkerEdgeColor','b',...
    'MarkerFaceColor','b');
hold on
    plot(Year,repmat(SharkP,size(Year)),'--g','Linewidth',5)
ax1=gca;
set(ax1,'Xtick',Year,'XColor','k')
set(ax1,'XtickLabel',{'1900-10','1910-20','1920-30','1930-40',...
        '1940-50','1950-60', '1960-70',...
        '1970-80','1980-90','1990-2000'},'FontSize',9)
ylabel('Probability of shark attack','FontSize',14)
xlabel('Year','FontSize',14)
title('Probability of Shark Attack in Florida','FontSize',14) 
figure(gcf)
pause
print -depsc -tiff -r300 sharkf4
hold off
clf

Odds=Shark./Pop./(1-(Shark./Pop));
h2=line(Year,Odds,'Color','b','Linewidth',3);
set(h2,'marker','o','markersize',6,...
    'MarkerEdgeColor','b',...
    'MarkerFaceColor','b');
hold on
    plot(Year,repmat(SharkP/(1-SharkP),size(Year)),'--g','Linewidth',5)
ax1=gca;
set(ax1,'Xtick',Year,'XColor','k')
set(ax1,'XtickLabel',{'1900-10','1910-20','1920-30','1930-40',...
        '1940-50','1950-60', '1960-70',...
        '1970-80','1980-90','1990-2000'},'FontSize',9)
ylabel('Odds of shark attack','FontSize',14)
xlabel('Year','FontSize',14)
title('Odds of Shark Attack in Florida','FontSize',14) 
figure(gcf)
pause
print -depsc -tiff -r300 sharkf5
hold off
clf

RelOdds=Odds./repmat(Odds(1),size(Odds))
h2=line(Year,RelOdds,'Color','b','Linewidth',3);
set(h2,'marker','o','markersize',6,...
    'MarkerEdgeColor','b',...
    'MarkerFaceColor','b');
hold on
    plot(Year,repmat(1,size(Year)),'--m','Linewidth',5)
ax1=gca;
set(ax1,'Xtick',Year,'XColor','k')
set(ax1,'XtickLabel',{'1900-10','1910-20','1920-30','1930-40',...
        '1940-50','1950-60', '1960-70',...
        '1970-80','1980-90','1990-2000'},'FontSize',9)
ylabel('Odds of shark attack relative to 1900-1910','FontSize',14)
xlabel('Year','FontSize',14)
title('Relative Odds of Shark Attack in Florida','FontSize',14) 
figure(gcf)
pause
print -depsc -tiff -r300 sharkf6
hold off

 [Tau,prob]=kendall(Year,RelOdds)
 % Probability=0.0143