function [betamat,Emat]= gs_switch(X,X1,X2,X3,k1,k2,k3,k,m,n,r1,r2,mu0,V0,...
    betadraw,v0,S0,Edraw,f,g,nf,ng,l0,r0,lpos,rpos,ydraw)

R = r1+r2;
betamat=zeros(k,r2); 
Emat=zeros((m*(m+1)/2)-1,r2); %collects unique elements of E as a vector
% recall first element is one, so we can skip a row; 

huncount=100;
% start main loop
i=1;
for i=1:R
   
   %%%%%%%%%%%%%%%%%%%%%%
   % draw beta
   %%%%%%%%%%%%%%%%%%%%%
   
   int1=0; 
   int2=0;
   j=1;
   for j=1:n
       Xi=cell2mat(X(j)); % note: round brackets needed for this
       int1=int1+Xi'*inv(Edraw)*Xi;
       int2=int2+Xi'*inv(Edraw)*ydraw(j,:)';
   end
   V1=inv(inv(V0)+int1);
   mu1=V1*(inv(V0)*mu0+int2);
   betadraw=mvnrnd(mu1,V1)';
   if i>r1
       betamat(:,i-r1)=betadraw;
   end
 
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % draw E
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  v1=v0+n;
  yr= reshape(ydraw',3*n,1);%stacking of yi pairs - 3xn by 1
  Xr=cell2mat(X); %2xn by k
  int3=yr-Xr*betadraw;
  int3=reshape(int3,m,n)';%n x 3
  S1=S0+(int3'*int3); 
  Edraw=iw_sig11(v1,inv(S1),1);
  if i>r1
       int=var_vec(Edraw);
       Emat(:,i-r1)=int(2:end);
   end
 
  %%%%%%%%%%%%%%%%%%%
  % draw ystar
  %%%%%%%%%%%%%%%%%%%
  sig11=1; %"1", of course
  sig12=Edraw(2,1);
  sig13=Edraw(3,1);
  sig22=Edraw(2,2);
  sig23=Edraw(2,3);
  sig33=Edraw(3,3);
  
  rho12=sig12/(sqrt(sig11*sig22));
  rho13=sig13/(sqrt(sig11*sig33));
  
  E2r=[sig12;sig23];
  E3r=[sig13;sig23];
  E12=[sig11 sig12;sig12 sig22];
  E13=[sig11 sig13;sig13 sig33];
  
  beta1=betadraw(1:k1);
  beta2=betadraw(k1+1:k1+k2);
  beta3=betadraw(k1+k2+1:end);
  
  %Case 1: y1i=0
  %%%%%%%%%%%%%%%%%%%
   y2=ydraw(:,2); 
  % STEP 1: draw y1star for zero amounts, conditional on y2star (observed)
  mu=X1*beta1+(sig12/sig22)*(y2-X2*beta2);%conditional mean
  mu(g)=[]; % eliminate "1" cases, so now mu is length f
  yint=tnormrnd(mu,sqrt(sig11*(1-rho12^2)),l0,r0);
  ydraw(f,1)=yint;
  %STEP 2: draw y3star from its bi-conditional (on y1star,y2star) for "0"
  %cases
  y1=ydraw(:,1);
  Ec=E3r'*inv(E12);
  mu=X3*beta3+Ec(1)*(y1-(X1*beta1))+Ec(2)*(y2-(X2*beta2));
  Vcon=sig33-E3r'*inv(E12)*E3r;
  mu(g)=[];%eliminate "1" cases
  yint=normrnd(mu,sqrt(Vcon));
  ydraw(f,3)=yint;
  
  %Case 2: y1i=1
  %%%%%%%%%%%%%%%%%%%
  y3=ydraw(:,3);
  % STEP 1: draw y1star for nonzero amounts, conditional on y3star (observed)
  mu=X1*beta1+(sig13/sig33)*(y3-X3*beta3);%conditional mean
  mu(f)=[]; % eliminate "0" cases, so now mu is length f
  yint=tnormrnd(mu,sqrt(sig11*(1-rho13^2)),lpos,rpos);
  ydraw(g,1)=yint;
  %STEP 2: draw y2star from its bi-conditional (on y1star,y3star) for "1"
  %cases
  y1=ydraw(:,1);
  Ec=E2r'*inv(E13);
  mu=X2*beta2+Ec(1)*(y1-(X1*beta1))+Ec(2)*(y3-(X3*beta3));
  Vcon=sig22-E2r'*inv(E13)*E2r;
  mu(f)=[];%eliminate "0" cases
  yint=normrnd(mu,sqrt(Vcon));
  ydraw(g,2)=yint;

  
  % so the only elements of ydraw that are never updated are ydraw(f,2)and
  % ydraw(g,3);
     
  if i== huncount
      i
      huncount=huncount+100; 
  end
      
  
        
    
end


