- make sure you use hold on
- Include a marker in your plot() options, such as plot(x, y, '-*')
Program Is not Plotting Correctly
1 visualización (últimos 30 días)
Mostrar comentarios más antiguos
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
0 comentarios
Respuestas (1)
Walter Roberson
el 22 de Oct. de 2023
If you are using plot() inside a loop but nothing is showing up on the plot, then:
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 comentarios
Walter Roberson
el 23 de Oct. de 2023
Editada: Walter Roberson
el 23 de Oct. de 2023
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))
markov_chain(sym(2e-4), sym(1e-5), sym(0.1))
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
el 23 de Oct. de 2023
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))
?
Ver también
Categorías
Más información sobre Markov Chain Models en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!