function [LL_all, G_sum, H] = HP_gi_mle(x,Xfi,Xri,yi,k1,k2,k,gdraw,Edraw,betadraw)

beta=betadraw;
gidraw=x;
lam=Xfi*beta + Xri*gidraw; %J by 1
int=yi.*lam-exp(lam);
int=sum(int);

LL_all = -(log(mvnpdf(gidraw,gdraw,Edraw))+int);

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

if nargout>1
   int = sum(repmat((yi-exp(lam)),1,k2).*Xri)';%
   G_sum=-(-inv(Edraw)*(gidraw-gdraw) + int); % fminunc reverts signs
end

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

