function [beta1mat,beta2mat,sig11mat,sig22mat,cpmat,Zmat]=...
gs_2cmm_differentMeans(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;
sig11mat=zeros(1,r2); %will collect draws of sig1;
sig22mat=zeros(1,r2); %will collect draws of sig2;
cpmat=zeros(G,r2); %will collect draws of regime probabilities
Zmat=zeros(n,G);%will collect draws of regime assignments
sig11draw=sig2draw;
sig22draw=sig2draw;
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
n1=length(f1); %number of obs. in regime 1
n2=length(f2);
% 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
%%%%%%%%%%%%%%%%%%
% 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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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)); %careful - normpdf requires std, not var!
int2=cpdraw(2)*normpdf(y,X*beta2draw,sqrt(sig22draw));
intmat=[int1 int2];
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
% update membership assigments and sizes (n1, n2)
f1=find(Zdraw(:,1)==1); %flag obs. that fall into 1st regime
f2=find(Zdraw(:,2)==1); %flag obs. that fall into 2nd regime
n1=length(f1); %number of obs. in regime 1
n2=length(f2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Draw regime probabilities
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
acp1=acp0+[n1;n2];
cpdraw=dirrnd(acp1,1);
if i>r1
cpmat(:,i-r1)=cpdraw;
end
if i== thoucount
i
thoucount=thoucount+1000;
end
end