function [betamat,gmat,Emat,cmat,gimat,acceptmat]...
    =gs_HOP_normal(y,Ximat,yimat,Xgmat,ygmat,k,kf,kr,P,G,Gn,N,obs,nc,r1,r2,mub0,Vb0,...
    mug0,Vg0,gdraw,phi0,gam0,Edraw,eta0,kappa0,del2draw,...
    stdc1,stdc2,stdc3,stdc4,stdc5,stdc6,Vc,f1,f2,f3,f4,f5,f6,f7,f8,f9,...
    l1,l2,l3,l4,l5,l6,l7,l8,l9,c1draw,c2draw,c3draw,c4draw,c5draw,c6draw);

ydraw=y; % starting draws for latent data
c1count=0;
c2count=0;
c3count=0;
c4count=0;
c5count=0;
c6count=0;


R = r1+r2;

betamat=zeros(kf,r2);
gmat=zeros(kr,r2);
Emat=zeros(kr,r2);
cmat=zeros(nc,r2);
gimat=zeros(kr*N,r2);


thoucount=1000;
% start main loop
i=1;
for i=1:R
      
   %%%%%%%%%%%%%%%%%%%%%%%%%%%
   % 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=del2draw*eye(P)+Xir*(Edraw)*Xir';
        int1=int1+((Xif'/Vi)*Xif);
        int2=int2+((Xif'/Vi)*(yi-Xir*gdraw));
   end
    V1=inv(inv(Vb0)+int1);
    mu1=V1*(inv(Vb0)*mub0+int2);
    betadraw=mvnrnd(mu1,V1)';
   if i>r1
       betamat(:,i-r1)=betadraw/sqrt(del2draw); %since we're re-parameterizing
   end
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % draw g (uncond. on gamg and gamgi); 
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   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=del2draw*eye(P)+Xir*(Edraw)*Xir';
        int1=int1+((Xir'/Vi)*Xir);
        int2=int2+((Xir'/Vi)*(yi-Xif*betadraw));
    end
    V1=inv(inv(Vg0)+int1);
    mu1=V1*(inv(Vg0)*mug0+int2);
    gdraw=mvnrnd(mu1,V1)';
   if i>r1
       gmat(:,i-r1)=gdraw/sqrt(del2draw); %since we're re-parameterizing
   end
      
   %%%%%%%%%%%%%%%
   % draw gamgi's
   %%%%%%%%%%%%%%%
   gamgimat=[];
   Xf=[];
   Tg=[];
   for j=1:N;
        yi=yimat{j};
        Xi=Ximat{j};
        Xif=Xi(:,1:kf);
        Xir=Xi(:,kf+1:end);
        V1=inv(inv(Edraw)+(1/del2draw)*(Xir'*Xir));
        mu1=V1*(inv(Edraw)*gdraw+(1/del2draw)*Xir'*(yi-Xif*betadraw));
        draw=mvnrnd(mu1,V1)';
        gamgimat=[gamgimat draw];% will be kr by N
        %the following is useful for draws of del2 below
        Xf=[Xf;Xif]; 
        intdel2=Xir*draw;
        Tg=[Tg;intdel2];
   end
     if i>r1
       gimat(:,i-r1)=reshape(gamgimat,kr*N,1);% no need to reparameterize
       % just needed for Chib routine
  end

   %%%%%%%%%%%%
   % draw del2
   %%%%%%%%%%%%
   eta1=(obs+2*eta0)/2;
   kappa1=(1/2)*((ydraw-Xf*betadraw-Tg)'*(ydraw-Xf*betadraw-Tg)+2*kappa0);
   del2draw=1/gamrnd(eta1,1/kappa1);
   
   
   %%%%%%%%%
   % draw E
   %%%%%%%%%
   int=(gamgimat'-repmat(gdraw',N,1));%N by kr
   phi1=(N+2*phi0)/2;
   for ii=1:kr
       gam1=(1/2)*(int(:,ii)'*int(:,ii)+2*gam0);
       Edraw(ii,ii)=1/gamrnd(phi1,1/gam1);
   end
   
   if i>r1
       Emat(:,i-r1)= diag(Edraw/del2draw);
   end
   
    
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Draw thresholds and latent data
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   int=Xf*betadraw+Tg;%n by 1
   int1=int(f1);
   int2=int(f2);
   int3=int(f3);
   int4=int(f4);
   int5=int(f5);
   int6=int(f6);
   int7=int(f7);
   int8=int(f8);
   int9=int(f9);
   
      % c1:  
   %%%%%%%%%%%%%%%%%%%%%%%%
   c1old=c1draw;
   c1can=tnormrnd(c1old,stdc1,0,c2draw);
   
   % relevant stuff for y=2 cases
   A=normcdf((c1old-int2)/sqrt(del2draw)); % will f1 by 1
   B=normcdf((-int2)/sqrt(del2draw));
   % relevant stuff for y=3 cases
   C=normcdf((c2draw-int3)/sqrt(del2draw));
   D=normcdf((c1old-int3)/sqrt(del2draw));
   logold=log(tnormpdf(0,sqrt(Vc),0,1,c1old))+sum(log(A-B)) + sum(log(C-D))...
       -log(tnormpdf(c1can,stdc1,0,c2draw,c1old));
  
   % relevant stuff for y=2 cases
   A=normcdf((c1can-int2)/sqrt(del2draw)); % will f1 by 1
   B=normcdf((-int2)/sqrt(del2draw));
   % relevant stuff for y=3 cases
   C=normcdf((c2draw-int3)/sqrt(del2draw));
   D=normcdf((c1can-int3)/sqrt(del2draw));
   logcan=log(tnormpdf(0,sqrt(Vc),0,1,c1can))+sum(log(A-B)) + sum(log(C-D))...
       -log(tnormpdf(c1old,stdc1,0,c2draw,c1can));
   
   logacc = logcan -logold; 
  if log(rand)<logacc;
      c1draw = c1can;
      c1count = c1count +1;
  end

   % c2:  
   %%%%%%%%%%%%%%%%%%%%%%%%
   c2old=c2draw;
   c2can=tnormrnd(c2old,stdc2,c1draw,c3draw);
   
   % relevant stuff for y=3 cases
   A=normcdf((c2old-int3)/sqrt(del2draw)); % will f1 by 1
   B=normcdf((c1draw-int3)/sqrt(del2draw));
   % relevant stuff for y=4 cases
   C=normcdf((c3draw-int4)/sqrt(del2draw));
   D=normcdf((c2old-int4)/sqrt(del2draw));
   logold=log(tnormpdf(0,sqrt(Vc),0,1,c2old))+sum(log(A-B)) + sum(log(C-D))...
       -log(tnormpdf(c2can,stdc1,c1draw,c3draw,c2old));
  
   % relevant stuff for y=3 cases
   A=normcdf((c2can-int3)/sqrt(del2draw)); % will f1 by 1
   B=normcdf((c1draw-int3)/sqrt(del2draw));
   % relevant stuff for y=4 cases
   C=normcdf((c3draw-int4)/sqrt(del2draw));
   D=normcdf((c2can-int4)/sqrt(del2draw));
   logcan=log(tnormpdf(0,sqrt(Vc),0,1,c2can))+sum(log(A-B)) + sum(log(C-D))...
       -log(tnormpdf(c2old,stdc1,c1draw,c3draw,c2can));
   
   
   logacc = logcan -logold; 
  if log(rand)<logacc;
      c2draw = c2can;
      c2count = c2count +1;
  end
  
   % c3:  
   %%%%%%%%%%%%%%%%%%%%%%%%
   c3old=c3draw;
   c3can=tnormrnd(c3old,stdc3,c2draw,c4draw);
   % relevant stuff for y=4 cases
   A=normcdf((c3old-int4)/sqrt(del2draw)); % will f1 by 1
   B=normcdf((c2draw-int4)/sqrt(del2draw));
   % relevant stuff for y=5 cases
   C=normcdf((c4draw-int5)/sqrt(del2draw));
   D=normcdf((c3old-int5)/sqrt(del2draw));
   logold=log(tnormpdf(0,sqrt(Vc),0,1,c3old))+sum(log(A-B)) + sum(log(C-D))...
       -log(tnormpdf(c3can,stdc1,c2draw,c4draw,c3old));
  
   
   % relevant stuff for y=4 cases
   A=normcdf((c3can-int4)/sqrt(del2draw)); % will f1 by 1
   B=normcdf((c2draw-int4)/sqrt(del2draw));
   % relevant stuff for y=5 cases
   C=normcdf((c4draw-int5)/sqrt(del2draw));
   D=normcdf((c3can-int5)/sqrt(del2draw));
   logcan=log(tnormpdf(0,sqrt(Vc),0,1,c3can))+sum(log(A-B)) + sum(log(C-D))...
       -log(tnormpdf(c3old,stdc1,c2draw,c4draw,c3can));
   
   logacc = logcan -logold; 
  if log(rand)<logacc;
      c3draw = c3can;
      c3count = c3count +1;
  end
  
   % c4:  
   %%%%%%%%%%%%%%%%%%%%%%%%
   c4old=c4draw;
   c4can=tnormrnd(c4old,stdc4,c3draw,c5draw);
   % relevant stuff for y=5 cases
   A=normcdf((c4old-int5)/sqrt(del2draw)); % will f1 by 1
   B=normcdf((c3draw-int5)/sqrt(del2draw));
   % relevant stuff for y=6 cases
   C=normcdf((c5draw-int6)/sqrt(del2draw));
   D=normcdf((c4old-int6)/sqrt(del2draw));
   logold=log(tnormpdf(0,sqrt(Vc),0,1,c4old))+sum(log(A-B)) + sum(log(C-D))...
       -log(tnormpdf(c4can,stdc1,c3draw,c5draw,c4old));
  
 
   % relevant stuff for y=5 cases
   A=normcdf((c4can-int5)/sqrt(del2draw)); % will f1 by 1
   B=normcdf((c3draw-int5)/sqrt(del2draw));
   % relevant stuff for y=6 cases
   C=normcdf((c5draw-int6)/sqrt(del2draw));
   D=normcdf((c4can-int6)/sqrt(del2draw));
   logcan=log(tnormpdf(0,sqrt(Vc),0,1,c4can))+sum(log(A-B)) + sum(log(C-D))...
       -log(tnormpdf(c4old,stdc1,c3draw,c5draw,c4can));
   
   logacc = logcan -logold; 
  if log(rand)<logacc;
      c4draw = c4can;
      c4count = c4count +1;
  end
  
   % c5:  
   %%%%%%%%%%%%%%%%%%%%%%%%
   c5old=c5draw;
   c5can=tnormrnd(c5old,stdc5,c4draw,c6draw);
   % relevant stuff for y=6 cases
   A=normcdf((c5old-int6)/sqrt(del2draw)); % will f1 by 1
   B=normcdf((c4draw-int6)/sqrt(del2draw));
   % relevant stuff for y=7 cases
   C=normcdf((c6draw-int7)/sqrt(del2draw));
   D=normcdf((c5old-int7)/sqrt(del2draw));
   logold=log(tnormpdf(0,sqrt(Vc),0,1,c5old))+sum(log(A-B)) + sum(log(C-D))...
       -log(tnormpdf(c5can,stdc1,c4draw,c6draw,c5old));
  
   
   % relevant stuff for y=6 cases
   A=normcdf((c5can-int6)/sqrt(del2draw)); % will f1 by 1
   B=normcdf((c4draw-int6)/sqrt(del2draw));
   % relevant stuff for y=7 cases
   C=normcdf((c6draw-int7)/sqrt(del2draw));
   D=normcdf((c5can-int7)/sqrt(del2draw));
   logcan=log(tnormpdf(0,sqrt(Vc),0,1,c5can))+sum(log(A-B)) + sum(log(C-D))...
       -log(tnormpdf(c5old,stdc1,c4draw,c6draw,c5can));
   
   logacc = logcan -logold; 
  if log(rand)<logacc;
      c5draw = c5can;
      c5count = c5count +1;
  end
  
   % c6:  
   %%%%%%%%%%%%%%%%%%%%%%%%
   c6old=c6draw;
   c6can=tnormrnd(c6old,stdc6,c5draw,1);
   % relevant stuff for y=7 cases
   A=normcdf((c6old-int7)/sqrt(del2draw)); % will f1 by 1
   B=normcdf((c5draw-int7)/sqrt(del2draw));
   % relevant stuff for y=8 cases
   C=normcdf((1-int8)/sqrt(del2draw));
   D=normcdf((c6old-int8)/sqrt(del2draw));
   logold=log(tnormpdf(0,sqrt(Vc),0,1,c6old))+sum(log(A-B)) + sum(log(C-D))...
       -log(tnormpdf(c6can,stdc1,c5draw,1,c6old));
  
   
   % relevant stuff for y=7 cases
   A=normcdf((c6can-int7)/sqrt(del2draw)); % will f1 by 1
   B=normcdf((c5draw-int7)/sqrt(del2draw));
   % relevant stuff for y=8 cases
   C=normcdf((1-int8)/sqrt(del2draw));
   D=normcdf((c6can-int8)/sqrt(del2draw));
   logcan=log(tnormpdf(0,sqrt(Vc),0,1,c6can))+sum(log(A-B)) + sum(log(C-D))...
       -log(tnormpdf(c6old,stdc1,c5draw,1,c6can));
   
   
   logacc = logcan -logold; 
  if log(rand)<logacc;
      c6draw = c6can;
      c6count = c6count +1;
  end

  if i>r1
       cmat(:,i-r1)=(1/sqrt(del2draw))*[c1draw;c2draw;c3draw;c4draw;...
           c5draw;c6draw;1];
  end
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % latent data 
  
  ydraw(f1)=tnormrnd(int1,sqrt(del2draw),...
      -inf*ones(l1,1),zeros(l1,1));
  ydraw(f2)=tnormrnd(int2,sqrt(del2draw),...
      zeros(l2,1),c1draw*ones(l2,1));
  ydraw(f3)=tnormrnd(int3,sqrt(del2draw),...
      c1draw*ones(l3,1),c2draw*ones(l3,1));
  ydraw(f4)=tnormrnd(int4,sqrt(del2draw),...
      c2draw*ones(l4,1),c3draw*ones(l4,1));
  ydraw(f5)=tnormrnd(int5,sqrt(del2draw),...
      c3draw*ones(l5,1),c4draw*ones(l5,1));
  ydraw(f6)=tnormrnd(int6,sqrt(del2draw),...
      c4draw*ones(l6,1),c5draw*ones(l6,1));
  ydraw(f7)=tnormrnd(int7,sqrt(del2draw),...
      c5draw*ones(l7,1),c6draw*ones(l7,1));
  ydraw(f8)=tnormrnd(int8,sqrt(del2draw),...
      c6draw*ones(l8,1),ones(l8,1));
  ydraw(f9)=tnormrnd(int9,sqrt(del2draw),...
      ones(l9,1),inf*ones(l9,1));
    
  % rejuvenate all elements associated with y:
  yimat=mat2cell(ydraw,P*ones(N,1),1);
  ygmat=mat2cell(ydraw,(P*Gn)*ones(G,1),1);
  
     
  if i== thoucount
      i
      thoucount=thoucount+1000; 
  end
      
  acceptmat=(1/r2)*[c1count;c2count;c3count;c4count;c5count;c6count];
    
end




