Looping for Markov Chain

3 ビュー (過去 30 日間)
Wendell
Wendell 2012 年 10 月 25 日
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)

回答 (1 件)

Kye Taylor
Kye Taylor 2012 年 10 月 25 日
編集済み: Kye Taylor 2012 年 10 月 25 日
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)')
  1 件のコメント
Wendell
Wendell 2012 年 10 月 26 日
Kye I appreciate very much the help...glad to have someone willing to share their knowledge and skill...

サインインしてコメントする。

カテゴリ

Help Center および File ExchangeMarkov Chain Models についてさらに検索

タグ

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by