function state = sample_markov_chain(initial_prob, trans, len, numex)%% Generate random sequences from a Markov chain.% state = sample_markov_chain(initial_prob, trans, len, numex)%% Each row of state is a sequence of length len.% The number of rows, numex, is 1 if omitted.if nargin==3  numex = 1;endstate = zeros(numex,len);S = length(initial_prob); % num statesR = rand(len*numex,1);rcount = 1;% Compute all the cdf's ahead of timeinit_cumprob = zeros(1,S+1);init_cumprob(1) = 0;for j=1:S  init_cumprob(j+1) = init_cumprob(j) + initial_prob(j);endcumprob = zeros(S, S+1);for i=1:S  cumprob(i,:) = [0 cumsum(trans(i,:))];endfor ex=1:numex  t = 1;    % j = sample(init_cumprob, R, rcount);  r = R(rcount);  rcount = rcount + 1;  for j=1:S    if (r >= init_cumprob(j)) & (r < init_cumprob(j+1))      break;    end  end    state(ex,t) = j;end  for ex=1:numex  for t=2:len    i = state(ex,t-1);    % j = sample(cumprob(i,:), R, rcount);    r = R(rcount);    rcount = rcount + 1;    for j=1:S      if (r >= cumprob(i,j)) & (r < cumprob(i,j+1))	break;      end    end        state(ex,t) = j;  endend