Chapter 4: negative binomial prior predictive and sampling

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% NegBinom Prior Predictive (F4.2) %%%

S=12; Nlist=[S:1200]'*ones(1,301); thetas=ones(size(S,1),1)*linspace(0,1,301);

%%% F4.4a %%%

Lnow=lnfactorial(Nlist-1)-(lnfactorial(S-1)+lnfactorial(Nlist-S))+S.*log(thetas)+(Nlist-S).*log(1-thetas); Lnow(isnan(Lnow))=0; Lnow=Lnow-max(Lnow); %nan generated when p=1

Lnow(isnan(Lnow))=0;

figure(5); clf; mesh(thetas(1,2:end),Nlist(:,1),exp(Lnow(:,2:end))); view(-9.8,67.6)

xlabel('\theta (rate)','FontName','Arial','FontSize',13)

ylabel('N (count)','FontName','Arial','FontSize',13)

zlabel('p','FontName','Arial','FontSize',13)

%%% F4.2b %%%

pri=-.5*(log(thetas(1,:))+log(1-thetas(1,:))); pri(isinf(pri))=nan; pri=pri-max(max(pri)); pri(isnan(pri))=0;

figure(6); clf; plot(thetas,exp(pri),'ko','MarkerFaceColor','k','MarkerSize',8); r=axis; axis([r(1:2) 0 1]); box off

xlabel('\theta (rate)','FontName','Arial','FontSize',13)

ylabel('p','FontName','Arial','FontSize',13)

%%% F4.2c %%%

Lnow=logsum(Lnow+ones(size(Nlist,1),1)*pri,2,diff(thetas(1,1:2)));

figure(7); clf; plot(Nlist(:,1),exp(Lnow),'ko','MarkerFaceColor','k','MarkerSize',8);

box off; ylabel('p','FontName','Arial','FontSize',13)

xlabel('N (count)','FontName','Arial','FontSize',13)

%%% F4.2d %%%

reps=101;

tdraw=rand(reps,1);

cumpdrawn=tdraw;

tevald=Nlist(:,1);

Bevald=exp(Lnow);

sumBevald=sum(Bevald);

cumpevald=cumsum(Bevald./sumBevald);

cumind=nan(1,reps);

for nbr=1:reps, cumind(nbr)=findnearestN(cumpevald,cumpdrawn(nbr),1); end

Bdraw=tevald(cumind);

figure(8); clf; hold on

[N X]=hist(Bdraw,Nlist(:,1));

bar(X,N,'FaceColor',.25*[1 1 1])

box off; axis([min(min(Nlist))-.5 max(max(Nlist))+.5 0 3])

xlabel('N (count)','FontName','Arial','FontSize',13)