function [betamat,gmat,gimat,Emat,acceptmat]=gs_HPpar(y,ycell,X,Xf,Xr,Xfcell,Xrcell,k1,k2,k,N,J,n,...
    mub0,Vb0,taub,vb,betadraw,...
    mug0,Vg0,taug,vg,v0,S0,Edraw,...
    r1,r2);


R = r1+r2;
thoucount=1000;
betamat=zeros(k1,r2);
gmat=zeros(k2,r2);
Emat=zeros(k2*(k2+1)/2,r2);
gimat=zeros(N*k2,1);

betacount=0;
gcountmat=zeros(R,1);

gdraw=mvnrnd(mug0,Vg0)';
gidraw=mvnrnd(gdraw,Edraw,N);%N by 2


% start main loop
for i=1:R
i
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% draw gi's
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
j=1;
gcount=zeros(N,1);
parfor j=1:N
    giold=gidraw(j,:)'; %2 by 1
    Xfi=Xfcell{j};
    Xri=Xrcell{j};
    yi=ycell{j};
    x0=giold;
    options=optimset('LargeScale','on','display','off','TolFun',0.0001,'TolX',0.0001,...
    'GradObj','on', 'Hessian','on','DerivativeCheck','off');
    [x,fval,exitflag,output,G_sum,H]=fminunc('HP_gi_mle',x0,options,...
    Xfi,Xri,yi,k1,k2,k,gdraw,Edraw,betadraw);
    tmean = x; %
    tvcov = taug*inv(H); 
    
    %take a candidate draw
    %Use Matlab's approach to draw from the mvt
    [Rho,stdV] = corrcov(tvcov);
    gican=(tmean'+mvtrnd(Rho,vg)*diag(stdV))'; 
        
    Xb=Xfi*betadraw+Xri*gican;
    int=sum(yi.*Xb-exp(Xb));
    logcangi = log(mvnpdf(gican,gdraw,Edraw))+ int -...
    log(mvtpdf_klaus(gican,tmean,tvcov,vg));
    
    Xb=Xfi*betadraw+Xri*giold;
    int=sum(yi.*Xb-exp(Xb));
    logoldgi = log(mvnpdf(giold,gdraw,Edraw))+ int -...
    log(mvtpdf_klaus(giold,tmean,tvcov,vg));
    
    logacc=logcangi-logoldgi;
    if log(rand)<logacc;
        gidraw(j,:)=gican';%N by J
        gcount(j) = 1;
    end
end
gcountmat(i)=sum(gcount)/N; %take mean over all panels
if i>r1
    gimat(:,i-r1)=reshape(gidraw,N*k2,1); 
    %stack all N first RC's , then all N second RCs
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% draw beta
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% MLE routine to get mean and VCOV for proposal function
betaold=betadraw;
x0 = betadraw;
options=optimset('LargeScale','on','display','off','TolFun',0.0001,'TolX',0.0001,...
    'GradObj','on', 'Hessian','on','DerivativeCheck','off');

[x,fval,exitflag,output,G_sum,H]=fminunc('HP_beta_mle',x0,options,Xf,Xr,y,k1,k2,k,...
    mub0,Vb0,Edraw,gidraw,N,J);
tmean = x;
tvcov = taub*inv(H); 

%take a candidate draw
%Use Matlab's approach to draw from the mvt
[Rho,stdV] = corrcov(tvcov);
betacan=(tmean'+mvtrnd(Rho,vb)*diag(stdV))'; 

int=mat2cell(gidraw,ones(N,1),k2);
intmat=cell2mat(cellfun(@(x) repmat(x,J,1),int,'UniformOutput',0));%replicate RCs

Xb=Xf*betacan+sum(Xr.*intmat,2); %n by 1
int=y.*Xb-exp(Xb);
int=squeeze(reshape(int',1,J,N))';% N by J like Scott would use; works fine
int=sum(int,2); % sum over sites for each person
int=sum(int);
logcanbeta = log(mvnpdf(betacan,mub0,Vb0))+ int -...
    log(mvtpdf_klaus(betacan,tmean,tvcov,vb));

Xb=Xf*betaold+sum(Xr.*intmat,2); %n by 1
int=y.*Xb-exp(Xb);
int=squeeze(reshape(int',1,J,N))';% N by J like Scott would use; works fine
int=sum(int,2); % sum over sites for each person
int=sum(int);
logoldbeta = log(mvnpdf(betaold,mub0,Vb0))+ int -...
    log(mvtpdf_klaus(betaold,tmean,tvcov,vb));

% log of aceptance probability
logacc = logcanbeta-logoldbeta;

% acceptance rule
lr=log(rand);
if lr<logacc;
    betadraw=betacan;
    betacount = betacount +1;
end

if i>r1
    betamat(:,i-r1)=betadraw; 
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% draw g (mean of RCs)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
V1=inv(inv(Vg0)+N*inv(Edraw));
mu1=V1*(inv(Vg0)*mug0+inv(Edraw)*sum(gidraw',2));
gdraw=mvnrnd(mu1,V1)';% k2 by 1

if i>r1
    gmat(:,i-r1)=gdraw; 
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% draw E (VCOV of RCs)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

int=gidraw-repmat(gdraw',N,1);%Nxk2
v1=v0+N;
S1=S0+(int'*int);% k2 by k2
Edraw = iwishrnd(S1,v1);
 
Estore=var_vec(Edraw); % collect unique elements of Edraw into k2(k2+1)/2 by 1 vector


if i>r1
    Emat(:,i-r1)=Estore; 
end

     
if i== thoucount
    i
    thoucount=thoucount+1000; 
end




  
end; %end main loop

acceptmat=[betacount/R; mean(gcountmat)]  





