**You are now following this question**

- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.

# Simulating a Continuous time markov chain

131 views (last 30 days)

Show older comments

##### 0 Comments

### Accepted Answer

Dana
on 27 Aug 2020

See here:

##### 20 Comments

Nazer Hdaifeh
on 27 Aug 2020

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

Dana
on 27 Aug 2020

Dana
on 28 Aug 2020

Nazer Hdaifeh
on 30 Aug 2020

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
on 31 Aug 2020

neeti gupta
on 15 Jan 2021

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
on 15 Jan 2021

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
on 15 Jan 2021

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

neeti gupta
on 16 Jan 2021

I am working with MATLAB R2015a

neeti gupta
on 16 Jan 2021

Dana
on 18 Jan 2021

"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
on 18 Jan 2021

"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
on 19 Jan 2021

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
on 19 Jan 2021

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
on 19 Jan 2021

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
on 2 Feb 2021

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
on 2 Feb 2021

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
on 2 Feb 2021

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
on 11 Sep 2021

Hi Susman,

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

Thanks

### More Answers (1)

Dana
on 3 Feb 2021

Edited: Dana
on 3 Feb 2021

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.

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!**An Error Occurred**

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list:

## How to Get Best Site Performance

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

### Americas

- América Latina (Español)
- Canada (English)
- United States (English)

### Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)

### Asia Pacific

- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)