Simulating a Continuous time markov chain

Question

Nazer Hdaifeh 2020-8-27

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/585344-simulating-a-continuous-time-markov-chain

评论： Basma Bargal 2021-9-11

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?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Dana 2020-8-27

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/585344-simulating-a-continuous-time-markov-chain#answer_485915

See here:

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

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

20 个评论
显示 18更早的评论隐藏 18更早的评论

Dana 2021-1-15

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.

Dana 2021-1-19

在 MATLAB Online 中打开

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.

Dana 2021-2-2

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

在 MATLAB Online 中打开

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;
-4.40770000000000	0.829800000000000	1.44820000000000	2.12970000000000;
0.672400000000000	-2.56400000000000	0	1.89160000000000;
0.413900000000000	0	-2.53480000000000	2.12090000000000;
0	0	0	0]

Basma Bargal 2021-9-11

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 2021-2-3

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/585344-simulating-a-continuous-time-markov-chain#answer_613757

编辑：Dana 2021-2-3

在 MATLAB Online 中打开

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.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

susman 2021-2-3

Great that makes sense. Thank you so much!

请先登录，再进行评论。

Simulating a Continuous time markov chain

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

20 个评论
显示 18更早的评论隐藏 18更早的评论

更多回答（1 个）

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

Community Treasure Hunt

Simulating a Continuous time markov chain

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

20 个评论 显示 18更早的评论隐藏 18更早的评论

更多回答（1 个）

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

20 个评论
显示 18更早的评论隐藏 18更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论