Chapter 2: multinomial sampling-likelihood surface

The multinomial distribution only differs from the binomial distribution when there are more than two count-types (as would occur when tracking the occurrence of multiple behavior types, such as the set of symptoms that correspond to a particular disease). This creates a bit of a problem for us, because we can only represent a 2D probability surface on the page, so adding count-dimensions in addition to what was needed to represent the binomial in Fig. 2.5 is not possible. We can, however, circumvent the issue by setting some of the additional variables constant. To create the surface in Fig. 2.8a, we limit the total number of trials to 44, and the first rate constant to 0.2.

%set constants

THETA=[.2 0 0]; N=22; Nmax=44; Htheta=.4125;

thetadeltsize=.025;

sumtheta=sum(THETA); offtheta=THETA(find(THETA~=0));

theta=[thetadeltsize/2:thetadeltsize:1-sumtheta-thetadeltsize/2]';

thetacomp=1-theta-sumtheta;

Ns=0:Nmax; Nstr='Ns'; Xstr='X{1}'; Istr='i(:,1)'; INstr='i(n,1)';

for n=2:length(THETA),

Nstr=[Nstr ',Ns'];

Xstr=[Xstr ' X{' num2str(n) '}'];

Istr=[Istr ' i(:,' num2str(n) ')'];

INstr=[INstr ',i(n,' num2str(n) ')']; end

eval(['[' Xstr ']=meshgrid(' Nstr ');']); %allows us to compute meshgrid for different numbers of multinomial counts

Z=X{1}; for n=2:length(THETA), Z=Z+X{n}; end

[ijk]=find(Z~=Nmax); [i j k]=ind2sub(size(Z),ijk);

for n=1:length(ijk), X{1}(i(n),j(n),k(n))=-1; end

eval('NX=size(X{1}); i=[]; ijk=[];');

THETA=[theta thetacomp];

for n=1:length(offtheta), THETA(:,end+1)=offtheta(n)*ones(size(theta)); end

for Nrun=Ns, eval('i=[];')

[ijk]=find(X{1}==Nrun);

eval(['[' Istr ']=ind2sub(NX,ijk);']);

MDtemp=[];

for n=1:length(ijk),

for nn=1:size(i,2), eval(['Nnow(1,nn)=X{nn}(' INstr ');']); end

MDtemp(n,:)=mnpdf(repmat(Nnow,[length(theta) 1]),THETA); end

MD(Nrun+1,:)=sum(MDtemp,1); end

check1=sum(MD,2);

figure(1); hold on

[X1 X2]=meshgrid(theta,Ns);

mesh(X1,X2,MD); colormap(bone); brighten(-1); view(-32,32);

axis([0 1 0 Nmax 0 .35])

delttheta=abs(theta-Htheta); deltN=abs(Ns-N);

indTheta=min(find(delttheta==min(delttheta)));

indN=min(find(deltN==min(deltN)));

plot3(X1(:,indTheta),X2(:,indTheta),max(max(MD))*.005+MD(:,indTheta),'ko','MarkerFaceColor','w','MarkerSize',8,'LineWidth',1.4)

plot3(X1(:,indTheta),X2(:,indTheta),zeros(size(X1(:,indTheta))),'k-','LineWidth',2)

plot3(X1(indN,:),X2(indN,:),max(max(MD))*.005+MD(indN,:),'ko','MarkerFaceColor','k','MarkerSize',8,'LineWidth',1.4)

plot3(X1(indN,:),X2(indN,:),zeros(size(X1(indN,:))),'k-','LineWidth',2)

xlabel('\theta_1','FontName','Times','FontSize',14)

ylabel('N_1','FontName','Times','FontSize',14)

zlabel('probability','FontName','Times','FontSize',14)

figure(2); hold on

aleph=sum(MD(:,indTheta));

for n=1:length(Ns),

plot([Ns(n) Ns(n)],[0 MD(n,indTheta)/aleph],'k:','LineWidth',1.7); end

plot(Ns,MD(:,indTheta)/aleph,'ko','MarkerFaceColor','w','MarkerSize',8,'LineWidth',1.9)

plot(indN-1,MD(indN,indTheta)/aleph,'ko','MarkerFaceColor','k','MarkerSize',8,'LineWidth',1.9)

xlabel('N_1','FontName','Times','FontSize',14)

r=axis; mprob=r(4);

figure(3); hold on

aleph=trapz(MD(indN,:));

for n=1:length(theta),

plot([theta(n) theta(n)],[0 MD(indN,n)/aleph],'k:','LineWidth',1.7); end

plot(theta,MD(indN,:)/aleph,'ko','MarkerFaceColor','k','MarkerSize',8,'LineWidth',1.9)

xlabel('\theta','FontSize',16)

r=axis;

mprob=max([mprob r(4)]);

axis([0 1 0 mprob])

figure(2); r=axis; axis([r(1) r(2) r(3) mprob])

The resulting figure is essentially a 2D ‘slice’ through the full multinomial probability distribution.