There seems to be many followup questions, it may be worth discussing the problem in some depth, how you might attack it in MATLAB. So lets start out with a discussion of such a Markov process, and how we would work with it. First, create a simple markov process. I'm not feeling terribly creative right now, so lets just pick something simple, thus a 5x5 transition matrix.
T = triu(rand(5,5),-1);
T = T./sum(T,2)
T =
0.17362 0.0029508 0.33788 0.19802 0.28752
0.059036 0.16812 0.29036 0.36644 0.11604
0 0.15184 0.30054 0.275 0.27262
0 0 0.42829 0.30672 0.26499
0 0 0 0.28417 0.71583
I've arbitrarily chosen a matrix with no absorbing state, but one where at each step, the process will tend to push the state to a higher number, but some chance you can climb backwards too.
If we look at the matrix above, if you are in state 5, with probability 0.71583 you will stay in state 5, but 28% of the time, you will drop back to state 4, etc. Next, consider a vector that describes your current state. Suppose we start out in state 1 (thus initially, 100% of the time, we are in state 1.)
X = [1 0 0 0 0];
After one step of the process, we don't know what state we will be in, but we can determine the probability that we will lie in any given state. That is found simply using a matrix multiplication.
% After one time step
X = X*T
X =
0.17362 0.0029508 0.33788 0.19802 0.28752
% After two time steps
>> X = X*T
X =
0.030319 0.052314 0.24588 0.27082 0.40066
% After three time steps
>> X = X*T
X =
0.0083525 0.04622 0.21532 0.28971 0.44039
% After four time steps
>> X = X*T
X =
0.0041788 0.040491 0.20504 0.29181 0.45848
Gradually, the probability grows that we will lie in state 5 MOST of the time. So after 100 time steps, the probability keeps growing that we lie in state 5.
X100 = [1 0 0 0 0]*T^100
X100 =
0.0025378 0.035524 0.19457 0.29167 0.4757
Does a steady state prediction of the long term state of this process exist? Yes. We can get that from the eigenvector of T', corresponding to the eigenvalue 1. That would be the vector V1, such that V1*T = V1.
eig(T')
ans =
1 + 0i
0.48873 + 0i
0.16513 + 0i
0.0054905 + 0.13467i
0.0054905 - 0.13467i
[V,D] = eig(T');
V1 = V(:,1)';
V1 = V1/sum(V1)
V1 =
0.0025378 0.035524 0.19457 0.29167 0.4757
That last step was to normalize the eigenvector to sum to 1, so they could be viewed as probabilities. Remember that eig normalizes its vectors to have unit norm.
So over the long term, the probability is 47.57% that the process lies in state 5.
This is what we can learn about the long term behavior of that system. But how about simulating the process? So, instead of thinking about where we will be as this process goes to infinity, can we simulate a SINGLE instance of such a Markov chain? This is a very different thing, since it does not rely on eigenvalues, matrix multiplication, etc.
Now, we need to look at the rows of T. I'll build this as a random walk that runs for many time steps. At any time, I'll simulate the progress of the random walk as it proceeds from one state to the next state. Then I'll simulate 1000 such random walks, and see where we ended at the final step of that process.
The trick is to use the cumulative sum of T, along the rows of T. What does that do for us?
CT = cumsum(T,2)
CT =
0.17362 0.17657 0.51445 0.71248 1
0.059036 0.22716 0.51752 0.88396 1
0 0.15184 0.45239 0.72738 1
0 0 0.42829 0.73501 1
0 0 0 0.28417 1
Suppose we start out in state 1. The first row of CT is pertinent here. Generate a single random number (using rand). If that number is less than 0.17362, then we started in state 1, and will remain in state 1. If the number lies between 0.51445 and 0.71248, then we will have moved to state 4, etc. We should see how to simulate this process now. I'll simulate 10000 such random Markov processes, each for 1000 steps. I'll record the final state of each of those parallel simulations. (If I had the parallel tolbox, I suppose I could do this using a parfor loop, but not gonna happen for me.)
I'll use histcounts, but older MATLAB releases need to use histc. I'll append a column of zeros at the beginning of CT to make histcounts return the proper bin index.
CT = [zeros(size(T,1),1),cumsum(T,2)];
numberOfSims = 10000;
simTime = 1000;
initialState = 1;
currentState = repmat(initialState,[1,numberOfSims]);
for n = 1:simTime
r = rand(1,numberOfSims);
for m = 1:numberOfSims
% for each sim, where does r(m) lie, relative to the boundaries
% in the corresponding row of CT?
[~,~,currentState(m)] = histcounts(r(m),CT(currentState(m),:));
end
end
finalState = currentState;
This took a minute or so to run, but that was a fair amount of work. The first 10 such parallel simulations ended in these states:
finalState(1:10)
ans =
2 4 5 3 5 5 5 5 4 5
Now, how often did our process end in any given state? We can count that using accumarray.
finalStateFreq = accumarray(finalState',ones(numberOfSims,1))/numberOfSims
finalStateFreq =
0.0029
0.0336
0.1957
0.2891
0.4787
If I did a good job here, this should mesh well with the steady state predictions from before. Was the simulation time long enough? Surely yes. This is a small Markov process, with only 5 states. And 10K total sims will be entirely adequate to predict if the result matches the steady state predictions.
V1 =
0.0025378 0.035524 0.19457 0.29167 0.4757
That I got a consistent answer from the simulations with the steady state predictions suggests I did a good job in the simulation. As you would expect from a random simulation, it will only approach the theoretical frequency as the number of simulations goes to infinity. With a little extra effort, I suppose could even have estimated the variance of the counts in each bin, based on the variance of a binomial distribution.
