function [betamat,sig2mat,gammat]=...
    gs_SSVS_fishing(X,y,n,k1,k2,k,r1,r2,mu0,Vx,V0,t2,c2,p,betadraw,v0,tau0,sig2draw);

R = r1+r2;
betamat=zeros(k,r2); %will collect draws of beta;
sig2mat=zeros(1,r2); %will collect draws of sig2;

thoucount=1000; 

% start main loop
i=1;
for i=1:R
        
    %%%%%%%%%%%%%
    % draw betas 
    %%%%%%%%%%%%%
    V1=inv(inv(V0)+(1/sig2draw)*X'*X);
    mu1=V1*(inv(V0)*mu0 +(1/sig2draw)*X'*y);    
    betadraw=mvnrnd(mu1,V1)';    
    if i>r1
       betamat(:,i-r1)= betadraw;
    end     
      
    %%%%%%%%%%%%%%%%%%
    % draw sig2
    %%%%%%%%%%%%%%%%%%
    v1=(n+2*v0)/2;
    tau1=(1/2)*((y-X*betadraw)'*(y-X*betadraw)+2*tau0);
    sig2draw=1/gamrnd(v1,1/tau1); 

    if i>r1
       sig2mat(:,i-r1)=sig2draw;
    end
    
    %%%%%%%%%%%%%%%%%%
    % draw gammas
    %%%%%%%%%%%%%%%%%%
    
%     Using a loop:
%     gamdraw=zeros(k-1,1);
%     for j=1:(k-1)
%         num=p*normpdf(betadraw(j+1),0,sqrt(c2*t2));
%         denom=num+(1-p)*normpdf(betadraw(j+1),0,sqrt(t2));
%         prob=num/denom;
%         gamdraw(j)=.5*sign(prob-rand)+0.5; %follow KPT code - nice & elegant!
%         % sign returns a -1 if prob<rand, else a 1, so this makes
%         % gamdraw(j) either 0 or 1
%     end

    % Without loop:
    num=p*normpdf(betadraw(k1+1:k),0,sqrt(c2*t2));
    denom=num+(1-p)*normpdf(betadraw(k1+1:k),0,sqrt(t2));
    prob=num./denom;
    gamdraw=binornd(1,prob); 
        
    if i>r1
        gammat(:,i-r1)=gamdraw; 
    end
    
    %KEY: update prior VCOV for beta
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    V0=diag([Vx;(1-gamdraw)*t2+gamdraw*c2*t2]); %k by k
    
    
    
        
    if i== thoucount
        i
        thoucount=thoucount+1000;
    end
        
    
 end



