function [logpmarg,lp,llf,lpost]...
=gs_HNRM_nocovs_chib(y,Ximat,yimat,k,kf,kr,n,J,obs,r1,r2,muf0,Vf0,...
mur0,Vr0,phi0,gam0,eta0,tau0,betafmat,betarmat,betarimat,Emat,sig2mat);
betaristart=reshape(betarimat(:,1),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=diag(mean(Emat,2));
sig2mean=mean(sig2mat);
clear Emat;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Round 1: chib components for E
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
thoucount = 1000;
phi1=(n+2*phi0)/2;
for i=1:r2
betaridraw=reshape(betarimat(:,i),kr,n);
betardraw=betarmat(:,i);
int=betaridraw'-repmat(betardraw',n,1);%n by kr
mult=1;
for ii=1:kr
gam1=(1/2)*(int(:,ii)'*int(:,ii)+2*gam0);
mult=mult*invgampdf(Emean(ii,ii),phi1,gam1);
end
Echibmat(i)=mult;
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));
mult=1;
for ii=1:kr
mult=mult*invgampdf(Emean(ii,ii),phi0,gam0);
end
lE=log(mult);
lp=lbetaf+lbetar+lsig2+lE;
logpmarg=lp+llf-lpost;
'chib routine complete'