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'


