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);
for 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)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 lrr1
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)]