1
$\begingroup$

I am currently trying to recreate results in a paper where I solve for the periodic steady state of the following rate matrix:

enter image description here

Where epsilon is some force parameter, I have found the left and right eigenvectors and the eigenvalues of the matrix to diagonalize it, using these results to construct the actual transition matrix. The issue is that this transition matrix is quite complicated, even using //FullSimplify.

The equation that I am using to solve for the transition matrix is: enter image description here

Where $Pd$ is a $3x6$ identity matrix (it is an I, sorry for the bad notation), $U$ and $V$ are the right and left eigenvectors of the rate matrix, respectively, $N$ is the product of $U$ and $V$, $Λ$ is the diagonalized matrix (which is being propagated through time over time steps $τ$), and lastly, $M$ is a $6x3$ identity matrix coupled with the probability of the outgoing bit value of the system of interest.

My question is, once I get $T$, can I apply Eigensystem[T] to it and receive the periodic steady state for the preferred eigenvalue of $1$? The problem is when I solve for the eigenvalues of $T$, I don't get an eigenvalue of $1$ unless I start plugging numbers in for $P0$ and $P1$, but I want to preserve these probabilities because the whole point is to solve for their fluctuations over time. Does this mean that the transition matrix I solved for is wrong? Or is there a system of functions that deal with Markov processes in Mathematica that I am missing? Sorry for the longwinded question; I am a Mathematica beginner.

Improve this question
$\endgroup$
5
  • 1
    $\begingroup$ It seems that there is a lot of variation across disciplines in the terminology people use to describe stochastic processes. Could you add some references to your question to help more people understand it? e.g. what do you mean by "periodic steady state"? $\endgroup$
    – Chris K
    Commented Feb 20 at 18:45
  • $\begingroup$ Welcome to Mathematica.SE! I suggest the following: 1) Take the tour and check the faqs. 2) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! 3) As you receive help, try to give it too, by answering questions in your area of expertise. $\endgroup$
    – Chris K
    Commented Feb 20 at 18:45
  • $\begingroup$ Hi Chris, by periodic steady state, I mean solving for the eigenvector of the transition matrix T with eigenvalue 1. $\endgroup$
    – magg13__
    Commented Feb 22 at 16:07
  • $\begingroup$ I'm not sure if this is exactly what you want, but couldn't you get the same steady state by solving for the eigenvector of the rate matrix with eigenvalue 0? $\endgroup$
    – Chris K
    Commented Feb 22 at 20:02
  • $\begingroup$ Thanks for the response. Yes this makes sense, but the problem is that the eigenvalues for the transition matrix T are neither 1 nor 0. Is this because of the symbolic evaluation? In the paper, they were able to preserve the symbols and still solve for the periodic steady state so I'm not sure where I went wrong... $\endgroup$
    – magg13__
    Commented Feb 24 at 1:21

1 Answer 1

Reset to default
0
$\begingroup$

If I understand the problem correctly, you can get the stationary distribution straight from the transition rate matrix, as the eigenvector associated with the dominant, zero eigenvalue (conservation of probability mass). See e.g. section 4.4 of these notes.

a := {
  {-1, 1, 0, 0, 0, 0},
  {1, -2, 1, 0, 0, 0},
  {0, 1, -2 + e, 1 + e, 0, 0},
  {0, 0, 1 - e, -2 - e, 1, 0},
  {0, 0, 0, 1, -2, 1},
  {0, 0, 0, 0, 1, -1}
}

N[Eigenvalues[a]]
(* {-3.73205, -3., -2., -1., -0.267949, 0.} *)

Looks like the last eigenvector is the one we want.

Eigenvectors[a][[-1]]
(* {-((1 + e)/(-1 + e)), -((1 + e)/(-1 + e)), -((1 + e)/(-1 + e)), 1, 1, 1} *)

A numerical example:

e = 0.1;
(* normalize eigenvector to sum to one *)
p = Eigenvectors[a][[-1]]/Total[Eigenvectors[a][[-1]]]
ListPlot[p, PlotRange -> {0, All}]
(* {0.183333, 0.183333, 0.183333, 0.15, 0.15, 0.15} *)

enter image description here

You could probably also use the built-in ContinuousMarkovProcess but I find the syntax inconvenient.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.