function [bdraw,sigdraw,pdraw,spsum,spssq]=mcexGibbs(y,x,b0,sig0,p0,nit)
% function [bdraw,sigdraw,pdraw,spsum,spssq]=mcexGibbs(y,x,b0,sig0,p0,nit)
% makes nit complete Gibbs draws from b,sig,p, and cumulates over these draws the sum and
% the sum of squares of the state-probabilities (from which their mean and standard deviation
% can be computed).  This is just for the case where b is a column vector, not for cases where b
% varies with the state.  (mcbdraw would have to be modified as well to make that work.)
T=size(y,1);
bdraw=zeros(nit,size(b0,1));
sigdraw=zeros(nit,size(sig0,2));
pdraw= zeros(nit,size(p0,1));
spsum=zeros(T,1);
spssq=zeros(T,1);
H=[p0(1) 1-p0(2);1-p0(1) p0(2)];
b=[b0 b0];
sig=sig0;
for git=1:nit
   if mod(git,100)==0, disp(git),end
   [lh,sp]=mclh(y,x,sig,b,H);
   sdraw=mcdraw(H,sp);
   spsum=spsum+sp(:,1);
   spssq=spsum+sp(:,1).^2;
   b=mcbdraw(y,x,sig,sdraw);
   bdraw(git,:)=b';
   b=[b b];
   sig=mcsdraw(y,x,b,sdraw);
   sigdraw(git,:)=sig;
   H=mcHdraw(sdraw,H);
   pdraw(git,:)=diag(H)';
end
