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