function [betafmat,betarmat,Emat,sig2mat,betarimat]...
    =gs_HNRM(y,Xf,Xr,Ximat,yimat,k,kf,kr,n,J,obs,r1,r2,muf0,Vf0,...
    mur0,Vr0,betardraw,v0,S0,Edraw,eta0,tau0,sig2draw);

R = r1+r2;
betafmat=zeros(kf,r2); 
betarmat=zeros(kr,r2);
betarimat=zeros(kr*n,1);
Emat=zeros(kr*(kr+1)/2,r2);
sig2mat=zeros(1,r2);

betaridraw=mvnrnd(betardraw,Edraw,n)'; %kr by n
% just to get things started

thoucount=1000;
% start main loop
i=1;
for i=1:R
        
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % draw fixed coefficients, conditional on betari 
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   int=reshape(betaridraw',n*kr,1)';
   int=repmat(int,J,1);
   brmat=reshape(int,n*J,kr);
   ywig=y-sum(Xr.*brmat,2);
   V1=inv(inv(Vf0)+(1/sig2draw)*Xf'*Xf);
   mu1=V1*(inv(Vf0)*muf0+(1/sig2draw)*Xf'*ywig);
   betafdraw=mvnrnd(mu1,V1)';
   if i>r1
       betafmat(:,i-r1)=betafdraw; 
   end
 
   %%%%%%%%%%%%%%%
   % draw betari's
   %%%%%%%%%%%%%%%
   betaridraw=[];
   for j=1:n;
        yi=yimat{j};
        Xi=Ximat{j};
        Xif=Xi(:,1:kf);
        Xir=Xi(:,kf+1:end);
        V1=inv(inv(Edraw)+(1/sig2draw)*Xir'*Xir);
        mu1=V1*(inv(Edraw)*betardraw+(1/sig2draw)*Xir'*(yi-Xif*betafdraw));
        draw=mvnrnd(mu1,V1)';
        betaridraw=[betaridraw draw];%will be kr by n
        % capture this for draws of sig2 below and Chib method
   end
  if i>r1
       betarimat(:,i-r1)=reshape(betaridraw,kr*n,1);% just needed for Chib routine
  end
      
  
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % draw betar (cond. on betari); 
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   ywig=y-Xf*betafdraw;
   V1=inv(inv(Vr0)+n*inv(Edraw));
   mu1=V1*(inv(Vr0)*mur0+inv(Edraw)*sum(betaridraw,2));
   betardraw=mvnrnd(mu1,V1)';
   if i>r1
       betarmat(:,i-r1)=betardraw; 
   end
    %%%%%%%%%%%%%%%%%%
    % draw sig2
    %%%%%%%%%%%%%%%%%%
    int=reshape(betaridraw',n*kr,1)';
    int=repmat(int,J,1);
    int=reshape(int,n*J,kr);
    
    betamat=[repmat(betafdraw',obs,1) int];
    int=sum(cell2mat(Ximat).*betamat,2);% Xf*bf + Xr*bri, stacked for all i
    
    eta1=(obs+2*eta0)/2;
    tau1=(1/2)*((y-int)'*(y-int)+2*tau0);
    sig2draw=1/gamrnd(eta1,1/tau1); 
    
    if i>r1
       sig2mat(:,i-r1)=sig2draw;
    end
    
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % draw E
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 v1=v0+n;
 int=betaridraw'-repmat(betardraw',n,1);%n by kr
 S1=S0+(int'*int); %equiv. to sum(bri-br)(bri-br)'
 Edraw=iwishrnd(S1,v1);
  
 if i>r1
      Emat(:,i-r1)=var_vec(Edraw);%collect unique elements of E as a vector
  end

    
    
        
    if i== thoucount
        i
        thoucount=thoucount+1000;
    end
        
    
 end


