function [LL_all, G_sum, H] = Poisson_beta_mle(x,X,y,k,mub0,Vb0)


beta=x;
lam=X*beta; %n by 1
int=sum(y.*lam-exp(lam));

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

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

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

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