function val=betabin(xy,data)
%
% (THIS IS THE CORRECTED VERSION)
% Computes log posterior of marginal density
% of hyperparameters in binomial/beta exchangeable model
%
% variables are log(m/(1-m)) and log(K)
% prior on (m, K) is 1/(K+1)^2
%
% data is [y n], where y are observed counts and n are sample sizes 

t1=xy(:,1);			% log(m/(1-m))
t2=xy(:,2);			% log(K)
b=exp(t2)./(1+exp(t1));
a=exp(t1+t2)./(1+exp(t1));

y=data(:,1); n=data(:,2);
k=length(y);

val=0*t1;

for i=1:k
	val=val+betaln(a+y(i),b+n(i)-y(i));
end
val=val-k*betaln(a,b);
val=val+t2-2*log(1+exp(t2));  % include jacobian and prior