function [betamat,sig2mat]= ...
    gs_normal_independent_keepall(X,y,n,k,r1,r2,mu0,V0,betadraw,v0,tau0,sig2draw);


% this version keeps all draws, i.e. does not discard the burn-ins

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

thoucount=1000; % shows every 1000th iteration in the command window
% you may want to change this to 100 for slower applications

% 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)';    

    betamat(:,i)= betadraw;
       
      
    %%%%%%%%%%%%%%%%%%
    % draw sig2
    %%%%%%%%%%%%%%%%%%
    v1=(n+2*v0)/2;
    tau1=(1/2)*((y-X*betadraw)'*(y-X*betadraw)+2*tau0);
    sig2draw=1/gamrnd(v1,1/tau1); 
    % Matlab defines the ig scale as 1/tau, thus the inversion for the last
    % term
    sig2mat(:,i)=sig2draw;
      
        
    if i== thoucount
        i
        thoucount=thoucount+1000;
    end
        
    
 end


