Chapter 4: binomial with sampling

%%% Binomial Prior Predictive (F4.1) %%

Sprior=20; Nprior=30;

N=30; Slist=[0:N]; thetas=linspace(0,1,301);

%%% a %%%

Snow=Slist'*ones(1,301); thetanow=ones(size(Slist,1),1)*thetas;

Lnow=lnfactorial(N)-(lnfactorial(Snow)+lnfactorial(N-Snow))+Snow.*log(thetanow)+(N-Snow).*log(1-thetanow); Lnow(isnan(Lnow))=0; Lnow=Lnow-max(Lnow); %nan generated when p=1

figure(5); clf; mesh(thetas(2:end-1),Slist,exp(Lnow(:,2:end-1))); view(-41,44.5)

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

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

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

%%% b %%%

pri=bpdf(Sprior,Nprior,thetas);

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

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

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

%%% c %%%

pri=pri-logsum(pri,1,diff(thetas(1:2)));

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

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

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

xlabel('S (count)','FontName','Arial','FontSize',13);

%%% d %%%

N=101;

cumpdrawn=rand(N,1);

Bevald=exp(Lnow);

cumpevald=cumsum(Bevald);

cumind=nan(1,N);

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

Bdraw=Slist(cumind);

figure(8); clf; hold on

[N X]=hist(Bdraw,Slist);

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

box off; axis([-.5 30.5 0 10])

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