function [logpmarg,lp,llf,lpost]...
    =gs_HNRM_chib(y,Ximat,yimat,k,kf,kr,n,J,obs,r1,r2,muf0,Vf0,...
    mur0,Vr0,v0,S0,eta0,tau0,betafmat,betarmat,betarimat,Emat,sig2mat);

betaristart=reshape(betarimat(:,end),kr,n)';%n by kr, each row for a given individual
betafstart=betafmat(:,end);
betarstart=betarmat(:,end);
Estart=Emat(:,end);
sig2start=sig2mat(end);

betaridraw=betaristart;
betafdraw=betafstart;
betardraw=betarstart;
Edraw=Estart;
sig2draw=sig2start;

betafchibmat=zeros(r2,1);
betarchibmat=zeros(r2,1);
Echibmat=zeros(r2,1);
sig2chibmat=zeros(r2,1);

betafmean=mean(betafmat,2);
betarmean=mean(betarmat,2);
Emean=vec_var(mean(Emat,2));
sig2mean=mean(sig2mat);


clear Emat;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Round 1: chib components for E
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
thoucount = 1000;

v1=v0+n; %posterior DoF - keep out of loop for speed
for i=1:r2
    betaridraw=reshape(betarimat(:,i),kr,n);
    betardraw=betarmat(:,i);
    int=betaridraw'-repmat(betardraw',n,1);%n by kr
    S1=S0+(int'*int); %equiv. to sum(bri-br)(bri-br)'

    Echibmat(i)=iwishpdf(Emean,v1,S1); 
    if i== thoucount
        i
        thoucount=thoucount+1000; 
    end 
end; %end main loop
'Chib round 1 done'
clear betarimat betarmat; %save memory

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Round 2: chib stuff for betaf
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
i=1;
thoucount=1000;
% start main loop
betaridraw=betaristart;
betafdraw=betafstart;
betardraw=betarstart;
sig2draw=sig2start;

% use Emean instead of Edraw

% start main loop
i=1;
for i=1:r2
     
   %%%%%%%%%%%%%%%%%%%%%%%%%%%
   % draw fixed coefficients 
   %%%%%%%%%%%%%%%%%%%%%%%%%%%
   int1=zeros(kf,kf);
   int2=zeros(kf,1);
   for j=1:n; 
        yi=yimat{j};
        Xi=Ximat{j};
        Xif=Xi(:,1:kf);
        Xir=Xi(:,kf+1:end);
        Vi=sig2draw*eye(J)+Xir*(Emean)*Xir';
        int1=int1+(Xif'*inv(Vi)*Xif);
        int2=int2+(Xif'*inv(Vi)*(yi-Xir*betardraw));
   end
    V1=inv(inv(Vf0)+int1);
    mu1=V1*(inv(Vf0)*muf0+int2);
    betafdraw=mvnrnd(mu1,V1)';
    betafchibmat(i)= mvnpdf(betafmean,mu1,V1);

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % draw betar (uncond. on betari); 
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   int1=zeros(kr,kr);
   int2=zeros(kr,1);
   for j=1:n;
        yi=yimat{j};
        Xi=Ximat{j};
        Xif=Xi(:,1:kf);
        Xir=Xi(:,kf+1:end);
        Vi=sig2draw*eye(J)+Xir*(Emean)*Xir';
        int1=int1+(Xir'*inv(Vi)*Xir);
        int2=int2+(Xir'*inv(Vi)*(yi-Xif*betafdraw));
    end
    V1=inv(inv(Vr0)+int1);
    mu1=V1*(inv(Vr0)*mur0+int2);
    betardraw=mvnrnd(mu1,V1)';
   
   %%%%%%%%%%%%%%%
   % draw betari's
   %%%%%%%%%%%%%%%
   betaridraw=[];
   for j=1:n;
        yi=yimat{j};
        Xi=Ximat{j};
        Xif=Xi(:,1:kf);
        Xir=Xi(:,kf+1:end);
        V1=inv(inv(Emean)+(1/sig2draw)*Xir'*Xir);
        mu1=V1*(inv(Emean)*betardraw+(1/sig2draw)*Xir'*(yi-Xif*betafdraw));
        draw=mvnrnd(mu1,V1)';
        betaridraw=[betaridraw draw];%will be kr by n
        % capture this for draws of sig2 below and Chib method
   end
      
    %%%%%%%%%%%%%%%%%%
    % draw sig2
    %%%%%%%%%%%%%%%%%%
    int=reshape(betaridraw',n*kr,1)';
    int=repmat(int,J,1);
    int=reshape(int,n*J,kr);
    
    betamat=[repmat(betafdraw',obs,1) int];
    int=sum(cell2mat(Ximat).*betamat,2);% Xf*bf + Xr*bri, stacked for all i
    
    eta1=(obs+2*eta0)/2;
    tau1=(1/2)*((y-int)'*(y-int)+2*tau0);
    sig2draw=1/gamrnd(eta1,1/tau1); 
    
    
      if i== thoucount
      i
      thoucount=thoucount+1000; 
  end
     
end
'Chib Round 2 done'

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Round 3: chib stuff for betar
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
i=1;
thoucount=1000;
% start main loop
betaridraw=betaristart;
sig2draw=sig2start;
% use Emean instead of Edraw
% use betafmean instead of betafdraw

% start main loop
i=1;
for i=1:r2
     
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % draw betar (uncond. on betari); 
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   int1=zeros(kr,kr);
   int2=zeros(kr,1);
   for j=1:n;
        yi=yimat{j};
        Xi=Ximat{j};
        Xif=Xi(:,1:kf);
        Xir=Xi(:,kf+1:end);
        Vi=sig2draw*eye(J)+Xir*(Emean)*Xir';
        int1=int1+(Xir'*inv(Vi)*Xir);
        int2=int2+(Xir'*inv(Vi)*(yi-Xif*betafmean));
    end
    V1=inv(inv(Vr0)+int1);
    mu1=V1*(inv(Vr0)*mur0+int2);
    betardraw=mvnrnd(mu1,V1)';
    betarchibmat(i)= mvnpdf(betarmean,mu1,V1);
    
   %%%%%%%%%%%%%%%
   % draw betari's
   %%%%%%%%%%%%%%%
   betaridraw=[];
   for j=1:n;
        yi=yimat{j};
        Xi=Ximat{j};
        Xif=Xi(:,1:kf);
        Xir=Xi(:,kf+1:end);
        V1=inv(inv(Emean)+(1/sig2draw)*Xir'*Xir);
        mu1=V1*(inv(Emean)*betardraw+(1/sig2draw)*Xir'*(yi-Xif*betafmean));
        draw=mvnrnd(mu1,V1)';
        betaridraw=[betaridraw draw];
   end
      
    %%%%%%%%%%%%%%%%%%
    % draw sig2
    %%%%%%%%%%%%%%%%%%
    int=reshape(betaridraw',n*kr,1)';
    int=repmat(int,J,1);
    int=reshape(int,n*J,kr);
    
    betamat=[repmat(betafmean',obs,1) int];
    int=sum(cell2mat(Ximat).*betamat,2);% Xf*bf + Xr*bri, stacked for all i
    
    eta1=(obs+2*eta0)/2;
    tau1=(1/2)*((y-int)'*(y-int)+2*tau0);
    sig2draw=1/gamrnd(eta1,1/tau1);
    
      
      if i== thoucount
      i
      thoucount=thoucount+1000; 
  end
     
end
'Chib round 3 done'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Round 4: chib stuff for sig2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
i=1;
thoucount=1000;
% start main loop
betaridraw=betaristart;
sig2draw=sig2start;

% use Emean instead of Edraw
% use betafmean instead of betafdraw
% use betarmean instead of betardraw

% start main loop
i=1;
for i=1:r2
       
   %%%%%%%%%%%%%%%
   % draw betari's
   %%%%%%%%%%%%%%%
   betaridraw=[];
   for j=1:n;
        yi=yimat{j};
        Xi=Ximat{j};
        Xif=Xi(:,1:kf);
        Xir=Xi(:,kf+1:end);
        V1=inv(inv(Emean)+(1/sig2draw)*Xir'*Xir);
        mu1=V1*(inv(Emean)*betarmean+(1/sig2draw)*Xir'*(yi-Xif*betafmean));
        draw=mvnrnd(mu1,V1)';
        betaridraw=[betaridraw draw];
   end
      
    %%%%%%%%%%%%%%%%%%
    % draw sig2
    %%%%%%%%%%%%%%%%%%
    int=reshape(betaridraw',n*kr,1)';
    int=repmat(int,J,1);
    int=reshape(int,n*J,kr);
    
    betamat=[repmat(betafmean',obs,1) int];
    int=sum(cell2mat(Ximat).*betamat,2);% Xf*bf + Xr*bri, stacked for all i
    
    eta1=(obs+2*eta0)/2;
    tau1=(1/2)*((y-int)'*(y-int)+2*tau0);
    sig2draw=1/gamrnd(eta1,1/tau1); 
    sig2chibmat(i)=invgampdf(sig2mean,eta1,tau1);
    
      
      if i== thoucount
      i
      thoucount=thoucount+1000; 
  end
     
end
'Chib round 4 done'

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Evaluate and add log-posterior 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lpost=log(mean(Echibmat))+log(mean(betafchibmat))+...
      log(mean(betarchibmat))+log(mean(sig2chibmat));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% evaluate the log-likelihood function
% at betafmean, betarmean, Emean - UNCONDITIONAL of betari's
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

llf=0;
j=1;% k - counter

for j=1:n;
    Xi=cell2mat(Ximat(j)); 
    yi=cell2mat(yimat(j));
    Xif=Xi(:,1:kf);
    Xir=Xi(:,kf+1:end);
    Vi=sig2mean*eye(J)+Xir*(Emean)*Xir';
    mui=Xif*betafmean+Xir*betarmean;
    llfi=log(mvnpdf(yi,mui,Vi)); %individual LH component
    llf=llf+llfi; % add up over individuals
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Evaluate and add log-priors
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

lbetaf=log(mvnpdf(betafmean',muf0',Vf0));
lbetar=log(mvnpdf(betarmean',mur0',Vr0));
lsig2=log(invgampdf(sig2mean,eta0,tau0));
lE=log(iwishpdf(Emean,v0,S0));

lp=lbetaf+lbetar+lsig2+lE;

logpmarg=lp+llf-lpost;
'chib routine complete'
  
  
  
  
  