function [betamat,cmat,acceptmat]=gs_oprobit(X,y,k,n,nc,r1,r2,mu0,V0,...
    betadraw,eta0,kappa0,del2draw,Vc,stdc1,stdc2,c1draw,c2draw,...
    f0,f1,f2,f3,f4);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Data preparation for initial stuff that 
% (i) will not be updated in GS, or
% (ii) will eventually be updated, but is needed before that point
% in the sampler
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ydraw=y; % starting draws for latent data
c1count=0;
c2count=0;

R = r1+r2;

betamat=zeros(k,r2);
cmat=zeros(nc,r2);


thoucount=1000;
% start main loop
i=1;
for i=1:R
      
   %%%%%%%%%%%%
   % draw beta 
   %%%%%%%%%%%%
   V1=inv(inv(V0)+(1/del2draw)*(X'*X));
   mu1=V1*(inv(V0)*mu0+(1/del2draw)*(X'*ydraw));
   betadraw=mvnrnd(mu1,V1)';
   if i>r1
       betamat(:,i-r1)=betadraw/sqrt(del2draw); %since we're re-parameterizing
   end
 

   %%%%%%%%%%%%
   % draw del2
   %%%%%%%%%%%%
   eta1=(n+2*eta0)/2;
   kappa1=(1/2)*((ydraw-X*betadraw)'*(ydraw-X*betadraw)+2*kappa0);
   del2draw=1/gamrnd(eta1,1/kappa1);
   

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Draw thresholds and latent data
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
   % c1:  
   %%%%%%%%%%%%%%%%%%%%%%%%
   c1old=c1draw;
   
   % relevant stuff for y=2 cases
   X1=X(f1,:);
   A=normcdf((c1old-X1*betadraw)/sqrt(del2draw)); % will be f1 by 1
   B=normcdf((-X1*betadraw)/sqrt(del2draw));
   
   % relevant stuff for y=3 cases
   X2=X(f2,:);
   C=normcdf((c2draw-X2*betadraw)/sqrt(del2draw));
   D=normcdf((c1old-X2*betadraw)/sqrt(del2draw));
   
   logold=log(tnormpdf(0,sqrt(Vc),0,1,c1old))+sum(log(A-B)) + sum(log(C-D));
  
   c1can=tnormrnd(c1old,stdc1,0,c2draw);
   % relevant stuff for y=2 cases
   A=normcdf((c1can-X1*betadraw)/sqrt(del2draw)); % will f1 by 1
   B=normcdf((-X1*betadraw)/sqrt(del2draw));
   
   % relevant stuff for y=3 cases
   X2=X(f2,:);
   C=normcdf((c2draw-X2*betadraw)/sqrt(del2draw));
   D=normcdf((c1can-X2*betadraw)/sqrt(del2draw));
   
   logcan=log(tnormpdf(0,sqrt(Vc),0,1,c1can))+sum(log(A-B)) + sum(log(C-D));
   

   logacc = logcan -logold; 
  if log(rand)<logacc;
      c1draw = c1can;
      c1count = c1count +1;
  end



  % c2:  
   %%%%%%%%%%%%%%%%%%%%%%%%
   c2old=c2draw;
   
   % relevant stuff for y=3 cases
   X2=X(f2,:);
   A=normcdf((c2old-X2*betadraw)/sqrt(del2draw));
   B=normcdf((c1draw-X2*betadraw)/sqrt(del2draw));
   
   % relevant stuff for y=4 cases
   X3=X(f3,:);
   C=normcdf((1-X3*betadraw)/sqrt(del2draw));
   D=normcdf((c2old-X3*betadraw)/sqrt(del2draw));
   
   logold=log(tnormpdf(0,sqrt(Vc),0,1,c2old))+sum(log(A-B)) + sum(log(C-D));
  
   
   c2can=tnormrnd(c2old,stdc2,c1draw,1);
   
   % relevant stuff for y=3 cases
   X2=X(f2,:);
   A=normcdf((c2can-X2*betadraw)/sqrt(del2draw));
   B=normcdf((c1draw-X2*betadraw)/sqrt(del2draw));
   
   % relevant stuff for y=4 cases
   X3=X(f3,:);
   C=normcdf((1-X3*betadraw)/sqrt(del2draw));
   D=normcdf((c2can-X3*betadraw)/sqrt(del2draw));
   
   logcan=log(tnormpdf(0,sqrt(Vc),0,1,c2can))+sum(log(A-B)) + sum(log(C-D));
   
   
   logacc = logcan -logold; 
  if log(rand)<logacc;
      c2draw = c2can;
      c2count = c2count +1;
  end

  if i>r1
       cmat(:,i-r1)=[c1draw/sqrt(del2draw); c2draw/sqrt(del2draw) ;1/sqrt(del2draw)];
  end
  

  % latent data 
  
  ydraw(f0)=tnormrnd(X(f0,:)*betadraw,sqrt(del2draw),...
      -inf*ones(length(f0),1),zeros(length(f0),1));
  ydraw(f1)=tnormrnd(X(f1,:)*betadraw,sqrt(del2draw),...
      zeros(length(f1),1),c1draw*ones(length(f1),1));
  ydraw(f2)=tnormrnd(X(f2,:)*betadraw,sqrt(del2draw),...
      c1draw*ones(length(f2),1),c2draw*ones(length(f2),1));
  ydraw(f3)=tnormrnd(X(f3,:)*betadraw,sqrt(del2draw),...
      c2draw*ones(length(f3),1),ones(length(f3),1));
  ydraw(f4)=tnormrnd(X(f4,:)*betadraw,sqrt(del2draw),...
      ones(length(f4),1),inf*ones(length(f4),1));
     
  if i== thoucount
      i
      thoucount=thoucount+1000; 
  end
      
  
  acceptmat=[c1count/R c2count/R];     
    
end




