Program Is not Plotting Correctly

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

Answer 1

[Walter Roberson]

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.

[Walter Roberson]

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]

@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))

?