% Markov chain emulator

m = input('Enter the Markov chain matrix (or its name) -> ');


% check that the matrix makes sense

[r,c] = size(m);

if (r~=c)
	error('Matrix must be square');
end;

if any(any(m<0))
	error('Matrix must be nonnegative');
end;

s = abs(sum(m')-1);
if any(s > 10*eps)
	error('Matrix rows must sum to 1')
end;

n = r; % n is the size of m

nits = input('Enter number of HUNDREDS of iterations -> ');


counts = zeros(1,n);
psums = cumsum(m')';


state = 1;

for k=1:nits
	rlist = rand(1,100); % grab 100 numbers
	for i=1:100
		% record this state
		counts(state) = counts(state)+1;
	
		% compute the next state
		state = min(find(rlist(i)<psums(state,:)));
	end;
	bar(counts/sum(counts));
	title([num2str(k*100),' iterations'])
	xlabel('State');
	ylabel('Frequency');
	grid;
	pause(0.5);
end;