function [logpmarg,lp,llf,lpost]...
=gs_sur_chib(X,Xcell,y,ycell,k1,k2,k3,k,M,n,r1,r2,mu0,...
V0,v0,S0,betamat,Emat);
Echibmat=zeros(r2,1);
betamean=mean(betamat,2);
Emean=vec_var(mean(Emat,2));
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
betadraw=betamat(:,i);
int=y-X*betadraw;%MxN by 1
int=reshape(int,M,n)';%n x M
S1=S0+(int'*int);
Echibmat(i)=iwishpdf(Emean,v1,S1);
if i== thoucount
i
thoucount=thoucount+1000;
end
end; %end main loop
'Chib round 1 done'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Round 2: chib stuff for beta
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%use Emean instead of Edraw
int1=0;
int2=0;
for j=1:n
Xi=cell2mat(Xcell(j)); % note: round brackets needed for this
yi=cell2mat(ycell(j));
int1=int1+Xi'*inv(Emean)*Xi;
int2=int2+Xi'*inv(Emean)*yi;
end
V1=inv(inv(V0)+int1);
mu1=V1*(inv(V0)*mu0+int2);
betachibmat= mvnpdf(betamean,mu1,V1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Evaluate and add log-posterior
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lpost=log(mean(Echibmat))+log(mean(betachibmat));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% evaluate the log-likelihood function
% at betamean, Emean
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
llf=0;
j=1;% k - counter
for j=1:n;
Xi=cell2mat(Xcell(j)); % note: round brackets needed for this
yi=cell2mat(ycell(j));
mui=Xi*betamean;
llfi=log(mvnpdf(yi,mui,Emean)); %individual LH component
llf=llf+llfi; % add up over individuals
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Evaluate and add log-priors
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lbeta=log(mvnpdf(betamean',mu0',V0));
lE=log(iwishpdf(Emean,v0,S0));
lp=lbeta+lE;
logpmarg=lp+llf-lpost;
'chib routine complete'