function [betamat,sig2mat]=gs_tobit(X,y,k,n,r1,r2,mu0,V0,betadraw,v0,tau0,...
    sig2draw,f0,f1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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

R = r1+r2;

betamat=zeros(k,r2);
sig2mat=zeros(1,r2);
%%

thoucount=1000;
% start main loop
i=1;
for i=1:R
      
   %%%%%%%%%%%%
   % draw beta 
   %%%%%%%%%%%%
   V1=inv(inv(V0)+inv(sig2draw)*X'*X);
   mu1=V1*(inv(V0)*mu0+inv(sig2draw)*X'*ydraw);
   betadraw=mvnrnd(mu1,V1)';
   if i>r1
       betamat(:,i-r1)=betadraw;
   end
 
%%
   %%%%%%%%%%%%
   % draw sig2
   %%%%%%%%%%%%
   v1=(n+2*v0)/2;
   tau1=(1/2)*((ydraw-X*betadraw)'*(ydraw-X*betadraw)+2*tau0);
   sig2draw=1/gamrnd(v1,1/tau1);
   if i>r1
       sig2mat(:,i-r1)=sig2draw;
   end
   
  %%%%%%%%%%%%%
  % latent data - only needed for y=0 cases
  %%%%%%%%%%%%%
  
  ydraw(f0)= tnormrnd(X(f0,:)*betadraw,sqrt(sig2draw),...
      -inf*ones(length(f0),1),zeros(length(f0),1));
     
  if i== thoucount
      i
      thoucount=thoucount+1000; 
  end
       
    
end




