You are stimulating a discrete Markov process. You have a left-stochastic(*) matrix
X(j,i) gives the probability of transitioning from state
i to the state
j. There are a lot of states, and most of the state transition probabilities are zero, so to fit your system into memory, you have built
X as a sparse matrix.
You have the index of the current state,
i; you want to simulate the next step and randomly draw a value for the index the next state
Easy enough to begin with, we find the indices and probabilities of the potential next states like this:
nextStates = find(X(:,i));
weights = X(nextStates,i);
So now we need to generate a random draw from a discrete distribution, where we know the probability of each value of the distribution. Isn’t there a function that does that? Oh yes,
randsample. After looking through the help for randsample (a task that is more difficult than it sounds) we might write:
j = randsample(nextStates, 1, 1, full(weights));
Now a really easy question: How and why does this break? Answer below the fold. (hint: look at the category of the post)
(*) while almost all the mathematical literature uses right-stochastic matrices, in MATLAB the sparse matrices are mush slower if you use them that way.
Indeed, if your dynamics are such that some states have just one transition out of them, then
nextStates will have one element, and
weights will have just the value 1. The fun part comes in when
randsample arbitrarily changes its behavior in response to its first argument having fewer than two elements.
In effect, if you ask MATLAB to randomly draw single ball out of an urn containing one red ball, it comes back with an “error:”
"there are fewer than 'red' balls in this urn."
Actually, what it says is:
W must have length equal to N.
Which is no more helpful because you did supply a W that was equal in length to N.