Chapter 6: difference of binomial rates

Recalling the method of computing differences of probability coordinates described in Appendix C, we will use that routine to obtain the posterior over differences to compute the model likelihoods and evidence for the current problem.

clear

plotFLAG=0; %#ok<*UNRCH>

N=200; S=[80 71]; Ntics=101; lambda=2.6721; Tmax=1;

thetas=linspace(0,1,Ntics+2)'; thetas=thetas(2:end-1);

s=struct;

%compute 2D Jeffreys prior over differences of rates

jef=[(thetas.*(1-thetas)).^-.5 thetas]; [~, ~, s.priJ]=pSumDiff(jef,jef,[],1,1);

eval(s.priJ.evalprep('s.priJ')); priSDj=pSD; %extract pSD, icell, Dmat, Smat

%%% separate prior; I0 consistent with M1, Ip consistent with M2

Ip=Dmat>0; I0=Dmat==0; In=Dmat<=0;

priSDJ0=priSDj.*I0/sum(sum(priSDj.*I0)); pmax=max(max(priSDJ0));

priSDJp=priSDj.*Ip/sum(sum(priSDj.*Ip)); pmax(2)=max(max(priSDJp));

if plotFLAG, figure(3); clf; hold on

mesh(Dmat,Smat,priSDJ0); colormap(bone)

plotmat([Dmat(:),Smat(:),priSDJ0(:)],'k.'); r0=axis;

figure(4); clf; hold on

mesh(Dmat,Smat,priSDJp); colormap(bone)

plotmat([Dmat(:),Smat(:),priSDJp(:)],'k.'); rp=axis;

for f=3:4, figure(f); axis([r0(1:4) 0 .01]); view(-20,10); end, end

Stest=0:N;

EVtest=nan(size(Stest));

for Snow=0:200,

vbp1=[bpdf(Snow,N,thetas) thetas];

vbp2=[bpdf(S(2),N,thetas) thetas];

%%% compute likelihoods of differences and sums of binomial rates

[~,~, s.like]=pSumDiff(vbp1,vbp2,[],1,1);

eval(s.like.evalprep('s.like')); likeSD=pSD;

pSD=priSDJ0.*likeSD; pSDJ0=pSD;

eval(s.like.evals{2}); pDJ=pD(:,1); %compute pD from pSD, Dmat

pSD=priSDJp.*likeSD; pSDJp=pSD;

eval(s.like.evals{2}); pDJ=[pDJ pD]; %compute pD from pSD, Dmat

ML=sum(pDJ(:,1:2));

%odds computation

o=ML(:,2)./ML(:,1);

%%% evidence

EVtest(Snow+1)=10*log10(o);

if and(plotFLAG,Snow==S(1)),

EV=10*log10(o); %EV(1)=.73, for uniform over delta

figure(5); clf; hold on

mesh(Dmat,Smat,pSDJ0); colormap(bone); plot3(Dmat(:),Smat(:),pSDJ0(:),'k.');

plot3(str8n(Dmat(I0)),str8n(Smat(I0)),pSDJ0(I0),'o','MarkerFaceColor',.65*[1 1 1],'MarkerEdgeColor',.6*[1 1 1],'MarkerSize',4)

plotmat([pDJ(:,3) max(Smat(:))*ones(size(pDJ(:,3))) max(max(pSDJ0))*pDJ(:,1)/max(pDJ(:,1))],'ko','MarkerFaceColor','k','MarkerSize',10);

plot3([0 0],max(Smat(:))*[1 1], [0 max(max(pSDJ0))],'k-')

ind0=find(pDJ(:,3)~=0);

plotmat([pDJ(ind0,3) max(Smat(:))*ones(size(ind0)) zeros(size(ind0))],'ko','MarkerFaceColor',.65*[1 1 1],'MarkerSize',8);

a=6.5; b=48.5; view(a,b); r2=axis; r2(1:2)=[-.8 .8]; axis(r2)

figure(6); clf; hold on

mesh(Dmat,Smat,max(max(pSDJ0))*pSDJp/max(max(pSDJp))); colormap(bone);

plot3(Dmat(:),Smat(:),max(max(pSDJ0))*pSDJp(:)/max(max(pSDJp)),'k.');

plot3(str8n(Dmat(Ip)),str8n(Smat(Ip)),zeros(sum(Ip(:))),'o','MarkerFaceColor',.6*[1 1 1],'MarkerEdgeColor',.65*[1 1 1],'MarkerSize',4)

plotmat([pDJ(:,3) max(Smat(:))*ones(size(pDJ(:,3))) max(max(pSDJ0))*pDJ(:,2)/max(pDJ(:,2))],'ko','MarkerFaceColor','k','MarkerSize',8);

ind0=find(pDJ(:,3)<=0);

plotmat([pDJ(ind0,3) max(Smat(:))*ones(size(ind0)) zeros(size(ind0))],'ko','MarkerFaceColor',.65*[1 1 1],'MarkerSize',8);

view(-a,b); axis(r2)

figure(7); clf; stem([-1 1],ML,'ko','MarkerFaceColor','k','MarkerSize',15)

r=axis; axis([-1.5 1.5 r(3:4)]); box off; end, end

figure(8); clf; hold on

plot([0 N],[0 0],'k--')

plot(Stest(EVtest>0),EVtest(EVtest>0),'ko','MarkerFaceColor',.7*[1 1 1],'MarkerSize',11)

plot(Stest(EVtest<0),EVtest(EVtest<0),'ko','MarkerFaceColor',.4*[1 1 1],'MarkerSize',11)

r=axis;

plot((Stest(find(EVtest>0,1))-.5)*[1 1],r(3:4),'k:')

plot((Stest(Stest==S(2))-.5)*[1 1],r(3:4),'k:')

xlabel('possible recovery data','FontName','Arial','FontSize',14)

ylabel('evidence (dB)','FontName','Arial','FontSize',14)