function[amat,betamat,sig2mat,gammat,acceptcount]=...
    gs_MC3(X,y,k,r1,r2,gamstart,g)

R = r1+r2;
n=length(y);
ybar=mean(y);

betamat=zeros(k,r2); %will stack betas for each activity and iteration
sig2mat=zeros(1,r2); %will stack sig2s for each activity and iteration
gammat=zeros(k,r2); %will collect model vectors (0/1 indicators for 
% "variable is in model or not")
amat=zeros(1,r2);

gamdraw=gamstart;
thoucount=1000;
acceptcount=0;
% start main loop
Xg=X;
f=find(gamstart==0);
Xg(:,f)=[]; %eliminate corresponding columns in X

i=1;
for i=1:R
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % draw gamma (i.e. a model)
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
   %compute the loged kernel of the marginal likelihood for the old model
   gamold=gamdraw;
   Xo=Xg;
   ko=size(Xo,2);
   tauo=(1/2)*((g/(1+g))*(y'*y - (y'*Xo)*((Xo'*Xo)\(Xo'*y))) + ...
       (1/(1+g))*(y-ybar)'*(y-ybar)-(g/(1+g))*n*ybar^2);
   logmlold=(ko/2)*log(1/(1+g))-((n-1)/2)*log(tauo);
   
   %draw a new gamma
   int=gamold;
   u=randsample(k,1);%determine which position in gamma to change
   int(u)=abs(gamold(u)-1);%this changes 0 to 1 and vice versa
    
   if rand<1/(k+1)
       gamnew=gamold; %allow for the old model to be re-chosen w. prob 1/(k+1)
   else
       gamnew=int;
   end
   
   f=find(gamnew==0);
   Xn=X;
   Xn(:,f)=[];%eliminate corresponding columns in X
   kn=size(Xn,2);
   taun=(1/2)*((g/(1+g))*(y'*y - (y'*Xn)*((Xn'*Xn)\(Xn'*y))) + ...
       (1/(1+g))*(y-ybar)'*(y-ybar)-(g/(1+g))*n*ybar^2);
   logmlnew=(kn/2)*log(1/(1+g))-((n-1)/2)*log(taun);
 
   logacc = (logmlnew-logmlold);  

   if log(rand)<logacc;
       gamdraw = gamnew;
       Xg=Xn;
       acceptcount = acceptcount +1; 
   end
   
   if i>r1
        gammat(:,i-r1)=gamdraw; 
   end
   
   %Draw sig2's
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   tau1=(1/2)*((g/(1+g))*(y'*y - (y'*Xg)*((Xg'*Xg)\(Xg'*y))) + ...
       (1/(1+g))*(y-ybar)'*(y-ybar)-(g/(1+g))*n*ybar^2);
   v1=(n-1)/2;
   sig2draw=1./gamrnd(v1,1/tau1); 

   if i>r1
        sig2mat(:,i-r1)=sig2draw; 
   end
   
   %Draw beta
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   if max(gamdraw)>0
       mu1=(g/(g+1))*(Xg'*Xg)\(Xg'*y);
       V1=sig2draw*(g/(g+1))*inv(Xg'*Xg);
       betadraw=mvnrnd(mu1,V1)';

       %We need to fill in the zeros before collecting
       int=zeros(k,1);
       f=find(gamdraw==1);
       int(f)=betadraw;%this will place the betadraws in the correct positions
   else
       int=zeros(k,1);
   end
 
   if i>r1
        betamat(:,i-r1)=int; 
   end
   
      %Draw intercept (alpha)
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   adraw=ybar+sqrt(sig2draw)*randn;
 
   if i>r1
        amat(:,i-r1)=adraw; 
   end
   
   
  if i==thoucount
      thoucount=thoucount+1000;
      i
  end     

end