Looping for Markov Chain
1 view (last 30 days)
Show older comments
I need to run the following:
P = [0.942 0 0 0 0.058 0 0 0; 0 0.982 0 0 0.018 0 0 0; 0 0 0.999 0 0 0.001 0 0; 0 0 0 0.993 0 0.007 0 0; 0 0 0 0 0.987 0 0.013 0; 0 0 0 0 0 0.999 0.001 0; 0 0 0 0 0 0 0.972 0.028; 0 0 0 0 0 0 0 1]; x_0 = [1000; 800; 300; 400; 700; 200; 300; 0]; n = 100; for i = 1:n x = P*x_0; x(i) = P^(i)*x_0; end disp (x(i)) plot (x(i),n)
0 Comments
Answers (1)
Kye Taylor
on 25 Oct 2012
Edited: Kye Taylor
on 25 Oct 2012
You have all the pieces... its just awkward to organize (below I use cells) and display.. Try something like
n = 100;
P = cell(1,n);
P{1} = [0.942 0 0 0 0.058 0 0 0; 0 0.982 0 0 0.018 0 0 0; 0 0 0.999 0 0 0.001 0 0; 0 0 0 0.993 0 0.007 0 0; 0 0 0 0 0.987 0 0.013 0; 0 0 0 0 0 0.999 0.001 0; 0 0 0 0 0 0 0.972 0.028; 0 0 0 0 0 0 0 1];
x_0 = [1000; 800; 300; 400; 700; 200; 300; 0];
% I expect x_0 to be a probability
x_0 = x_0/sum(x_0);
% Create P, P^2, P^3,...
for i = 2:n
P{i} = P{1}*P{i-1};
end
g = @(A)A*x_0;
evolvedX = cellfun(g,P,'un',0); % applies P^k to x_0 for k = 1,2,3,...
plot([evolvedX{:}]')
xlabel('Time-step')
ylabel('Probability')
legend('x_0(1)','x_0(2)','x_0(3)','x_0(4)','x_0(5)','x_0(6)','x_0(7)','x_0(8)')
See Also
Categories
Find more on Markov Chain Models in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!