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