function [LL_all, G_sum, H] =...
    HP_beta_mle(x,Xf,Xr,y,k1,k2,k,mub0,Vb0,Edraw,gidraw,N,J)

int=mat2cell(gidraw,ones(N,1),k2);
intmat=cell2mat(cellfun(@(x) repmat(x,J,1),int,'UniformOutput',0));%replicate RCs

beta=x;
lam=Xf*beta + sum(Xr.*intmat,2); %n by 1
int=y.*lam-exp(lam);
int=squeeze(reshape(int',1,J,N))';% N by J like Scott would use; works fine
int=sum(int,2); % sum over sites for each person
int=sum(int);

LL_all = -(log(mvnpdf(beta,mub0,Vb0))+int);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Gradient
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargout>1
   int = sum(repmat((y-exp(lam)),1,k1).*Xf)';% k1 by 1
   G_sum=-(-inv(Vb0)*(beta-mub0)+int); % fminunc reverts signs
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Hessian
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
if nargout >2;
   int = -(repmat(exp(lam),1,k1).*Xf)'*Xf;% k1 by k1
   H = -(-inv(Vb0) + int);% fminunc reverts signs
 end