Program Is not Plotting Correctly

1 次查看（过去 30 天）
Jason
Jason 2023-10-22
评论： Walter Roberson 2023-10-23
Hello,
I have the code below that is calculating the state transition matrix and the state probability vector with given failure rates of l1 = 2×10-4 failures/hr and λ2 = 1×10-5 failures/hr and Δt = 0.1 hr. I am trying to generate a graph of the state probabilities but nothing is showing up on the plot.
Any help would be much appreciated!
function markov_chain(l1, l2, t) %input variables: lambda 1,2 and t
%defining the transition matrix from the given Markov Chain
P = zeros(6,6); %starting the 6x6 matrix
P(1,1) = (2*l2*t) + (2*l1*t)- (l1+l2)*t;
P(1,2) = (2*l1*t) - (l1 + 2*l2)*t;
P(2,1) = (2*l2*t) + (3*l1*t) - (2*l1+l2)*t;
P(2,2) = (3*l1*t) - (2*l1 + 2*l2)*t;
P(3,1) = (2*l2*t) - (3*l1 + l2) *t;
P(3,2) = (-3*l1 - 2*l2)*t;
P = P - diag(sum(P,2));
%F1 = (l1*t +l1*t);
%F2 = (l2*t + l2*t + l2*t;
%Steady State Probabilty Vector
% Compute EigenVector
% Find the EV corresponding to the '1' EV
% Normalize the eigenvector such that the sum of its elements equals 1,
% The transpose of the result is the steady state probability row-vector.
[X, Y] = eig(P');
steadyState = X(:,1)/sum(Y(:,1));
%State Probabilities
n = 100; %Stepping through 100 times
time = 0:t:(n*t); %time vector defined by input variable t as the increment
states = zeros(length(time), 6); % Initializing the state vector with all zeros for the length of time
states(1,:) = steadyState'; %filled with steady state prob vector
for i = 2:length(time)
states(i,:) = states (i-1,:) * expm(P*t); %state vector is filled by computing states * exponential matrix of P*t
end
%Plotting the probabilities over time
plot(time, states)
xlabel("time(s)")
ylabel("State Probability")
title("State Probability Time")
figure;
bar(steadyState);
xlabel('State');
ylabel('Steady-State Probability');
title('Steady-State Probability Vector');
end

请先登录，再回答此问题。

回答（1 个）

Walter Roberson
Walter Roberson 2023-10-22
If you are using plot() inside a loop but nothing is showing up on the plot, then:
  • make sure you use hold on
  • Include a marker in your plot() options, such as plot(x, y, '-*')
If you now get out a bunch of disconnected points, then the challenge is that plot() only joins finite points that are plotted in the same plot() call, so if you plot one point at a time then no line will be generated; the cure for that is to record all of the coordinates and plot() them at the end.
If you still do not get out anything on the graph, check nnz(~isfinite(x)), nnz(~isfinite(y)) [replacing x and y with your actual variables.] plot() cannot draw lines where there are infinite or nan values.
  2 个评论
Walter Roberson
Walter Roberson 2023-10-23
编辑：Walter Roberson 2023-10-23
Ran in:
Your system is numerically very sensitive, and somehow the first eigenvector is coming out all zero, which is leading to infinite and nan steady state, which destroys the rest of the calculations.
If you proceed with symbolic numbers you can get a plot -- but it is pretty dubious.
format long g
markov_chain((2e-4), (1e-5), (0.1))
P = -1.7999999999999997e-05 1.8e-05 0 0 0 0 2.0999999999999995e-05 -2.0999999999999995e-05 0 0 0 0 -5.9000000000000011e-05 -6.2000000000000016e-05 0.00012100000000000003 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 X = 0.75925660236529646 -0.70710678118654746 -0.40824829046386296 0 0 0 0.65079137345596849 0.70710678118654757 -0.40824829046386307 0 0 0 0 0 0.81649658092772592 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 Y = 0 0 0 0 0 0 0 -3.8999999999999999e-05 0 0 0 0 0 0 0.00012100000000000003 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
steadyState = 6×1
Inf Inf NaN NaN NaN NaN
ans =
6
ans =
0
ans =
0
ans =
606
markov_chain(sym(2e-4), sym(1e-5), sym(0.1))
P = -1.8e-05 1.8e-05 0 0 0 0 2.0999999999999999e-05 -2.0999999999999999e-05 0 0 0 0 -5.8999999999999998e-05 -6.2000000000000003e-05 0.000121 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 X = 0.75925660236529646 -0.70710678118654757 -0.40824829046386296 0 0 0 0.65079137345596849 0.70710678118654746 -0.40824829046386307 0 0 0 0 0 0.81649658092772592 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 Y = -3.3881317890172014e-21 0 0 0 0 0 0 -3.8999999999999999e-05 0 0 0 0 0 0 0.000121 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
steadyState = 6×1
1.0e+00 * -2.2409299568171e+20 -1.92079710584323e+20 0 0 0 0
ans =
0
ans =
0
ans =
0
ans =
0
function markov_chain(l1, l2, t) %input variables: lambda 1,2 and t
%defining the transition matrix from the given Markov Chain
P = zeros(6,6); %starting the 6x6 matrix
P(1,1) = (2*l2*t) + (2*l1*t)- (l1+l2)*t;
P(1,2) = (2*l1*t) - (l1 + 2*l2)*t;
P(2,1) = (2*l2*t) + (3*l1*t) - (2*l1+l2)*t;
P(2,2) = (3*l1*t) - (2*l1 + 2*l2)*t;
P(3,1) = (2*l2*t) - (3*l1 + l2) *t;
P(3,2) = (-3*l1 - 2*l2)*t;
P = P - diag(sum(P,2));
%F1 = (l1*t +l1*t);
%F2 = (l2*t + l2*t + l2*t;
%Steady State Probabilty Vector
% Compute EigenVector
% Find the EV corresponding to the '1' EV
% Normalize the eigenvector such that the sum of its elements equals 1,
% The transpose of the result is the steady state probability row-vector.
[X, Y] = eig(P');
fprintf('P = \n');
fprintf('%26.17g %26.17g %26.17g %26.17g %26.17g %26.17g\n', P.');
fprintf('X = \n');
fprintf('%26.17g %26.17g %26.17g %26.17g %26.17g %26.17g\n', X.');
fprintf('Y = \n');
fprintf('%26.17g %26.17g %26.17g %26.17g %26.17g %26.17g\n', Y.');
steadyState = X(:,1)/sum(Y(:,1))
%State Probabilities
n = 100; %Stepping through 100 times
time = 0:t:(n*t); %time vector defined by input variable t as the increment
states = zeros(length(time), 6); % Initializing the state vector with all zeros for the length of time
states(1,:) = steadyState'; %filled with steady state prob vector
PT = expm(P*t);
nnz(~isfinite(states))
nnz(~isfinite(PT))
for i = 2:length(time)
states(i,:) = states (i-1,:) * PT; %state vector is filled by computing states * exponential matrix of P*t
end
nnz(~isfinite(time))
nnz(~isfinite(states))
%Plotting the probabilities over time
plot(time, states)
xlabel("time(s)")
ylabel("State Probability")
title("State Probability Time")
figure;
bar(steadyState);
xlabel('State');
ylabel('Steady-State Probability');
title('Steady-State Probability Vector');
end
Walter Roberson
Walter Roberson 2023-10-23
@John D'Errico
If you have a bit of time, could you look at how the eig(P.') call is generating a diagonal matrix (Y here) that has an all-zero column? When I look around, I see various sources saying that by convention an eigenvector will never be all 0 ??
Mind you, even if the all-zero eigenvector is not possible, it seems plausible to me that the sum() of an given eigenvector could be exactly 0, so I wonder if
steadyState = X(:,1)/sum(Y(:,1))
would be better as
steadyState = X(:,1)/norm(Y(:,1))
?

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Markov Chain Models 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by