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