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