diagnosticfig.m

function diagnosticfig(N,dataname)

h=figure(1); %diagnostics plot
clf
set(h,'Position',[82 218 574 626]);


warning off;
lenN=length(N);
Nlogratios=log(N(2:lenN)./N(1:(lenN-1)));
tmp=1:length(Nlogratios);
Nlogratiosyrs=tmp(~isinf(Nlogratios) & ~isnan(Nlogratios)); %indices where Nlogratios is not zero; need for trend test
Nlogratios=Nlogratios(~isinf(Nlogratios) & ~isnan(Nlogratios)); %chuck cases where N is zero
warning on;
year = 1:lenN;

%Plot interpolated data
subplot(2,2,1)
plot(year,N,'-r');
xlabel('year');
ylabel('count');
title('Time series');
tmppos=get(gca,'Position');
tmppos(4)=.8*tmppos(4);
set(gca,'Position',tmppos);

%Test for trends in data
subplot(2,2,3)
plot(Nlogratiosyrs,Nlogratios,'or');
[a,b,c,d,stats]=regress(Nlogratios',[ones(length(Nlogratios),1) Nlogratiosyrs']);
trendpvalue=stats(3);
title(['p.value for trend test:' num2str(trendpvalue)])
xlabel('year');
ylabel('ln(N_{t+1}/N_t)');
tmppos=get(gca,'Position');
tmppos(4)=.8*tmppos(4);
set(gca,'Position',tmppos);

%check for normality
subplot(2,2,4)
normplot(Nlogratios)
tmppos=get(gca,'Position');
tmppos(4)=.8*tmppos(4);
set(gca,'Position',tmppos);

%Print out params on plot
subplot(2,2,2)
set(gca,'Visible','off','XLim',[0 1],'YLim',[0 1]);
kalman = get_kalman_ests(N);
[reml.mu aa reml.s2p reml.s2s] = remltrend(N');
str(1)={'Kalman Estimates'};
str(2)={' '};
str(3)={['\mu = ' num2str(kalman.mu,3) '  \sigma^2 _p= ' num2str(kalman.s2p,3)]};
str(4)={' '};
str(5)={['\sigma^2_{np} = ' num2str(kalman.s2s,3)]};
str(6)={' '};
str(7)={' '};
str(8)={'REML Estimates'};
str(9)={' '};
str(10)={['\mu = ' num2str(reml.mu,3) '  \sigma^2 _p= ' num2str(reml.s2p,3)]};
str(11)={' '};
str(12)={['\sigma^2_{np} = ' num2str(reml.s2s,3)]};
text(0,.5,str,'FontSize',8);

%Put title at the top of the whole graph
h=axes('Position',[0 0 1 1],'Visible','off');
text(.5,.98,dataname,'HorizontalAlignment','center','FontSize',16)

set(gcf,'PaperPosition',[.5 1 7.5 9]);