function [betamat]=gs_probit(X,y,k,n,r1,r2,mu0,V0,betadraw,f0,f1);

ydraw=y; % starting draws for latent data
R = r1+r2;
betamat=zeros(k,r2);


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

   
  %%%%%%%%%%%%%
  % latent data 
  %%%%%%%%%%%%%
  
  ydraw(f0)= tnormrnd(X(f0,:)*betadraw,1,...
      -inf*ones(length(f0),1),zeros(length(f0),1));
  %next is only difference to Tobit:
  ydraw(f1)= tnormrnd(X(f1,:)*betadraw,1,...
      zeros(length(f1),1),inf*ones(length(f1),1));
     
  if i== thoucount
      i
      thoucount=thoucount+1000; 
  end
       
    
end




