Probability vector from a Markov Transition Matrix

Can somebody tell me how to calculate the probability matrix from a Markov Transition Matrix which values are known? I'm not familiar with this calulations in Matlab.

This is the Markov Transition Matrix:

Note: All transition rates of the matrix are known values except for the variable "θpm"

Which script lines should I enter into Matlab to find the values of the P vector for a given "θpm"? And how can I plot a graphic "AbUn" vs. "θpm"? knowing that AbUn=P4+P8 and "P4" and P8" are variables from the probability vector (P).

Thanks a lot in advance

[William Rose]

@Daniel Caballero, You said you want to find the probability matrix. Do you mean the probability vector?

I wll assume that you mean the probability vector. It is good to recognize from the problem statement that you are being asked to find the equilibrium probability distribution, i.e. the distribution which, once attained, will stay the same. It is also helpful to recognize that this is an eigenvector problem: you are finding the eigenvector with the largest eigenvalue, for the state transition matrix (or Markov transition matrix) T. But is an issue: matrix T-I is singular. But there's a way out: since the probbailities sum to 1, eight probabiliities determine the 9th one uniquely. So we can reduce the rank of the matrix by one, and it won't be singular. This will take some algebra, because can't just delete one row and one column. You have to adjust the rates that remain. You do this by replacing p2 (if you choose p2) with 1-(p1+p3+p4+p5+p6+p7+p8+p9), and then do some stuff which I forget, but which you have probably learned. (I took a course on this in 1983.) I recommend that you delete p1, p2, p6, or p7, because the other probabilities are needed to compute AbUn and RelUn. I would try deleting p2, because its row and column have the fewest off-diagonal elements, so the algebra associated with deleting it may be simpler than for the other options.

As for the graph, once you do the work above, you will have a theoretical 8x8 matrix. Then you assume a value for θpm, and plug in all the other known values. Then you solve, ie. you find the eigenvector with the biggest eigenvalue. Then add up the probabilities p4+p8=AbUn from that eigenvector, and p3+p5+p9=RelUn, and save those values, associated with the θpm which you used. Then repeat with different values of θpm. When you have a sufficient set of θpm values and associated RelUn and AbUn, plot θpm vs RelUn and θpm vs. AbUn.

[Daniel Caballero]

@William Rose, first of all, thank you for your comment. You are right, I was talking about the probability vector (I have edited my question). I am wondering if there are some direct functions in Matlab to calculate the eigenvector as you mentioned and save all the algebra needed. In case there aren't (because of the singularity), do you know where to find concrete information about what will come after replacing p2 with 1-(p1+p3+p4+p5+p6+p7+p8+p9)? And how to proceed with the calculations in Matlab?

[William Rose]

@Daniel Caballero,

I actually didn;t know that T-I would be non-diagonalizable. I would plug in b]avlues for the constants (including thetaPM) and then try

[V,D]=eig(T);

That will return the matrix of eigenvectors, V, and the vector of eigenvalues, D. As I recall, Matlab orders them so the first element of D is the largest, and the first column of V is the corresponding eigenvector, the one you want. They may be normalized to unit length. But that is Euclidean length. You want the sum of the elements (not the sum of the squares of the elements) to be one, to be a probability vector, so multiply the eigenvector by an appropriate constant to make it so. If any elements of the eigenvector are negative then something went wrong.

Once you have T, an easy way to estimate the eigenvector, or check your work, is to apply T repeatdly to any initial vector, such as v0=[1 0 0 0 0 0 0 0 0 0]. v1=T*v0, v2=T*v1, v3=T*v2, and so on. After 10-20 tries, v should converge and stop changing. This is the eigenvector! It should match the one found with eig().

Did you notice that the equations for the diagonal elements, such as a22=1-Sn, are such that the sum of all the elements in each row will be one? That is a feaure of Markov transition matrices.

[Daniel Caballero]

@William Rose,

I tried to do what you told me and I got stucked as some of the elements of the eigenvector resulted in negative numbers. I can see what could be the problem. The elements of the transition matrix (T) aren't all positive numbers because I'm using the steady-state frequencies for the elements of the transition matrix T (which are in units: "operations per hour"). I think I need to convert those frequencies into probabilities so adding the elements of the same row of matrix T results into 1 (how can I convert them?. Please, correct me if I'm wrong or if there's something else.

I'm sending you my Matlab script and I'm attaching the paper that contains the problem I want to solve. I hope it's very clear for you. Thanks again

[William Rose]

@Daniel Caballero,

What you have requested is a considerable effort, at least for me. I could not find a textbook section for this topic that presented the material in the way that is most helpful to me and to you, so I wrote my own - see attached PDF. I think it will help you. I am also attaching a program that converts a matrix of rate constants, k, to the probability transition matrix, A, when the time step, tau, is specified. Then it finds the equilbrium probability distribution, pEq. You can adapt this to your situation.

Please provide the complete manuscript citation, so I can see the rest of it, including references. Thanks.

Specific notes for this model

Be sure to convert all rates to have the same time units in their denominator.

• What are the values of Fc, Fpp, Fst? I don't see them.

The rates in this model vary by over ten orders of magnitude: Sn=43200/hr, Fcc=1e-6/hr.

I think this makes the A matrix (or maybe A-I) ill conditioned, i.e. the matrix inverse can be very sensitive to small round-off errors. This may negatively affect the robustness of the computed equilibrium probability (pEq). The very large difference in time constants means it is a stiff system, which means it will be challenging to do numerical integration.

Since Sn (rate from state 2 to 6) and Sb (rate from state 4 to 8) are much faster than any other time constants in the model, any probability that gets to state 2 or 4 will immediately be transferred into state 6 or 8. Therefore I predict the equilbrium probability in 2 and 4 will be extremely low. Therefore I would rewrite the model without states 2 or 4, for computational simplicity. The rates that used to go into states 2 and 4 will now go directly into states 6 and 8. This will remove the two very fast rates. There remain the very slow rate Fc, Fst, Fpp, and Fcc=1e-6/hr (which equas about 1.1%/year,since 1 yr=8760 hr). Leave them in the model as is, for now.

[William Rose]

I updated my answer above at 9:40 pm EDT by replacing the PDF with a new version. The new PDF corrects some mistakes.

[Daniel Caballero]

@William Rose, I really appreciate your work and your feedback. I have to tell you that the script you sent me helped a lot since I got the result I was expecting. I didn't rewrite the model without the states 2 or 4 as you suggested me and I didn't have any problem as far as I know. I'm sending you the script I adapted, there you will see that the plot for the variables AbUn vs Ipm replicates the "Figure 2" of the paper problem. I'm also sending you the script that replicates de graphic of "Figure 3" and the complete paper as required.

The values for Fpp and Fst can be calculated from the equations shown in the first page of the attached PDF (paper problem.pdf). These variables depend on the ST factor which is 0 for the first excersice (Figure 2) as we assume a relay that does not have self-testing:

So, Fst=Fp*ST=0 and Fpp=Fp(1-ST)=Fp=1/MTBF=1/50 failures/hour which is 1/(50*8760) failure/hour. Finally, Fc is the component failures per year. As shown in the above image we have that Fc=2 failures/year which is 2/8760 failures/hour. These equations are also written in the Matlab script I'm sending you.

For these excercises, I see that it's necessary to use a tau=1 in order to replicate the same figures shown in the paper. Regarding the PDF explanation, can you please explain me what would be the the reason to have different values of (tau) for this specific exercise (paper) since it doesn't mention anything about it?

[William Rose]

@Daniel Caballero, you're welcome. Thank you for posting the full paper. The figure I posted, with ST=0, 50%, 100%, matches Fig.6 of the paper (which I had not seen when I made and posted the figure), except I did not include ST=90%. As for tau: I am not convinced that one must use tau=1 to match the paper results. My results appear to be a good match, and I used tau=1e-2. In an early version of the program, I saw small differences in the equilibrium probability when I changed tau, but as tau became small, the equilibrium probability converged to a stable level, as expected.

[Daniel Caballero]

Hello @William Rose, I still have one doubt if you could help me, regarding the transition rate matrix (k), why do we have to assume a value of zero (0) for all the values of kii of the nine states?. I know that with this values, we can get the same results as the paper. However, why does the same paper tell us that kii=1-ki (different from 0) where ki is the total rate of leaving the state i?

About it, do you agree to have a contact meeting for clarifying some related things and to present you another project I'm intereseted to develop with your possible help?

Kind regards,

Daniel Caballero

[William Rose]

Daniel,

I reviewed our previous discussion, above. It appears that I did not post the file I intended to post, in one of my past comments. Therefore I am posting now a document about finding the transition probability matrix. Please see if it answers your question. Please email me at [e-mail] to dicsuss possible follow-up.

[William Rose]

Daniel,

In the PDF document which I posted last night, I specifically did not define , because I did not need it to derive the matrix A. is the derived the rate per unit time of going from state i to j. Therefore is the "rate of going from state i to i", i.e. rate of doing nothing. Is it possible to define or measure the rate of doing nothing? I don't know. But I don't need to define A, so I am OK with not knowing.

Answer 1

[William Rose]

@Daniel Caballero, I used the code and ideas I posted in the earlier answer, with appropriate constants for your particular model, to try to reproduce Fig.2 of the paper you posted. It worked. Remember that Fig. 2 was Ab.Un. verus test interval, for a relay with no self test capability, i.e. ST=0. The red trace in the attached figure shows the result from the code, and it is an excellent or perfect match to Fig.2. of the paper.

I also calculated the results for a relay with 50% self test effectiveness and 100% self test effectiveness. Those results are also plotted. The results show that if a relay has 100% self test efficiency, then one should not inspect it at all, because inspecting it takes it out of service for about 2 hours, and does not cause more rapid fault detection.