function [beta1mat,beta2mat,beta3mat,sig11mat,sig22mat,sig33mat,cpmat,Zmat]=...
    gs_fmrm(X,y,k,n,G,r1,r2,mu0,V0,betadraw,v0,tau0,sig2draw,...
    acp0,cpdraw,Zdraw);

R = r1+r2;
beta1amat=zeros(k,r2); %will collect draws of beta1;
beta2amat=zeros(k,r2); %will collect draws of beta2;
beta3amat=zeros(k,r2); %will collect draws of beta3;

sig11mat=zeros(1,r2); %will collect draws of sig1;
sig22mat=zeros(1,r2); %will collect draws of sig2;
sig33mat=zeros(1,r2); %will collect draws of sig3;

cpmat=zeros(G,r2); %will collect draws of regime probabilities
Zmat=zeros(n,G);%will collect draws of regime assignments

sig11draw=sig2draw;
sig22draw=sig2draw;
sig33draw=sig2draw;

% beta1draw=betadraw;
% beta2draw=betadraw;
% beta3draw=betadraw;

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

    f1=find(Zdraw(:,1)==1); %flag obs. that fall into 1st regime
    f2=find(Zdraw(:,2)==1); %flag obs. that fall into 2nd regime
    f3=find(Zdraw(:,3)==1); %flag obs. that fall into 3rd regime
    
    n1=length(f1); %number of obs. in regime 1
    n2=length(f2);
    n3=length(f3);


% start main loop
i=1;
for i=1:R
     
      
    %%%%%%%%%%%%%
    % draw betas 
    %%%%%%%%%%%%%
       
    X1=X(f1,:);
    y1=y(f1,:);
        
    V1=inv(inv(V0)+(1/sig11draw)*X1'*X1);
    mu1=V1*(inv(V0)*mu0 +(1/sig11draw)*X1'*y1);    
    beta1draw=mvnrnd(mu1,V1)';    
    if i>r1
       beta1mat(:,i-r1)= beta1draw;
    end
    
    X2=X(f2,:);
    y2=y(f2,:);
        
    V1=inv(inv(V0)+(1/sig22draw)*X2'*X2);
    mu1=V1*(inv(V0)*mu0 +(1/sig22draw)*X2'*y2);    
    beta2draw=mvnrnd(mu1,V1)';    
    if i>r1
       beta2mat(:,i-r1)= beta2draw;
    end 
    
    X3=X(f3,:);
    y3=y(f3,:);
        
    V1=inv(inv(V0)+(1/sig33draw)*X3'*X3);
    mu1=V1*(inv(V0)*mu0 +(1/sig33draw)*X3'*y3);    
    beta3draw=mvnrnd(mu1,V1)';    
    if i>r1
       beta3mat(:,i-r1)= beta3draw;
    end 
    
  
    
    %%%%%%%%%%%%%%%%%%
    % draw sig2's
    %%%%%%%%%%%%%%%%%%
    
    v1=(n1+2*v0)/2;
    tau1=(1/2)*((y1-X1*beta1draw)'*(y1-X1*beta1draw)+2*tau0);
    sig11draw=1/gamrnd(v1,1/tau1); 
    % Matlab defines the ig scale as 1/tau, thus the inversion for the last
    % term
        if i>r1
            sig11mat(:,i-r1)=sig11draw;
        end
    
    v1=(n2+2*v0)/2;
    tau1=(1/2)*((y2-X2*beta2draw)'*(y2-X2*beta2draw)+2*tau0);
    sig22draw=1/gamrnd(v1,1/tau1); 
    % Matlab defines the ig scale as 1/tau, thus the inversion for the last
    % term
        if i>r1
            sig22mat(:,i-r1)=sig22draw;
        end
    
    v1=(n3+2*v0)/2;
    tau1=(1/2)*((y3-X3*beta3draw)'*(y3-X3*beta3draw)+2*tau0);
    sig33draw=1/gamrnd(v1,1/tau1); 
    % Matlab defines the ig scale as 1/tau, thus the inversion for the last
    % term
        if i>r1
            sig33mat(:,i-r1)=sig33draw;
        end
    
        
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Draw Z matrix (regime assignments)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % we will do this without a loop
    
    % compute vectors of multinomial probabilities for all observations
    int1=cpdraw(1)*normpdf(y,X*beta1draw,sqrt(sig11draw)); %n x 1, careful - normpdf requires std, not var!
    int2=cpdraw(2)*normpdf(y,X*beta2draw,sqrt(sig22draw)); %n x 1
    int3=cpdraw(3)*normpdf(y,X*beta3draw,sqrt(sig33draw)); %n x 1
    intmat=[int1 int2 int3];%n x G
    
    intsum=sum(intmat,2); % get denominator term for each obs.
    probmat=intmat./(repmat(intsum,1,G));
    
    Zdraw=mnldraws(probmat); %n by G
    
    %capture Zdraw to examine observation-specific assigment probabilities
        if i>r1
            Zmat=Zmat+Zdraw; %just keep adding up how often each obs. fell into a specific regime.
        end
    
        % IMPORTANT - update regime sized
    f1=find(Zdraw(:,1)==1); %flag obs. that fall into 1st regime
    f2=find(Zdraw(:,2)==1); %flag obs. that fall into 2nd regime
    f3=find(Zdraw(:,3)==1); %flag obs. that fall into 3rd regime
  
    n1=length(f1); %number of obs. in regime 1
    n2=length(f2); %number of obs. in regime 2
    n3=length(f3); %number of obs. in regime 1
        
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Draw regime probabilities
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
    acp1=acp0+[n1;n2;n3];
    cpdraw=dirrnd(acp1,1);
    
        if i>r1
            cpmat(:,i-r1)=cpdraw;
        end
        
    if i== thoucount
        i
        thoucount=thoucount+1000;
    end
        
    
 end


