Chapter 4: binomial prior predictive and sampling

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

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

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

%%% F4.1a %%%

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

figure(5); clf; mesh(thetas(1,2:end-1),Slist(:,1),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)

%%% F4.1b %%%

pri=zeros(1,size(thetas,2))-log(N+1);

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

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

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

%%% F4.1c %%%

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

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

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

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

%%% F4.1d %%%

BinTrials=N;

N=101;

tdraw=rand(N,1);

cumpdrawn=tdraw;

tevald=0:BinTrials;

Bevald=exp(Lnow);

sumBevald=sum(Bevald);

cumpevald=cumsum(Bevald./sumBevald);

cumind=nan(1,N);

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

Bdraw=tevald(cumind);

figure(8); clf; hold on

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

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

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

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