Simulating a Continuous time markov chain

I have a transition matrix Q of 5 states (5x5), reoccurrence is allowed. the final state is state five (Death) and initial state is State one (no disease) and the other states are the levels of the disease(lets say breast cancer). I want to simulate the markov chain using matlb ... any one can help me with that please?

Answer 1

[Dana]

See here:

https://www.mathworks.com/matlabcentral/answers/57961-simulating-a-markov-chain

Also, if you have the econometrics toolbox, you can use simulate.

[Nazer Hdaifeh]

Thank you for sharing, my question was about continuous time markov chain not discrete time MC.

[Dana]

Sorry, I missed that in the title. As there are different ways to characterize the same continuous-time Markov-chain (CTMC) process, can you describe exactly what Q is, and any other information you have? For example, one common way would be to have a vector of 5 exponential distribution parameters that govern the processes for leaving each of the 5 states, and then a transition probability matrix P, the ij-th element of which gives the probability of transitioning from state i to state j, conditional on the fact that the process is leaving state i (which means that the diagonal elements of P must all be 0). Is this what you have?

[Dana]

I updated some old code I had written to simulate continuous-time Markov chains (CTMC) and uploaded it to the File Exchange here. You might find this helpful. It takes, as the key input, the "transition intensity" matrix Q for the CTMC, where the diagonal elements are the negatives of the exponential parameters governing jumps out of each state, and the off-diagonals in a given row govern the relative likelihood of jumping to each of the other states in that row, conditional on a jump happening.

[Nazer Hdaifeh]

Thank you for the reply Dana,

the matrix I have is for a Markov chains where transitions occur in continuous time transition-rate matrix:

The non-diagonal elements of this matrix must be greater than 0. The diagonal elements, , are chosen as follows:

, for i = 1, . . . ,N ( to ensure the row sums are 0)

is the rate with which the CTMC, X(t), moves from state i to state j, given that X(t) is in state i. #

As a question if I may ask:

what is T in the funtion simCTMC(Q,T,nsim,[instt]) you made, as I understood it is the period until terminate the simulation? what if I need it to continue simulating until all patients are in the state 5 (death).

[Dana]

Yes, T is the period that the simulation is terminated. The code was not originally built to do something like simulate until an absorbing state is reached. I've made a minor modification so that you can do that now by setting T=Inf (you'll need to download the updated files from the file exchange to use this feature). FYI, this will only work if there is an absorbing state (obviously), and if an absorbing state will surely be reached at some point, regardless of the initial state (i.e., in every stationary distribution, there must be a state with probability 1 attached to it).

[neeti gupta]

Hello Dana

I am trying to use your code for simulating my CTMC with Q = [-0.08 0.02 0.02; 0.01 -0.03 0.02; 0.01 0.01 -0.02], T=[0, 200], nsim=20, instt=[1 0 0]. but getting the following error

Error using ./

Matrix dimensions must agree.

Error in simCTMC (line 52)

Pi = Pi ./sum(Pi,2);

Kindly help.. I will be highly grateful..

Neeti

[Dana]

Hi Neeti,

I can't reproduce your error because the inputs you gave are not valid:

1) The rows of Q need to sum to 0, but that's not the case for the first row in yours.

2) T should be a scalar not a vector. It's just the length of the simulation. So I'm guessing you want T=200.

3) nsim and instt are not both used. If you want to choose specific initial states, you pass instt as a vector where each element is the index of a state number between 1 and m, where m is the dimension of Q. The number of simulations nsim is then set to the number of elements of instt. If you were to pass in a value of nsim in this case, that value would be ignored.

If you don't specify an instt, then you need to specify nsim as the number of parallel simulations to run. The first state is then chosen randomly from the stationary distribution associated with Q.

So for example, if you pass instt = [1; 3; 2; 2];, then the code will do four parallel simulations (regardless of whether you also pass in some other value of nsim), where the first simulation starts with state 1, the second with state 3, and the last two with state 2.

Looking at your instt, you can see right away that there's a problem: you're telling the program to run 3 simulations, with the second 2 starting with state 0. But there's no state 0, there are only states 1, 2 and 3. So that won't work. In addition, you've also passed in nsim=20, but that's just being ingnored.

[neeti gupta]

Hello Dana,

I really appreciate your gesture.. You have helped me a lot to understand the terms in code. However, even after implementing corrections as per your suggestion I am getting the same error.

simulating Q = [-0.08 0.04 0.04; 0.01 -0.05 0.04; 0.01 0.01 -0.02] for T=200 (each simulation for 200 hrs), instt=[1;1;1;1;1;1;1;1;1;1] and nsim=10 as I want 10 simulation runs each starting with state 1.

error:

Error using ./

Matrix dimensions must agree.

Error in simCTMC (line 52)

Pi = Pi ./sum(Pi,2);

Kindly guide me where I am going wrong

Thanks

Neeti

[Dana]

Hi Neeti, I can't reproduce the error. Running the code with these inputs doesn't cause any problems for me. Which version of MATLAB are you using?

[neeti gupta]

I am working with MATLAB R2015a

[neeti gupta]

Please tell me whether the code works for any number of states or do we need to modify something if number of states changes.

[Dana]

"I am working with MATLAB R2015a"

Ah, I think that explains the error. To fix it, you can replace the offending line of code:

Pi = Pi./sum(Pi,2);

with:

Pi = bsxfun(@rdivide,Pi,sum(Pi,2));

The downside of this fix is that I can't promise the same problem won't pop up elsewhere in the code (though you should be able to come up with a similar fix if it does). Of course, the best solution is to update to a more recent version of MATLAB (R2016b or later in this case), but I understand that may not be possible for you.

[Dana]

"Please tell me whether the code works for any number of states or do we need to modify something if number of states changes."

I'm not sure exactly what you're asking. The number of states is determined based on the dimensions of the square Q matrix. So the number of states is changed if and only if you pass in a different-sized Q. And Q can be a square matrix of any size you want, so in that sense there's no restriction on how many states you can have.

[neeti gupta]

Dear Dana

I cant thank you enough for all the help you are providing. I have implemented the fix as per your suggestion.The previous error has been resolved by it. However, I am getting a new error down in the code:

Error using <=

Matrix dimensions must agree.

Error in simCTMC (line 148)

[~,newstt] = max(p<=cuPijmp,[],2);

Please see if you can help.

Thanks in advance

Neeti

[Dana]

Hi Neeti,

Yes, this is what I suspected might happen. It's the same fundamental problem as the other one, which has to do with implicit expansion when doing element-wise operations.

As an example, suppose A is 5x5 and B is 5x1, and I want to add the vector B to each column of A (so that the i-th column of the result is A(:,i)+B). Mathematically speaking, you can't do A+B since + here indicates element-wise addition, and that only works if A and B are the same size (which they're not here). As a result, in MATLAB R2016a and earlier, if you tried to do A+B you'd get an error about the matrix dimensions not agreeing.

You could of course work around this by, say, creating a 5x5 matrix C, and then looping through the columns to set C(:,i)=A(:,i)+B. Or you could replicate the vector B into a 5x5 matrix using repmat(B,1,5) and then add the result directly to A. But these are not very efficient ways to do this computationally. A better approach is to use the function bsxfun. In this example, you can do bsxfun(@plus,A,B) which tells MATALB to automatically "implicitly expand" B to the same size as A in a way that's mathematically equivalent to the repmat method, but in a much more computationally efficient way. Note that bsxfun works for many element-wise operations, not just + (you can see the list on the help page).

Starting in In MATLAB R2016b, the behavior provided by bsxfun now happens automatically when you use (many) element-wise operations. So in these more recent versions of MATLAB, A+B automatically does exactly what bsxfun(@plus,A,B) used to do. The code I wrote takes advantage of this more simple approach, but as a result is not backward-compatible with R2016a and earlier.

To bring this back around to your case, you can replace the offending line of code

[~,newstt] = max(p<=cuPijmp,[],2);

with

[~,newstt] = max(bsxfun(@le,p,cuPijmp),[],2);  

So here, the variables p and cuPijmp are not the same size: they have the same number of rows, but cuPijmp has more columns. In R2016b and later, this isn't a problem when applying the element-wise operation <=, since MATLAB automatically "replicates" p so that it has the same number of columns as cuPijmp, but in earlier versions of MATLAB (such as yours) it throws an error. You instead need to use the bsxfun approach.

[neeti gupta]

Hey Dana Thank you so much for your help. I am finally able to run the code and simulate my CTMC. Further with the help of your wonderful explaination I have developed a great understanding of code.

Thanks a lot..

[susman]

Dear Dana,

I have seen your code and I would say that it works perfectly well. It helped me to understand few aspects of CTMC, which were very confusing.

I have a question relevant to your code.

Question: "Does this code work in the following setting?"

I have a non-homogenous transition matrix that varies over time. e.g.

Q(t) = [q11 q12 q13; q21 q22 q33; q31 q32 q33]

and each element of Q(t) varies with time like q11 = 0.67 + 0.012*age. That means Q(t) will be different for each age.

Second, which algorithm did you use in your coding? any literature reference?

[Dana]

Hi susman,

> Question: "Does this code work in the following setting?"

No, the code isn't set up for anything like that unfortunately. Q has to be constant. I'm not sure exactly how you'd go about modifying it for a time-varying Q, as the time until the next switch is no longer a simple exponential random variable. There's probably a way to deal with this, but I don't know off-hand how.

> Second, which algorithm did you use in your coding? any literature reference?

Fundamentally, the code works as follows. For a given simulation, conditional on the system having jumped to state i at date t, first draw a value τ from the exponential distribution with parameter (the absolute value of the ith diagonal element of the Q matrix), which gives a random time until the next jump, i.e., the next time the system switches states will be at date . To determine which state the system switches into at , pick randomly from among the other states (i.e., from every state except i) using the relative intensities given in the non-diagonal elements of the ith row of Q to control the probabilities of switching into each state.

My reference is a textbook by J.R. Norris called Markov Chains, which is a nice undergraduate-level text that covers both discrete- and continuous-time Markov chains.

[susman]

Thank you Dana. I have to work with non-homogeneous continuous time markov chains. But as a first looking a homogeneous version gives a very clear picture through your code. Please let me know if you know any good reference for non-homogeneous CTMC especially the coding part.

Secondly, regarding your code, I have following Q matrix with one absorpbing state as death and I set the parameters as T = 65 and nsims = 1000. But my simulations remain a straight line and there are no transitions. All simullations remain in state zero from the very begining. However, if I remove the absorbing state, it works well. Can you explain the reason for that?

Q = [-2.79760000000000	0.628400000000000	0	0	2.16920000000000;
0	-4.40770000000000	0.829800000000000	1.44820000000000	2.12970000000000;
0	0.672400000000000	-2.56400000000000	0	1.89160000000000;
0	0.413900000000000	0	-2.53480000000000	2.12090000000000;
0	0	0	0	0]

[Basma Bargal]

Hi Susman,

Wondering whether you found any code or references on how to go about non-homogenous MArkov Chains in Matlab?

Thanks

Answer 2

[Dana]

I don't know offhand of any non-homogeneous CTMC references, sorry.

Regarding your case, this part of the help section regarding ths inputs of simCTMC.m is relevant:

%   nsim: number of simulations to run (only used if instt is not passed in)
%   instt: optional vector of initial states; if passed in, nsim = size of
%       instt; otherwise, nsim draws are made from the stationary
%       distribution of the Markov chain (if there are multiple stationary
%       distributions, an error is returned).

So if you don't pass in a vector of initial states in the variable instt, the program randomly draws nsim initial states from the stationary distribution of the MC. But when you have an absorbing state, the stationary distribution of the MC has probability 1 associated with the absorbing state, and 0 for every other state. If you draw the initial state randomly from this stationary distribution, you're always going to start with the absorbing state, and then of course you'll stay in that state forever, which obviously isn't very interesting.

Long story short, for cases like these you probably want to pick your own initial vector instt. How to pick it is up to you.

[susman]

Great that makes sense. Thank you so much!