Modeling a discrete probabilistic time evolution in matlab

Question

Good day, I am trying to model the following situation.

There is a system, which can be in three states. State M1, M2 and M3. From state M1, it can go into state 2 with probability p1. From state M2 it can go into state 1 with probability p2, and into state 3 with probability p3. From state M3 it can go into state 2 with probability p4.

Now my first attempt at tackling this was by working with a system of three coupled differential equations

equation

But this doesn't really seem to work, as when I try to plot the solutions I'm not really seeing the expected dynamics.

So instead I want to try and solve it in a discrete sense, but I honestly don't know where to begin. My guess is that one of the easier ways of doing so is by using matlab?

As for boundary conditions, I'd like to consider starting in M1.

Now I understand that I am really not giving you much, but I honestly don't have any idea where to begin. So my question basically boils down to how should I start tackling this problem, and is there any closely related problem that I can look up, to see how they handle that?

Kind regards

Edit: It might be useful to note that I am doing this in order to find out how long, on average, the system spends in each state. I want to eventually fit this to data, in order to find the probabilities.

Sounds like you have a [Markov chain](http://en.wikipedia.org/wiki/Markov_chain) - you could start by finding the state transition matrix. — bdecaf, Jan 14 '14 at 14:28
That is excellent! Thank you. I have the matrix. However, searching for Markov Chain Simulation Matlab looks.. daunting. — user129412, Jan 14 '14 at 14:57
when you have the data and want to find the probabilities, it is actually a easier problem. as bdecaf suggested, it is a Markov chain and you just need to calculate the conditionally probability (e.g. probablility of state one following state two) empirically. — Cici, Jan 14 '14 at 14:59
you are nearly done. Simulation is just repeatedly multiplying the matrix by itself (times initial state). As a tip - calculate the eigenvalues and eigenvectors. — bdecaf, Jan 14 '14 at 14:59
Hm, shouldn't I do something with drawing a random variable T times (if T is the total time in steps of 1) using rand(T,1), and then comparing it to the probabilities, or something of the sorts? — user129412, Jan 14 '14 at 15:10
So far we spoke about probabilites. What you suggest is simulating a *realization* of the process. Of course you can do that, it's just less efficient as the calculation by probabilities. — bdecaf, Jan 14 '14 at 15:15
Hm alright, I see what you mean. I suppose I Should have a bit more clear. What I am looking to do is actually simulating the state of the system, versus time. — user129412, Jan 14 '14 at 15:19

score 0 · Accepted Answer · edited May 23 '17 at 12:18

As said in the comments, this looks a lot like a Markov chain. In your case, the transition matrix M is given by

M = [1-p1, p2,      0
     p1,   1-p2-p3, p4
     0,    p3,      1-p4];

Note that all columns have to sum to 1. The current probability is given by the state vector, which is a column vector that also has to sum to 1. Suppose you start in state 1: state = [1; 0; 0], then the state vector after n transitions is M^n * state. See e.g. here for the math.

But, as you said, you want to know how to do a concrete simulation of such a process, not how the probabilities change over time. This can be done like so (untested):

state = 1; % start state

while not_stopped()
    do_stuff(state);
    transition_prob = M(:, state); % get the appropriate column
    state = weighted_randi(transition_prob)
end

Replace do_stuff by whatever you want to do for that state, for example recording the current state in some vector for later plotting. The function weighted_randi can be defined like so (inspired by this answer):

function R = weighted_randi(weights)
% returns an integer from the range 1:length(weights) with probability
% given by the weights

[~, R] = histc(rand(), cumsum([0; weights(:) ./ sum(weights)]));

score 0 · Answer 2 · edited May 23 '17 at 10:33

I tried working it out using Bas's method outlined above, and although I did get somewhere, I wasn't able to completely pull it off. However, with all of the comments I have received I was able to make my search more accurate, and I found this post

Using this, I was able to do exactly what I wanted with the code given by the main poster. I understand that it is much, much better to write your own code instead of copying someone else's, but I have made sure to understand every single step in there and for my purposes I am satisfied with this. I want to thank you again for all the help!

Modeling a discrete probabilistic time evolution in matlab

2 Answers2