function H=mcHdraw(sdraw,oldH)
% function H=mcHdraw(sdraw,oldH)
% produces a Metropolis-Hastings draw from the Dirichlet, modified 
% by the initial steady-state density, that is the posterior on the elements of H.
% The draw is made from the Dirichlet itself, and then a retain/repeat decision is
% made using the initial steady-state distribution as weighting function, so that the
% the target distribution becomes the actual posterior.
T=size(sdraw,1);
ns=size(oldH,1);
H=zeros(ns,ns);
for it=2:T
   H(sdraw(it),sdraw(it-1))=H(sdraw(it),sdraw(it-1))+1;
end
for is=1:ns
   for js=1:ns
      H(is,js)=sum(randn(2*(H(is,js)+1),1).^2)/2;
   end
end
Htot=sum(H);
badcolindex=find(Htot==0);
if ~isempty(badcolindex)
   disp('at bad H')
end
H(:,badcolindex)=ones(ns,length(badcolindex))/ns;
H=H./repmat(sum(H),ns,1);
[v,d]=eig(H);
j1=find(abs(1-diag(d))<100*eps);
pbar=v(:,j1)./repmat(sum(v(:,j1)),ns,1);
if size(pbar,2)>1
   pbar = mean(pbar,2);% in case of non-ergodic H matrix.
end
[v,d]=eig(oldH);
j1=find(abs(1-diag(d))<eps*100);
oldpbar=v(:,j1)./sum(v(:,j1));
if size(oldpbar,2)>1
   oldpbar = mean(oldpbar,2);% in case of non-ergodic H matrix.
end
if oldpbar(sdraw(1))==0
   r=1.1;
else
   r=pbar(sdraw(1))/oldpbar(sdraw(1));
end
u=rand;
if u>r
   H=oldH;
end
