function val=poigam(th,data)

% computes log posterior density of (theta_1, theta_2) for 
% Poisson/gamma exchangeable model
% theta_1 = log(alpha), theta_2 = log(mu)
% data = [z(i) e(i)]

z=data(:,1); e=data(:,2);
K=length(z);
zo=.53;  % hyperparameter value used in Albert paper

th1=th(:,1); th2=th(:,2);
al=exp(th1); mu=exp(th2);	

% compute likelihood
val=-K*gammaln(al)+K*al.*log(al./mu);
for i=1:K
   val=val+gammaln(al+z(i))-(al+z(i)).*log(al./mu+e(i));
end

% add prior and jacobian
val=val-2*log(al+zo)+th1;