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'


