function out=simulate_markov(P,n,initial_state)
%SIMULATE_MARKOV takes a markov transition matrix P and generates a
%simulation of the chain for n steps.  It begins the process in the state
%given by initial_state
out=zeros(1,n);
out(1)=initial_state;
for i=2:n
    p=rand();
    tmp=cumsum(P(out(i-1),:));
    out(i)=find(p<tmp,1,'first');
end