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]);