function[logpmarg,lp,llf,lpost]=gs_Poisson_rwc_chib(y,X,k,N,mub0,Vb0,c,Vc,r1,r2,betamat);
betamean=mean(betamat,2); %k x 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% evaluate posterior
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nummat=zeros(r2,1); %will collect numerator stuff
denommat=zeros(r2,1);%will collect denominatory stuff
%%%%%%%%%%%%%%%%%%%%%
% numerator stuff
%%%%%%%%%%%%%%%%%%%%%%
% two components to evaluate:
% 1) the proposal density of betamean, conditional on some draw of beta
% 2) the acceptance probability of going from some beta to betamean
for i=1:r2
betacan=betamean;
betadraw=betamat(:,i);%pick a draw form original GS
lam=X*betacan; %n by 1
int=sum(y.*lam-exp(lam));
logcanbeta = log(mvnpdf(betacan,mub0,Vb0))+int;
lam=X*betadraw; %n by 1
int=sum(y.*lam-exp(lam));
logdrawbeta = log(mvnpdf(betadraw,mub0,Vb0))+int;
% log of aceptance probability
logacc = logcanbeta-logdrawbeta;
nummat(i)=(mvnpdf(betamean,betadraw,c*Vc))*exp(logacc);
end;
'numerator done'
%%%%%%%%%%%%%%%%%%%%%
% denominator stuff
%%%%%%%%%%%%%%%%%%%%%%
% draw new betas, given betamean
% for each draw, obtain and evaluate acceptance probability of moving from betamean to
% the new draw
betadraw=betamean;
for i=1:r2
betacan=mvnrnd(betadraw, c*Vc)'; %draw candidate beta
lam=X*betacan; %n by 1
int=sum(y.*lam-exp(lam));
logcanbeta = log(mvnpdf(betacan,mub0,Vb0))+int;
lam=X*betadraw; %n by 1
int=sum(y.*lam-exp(lam));
logdrawbeta = log(mvnpdf(betadraw,mub0,Vb0))+int;
% log of aceptance probability
logacc = logcanbeta-logdrawbeta;
denommat(i)=exp(logacc);
end
'denominator done'
%evaluate posterior
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lpost = log(mean(nummat)/mean(denommat));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%
% evaluate prior
%%%%%%%%%%%%%%%%%%%%%%%%
lp=log(mvnpdf(betamean',mub0',Vb0));
%%%%%%%%%%%%%%%%%%%%%%%
% evaluate log-LHF
%%%%%%%%%%%%%%%%%%%%%%%
mux=exp(X*betamean); %get n x 1 vector of "lambdas" for each observation
llf=sum(log(poisspdf(y,mux)));
logpmarg=lp+llf-lpost;
'chib routine complete'