Chapter 2: poisson sampling-likelihood surface

The Poisson distribution has a naturally 2D sampling-likelihood surface, so we need not make any assumptions beyond what is necessary to define the Poisson. The 2D sampling-likelihood surface requires us to compute Poisson probabilities of a range of combinations of observed count (data) and theoretical count rates.

function [X Y PD]=BuildLikelihoodPoisson(ThetaLimits,Nmax,Hlambda,Hn)

if nargin==0, thetadeltsize=1; ThetaLimits=[0 50]; Nmax=round(ThetaLimits(2)*1.2); Hlambda=20; Hn=30; end

ThetaMIN=ThetaLimits(1); ThetaMAX=ThetaLimits(2);

Lambda=ThetaMIN+thetadeltsize/2:thetadeltsize:ThetaMAX-thetadeltsize/2; N=0:Nmax;

[X Y]=meshgrid(Lambda,N); PD=[];

for n=1:length(Lambda),

PD(:,n)=ppdf(N,Lambda(n)); end

figure(1); hold on

mesh(X,Y,PD); colormap(bone); brighten(-1); view(-33,40);

deltlambda=abs(Lambda-Hlambda); deltN=abs(N-Hn);

indLambda=min(find(deltlambda==min(deltlambda)));

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

plot3(X(:,indLambda),Y(:,indLambda),max(max(PD))*.005+PD(:,indLambda),'ko','MarkerFaceColor','w','MarkerSize',8,'LineWidth',1.4);

plot3(X(:,indLambda),Y(:,indLambda),zeros(size(Y(:,indLambda))),'k-','LineWidth',2);

plot3(X(indN,:),Y(indN,:),max(max(PD))*.005+PD(indN,:),'ko','MarkerFaceColor','k','MarkerSize',8,'LineWidth',1.4);

plot3(X(indN,:),Y(indN,:),zeros(size(Y(indN,:))),'k-','LineWidth',2);

xlabel('\lambda','FontName','Times','FontSize',14);

ylabel('N','FontName','Times','FontSize',14);

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

r=axis; axis([r(1:5) .2]);

figure(2); hold on

aleph=trapz(PD(:,indLambda));

for n=1:length(N), plot([N(n) N(n)],[0 PD(n,indLambda)/aleph],'k:','LineWidth',1.7); end

plot(N,PD(:,indLambda)/aleph,'ko','MarkerFaceColor','w','MarkerSize',8,'LineWidth',1.9);

plot(Hn,PD(indN,indLambda)/aleph,'ko','MarkerFaceColor','k','MarkerSize',8,'LineWidth',1.9);

xlabel('N','FontName','Times','FontSize',14);

r=axis; mprob=r(4);

figure(3); hold on

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

for n=1:length(Lambda), plot([Lambda(n) Lambda(n)],[0 PD(indN,n)/aleph],'k:','LineWidth',1.7); end

plot(Lambda,PD(indN,:)/aleph,'ko','MarkerFaceColor','k','MarkerSize',8,'LineWidth',1.9);

xlabel('\lambda','FontSize',16);

r=axis;

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

axis([r(1:3) mprob]);

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