function[betamat,acceptmat]=gs_Poisson_ic(y,X,k,N,mub0,Vb0,tau,v,betadraw,r1,r2);

R = r1+r2;
thoucount=1000;
betamat=zeros(k,r2);
betacount=0;

% start main loop
for i=1:R
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% draw beta
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% MLE routine to get mean and VCOV for proposal function
x0 = betadraw;

options=optimset('Algorithm','trust-region','display','off','TolFun',0.0001,'TolX',0.0001,...
     'GradObj','on', 'Hessian','on','DerivativeCheck','off');


[x,fval,exitflag,output,G_sum,H]=fminunc('Poisson_beta_mle',x0,options,X,y,k,mub0,Vb0);
tmean = x;
tvcov = tau*inv(H); % here enters the tuning parameter tau; note: Hessian coming
% out of fminunc is already pos. def.

%take a candidate draw
%Use Matlab's approach to draw from the mvt
[Rho,stdV] = corrcov(tvcov);
betacan=(tmean'+mvtrnd(Rho,v)*diag(stdV))'; 

lam=X*betacan; %n by 1
int=sum(y.*lam-exp(lam));
logcanbeta = log(mvnpdf(betacan,mub0,Vb0))+int - log(mvtpdf_klaus(betacan,tmean,tvcov,v));

lam=X*betadraw; %n by 1
int=sum(y.*lam-exp(lam));
logdrawbeta = log(mvnpdf(betadraw,mub0,Vb0))+int - log(mvtpdf_klaus(betadraw,tmean,tvcov,v));

% log of aceptance probability
logacc = logcanbeta-logdrawbeta;

% acceptance rule
lr=log(rand);
if lr<logacc;
    betadraw=betacan;
    betacount = betacount +1;
end

if i>r1
    betamat(:,i-r1)=betadraw; %since we're re-parameterizing
end
     
if i== thoucount
    i
    thoucount=thoucount+1000; 
end
      
  
end; %end main loop

acceptmat=betacount/R;     





