Matlab doesn’t know how to draw one ball out of an urn containing one ball.
You are stimulating a discrete Markov process. You have a left-stochastic(*) matrix X where 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 j.
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.
Drew Frank said,
July 23, 2010 at 4:50 pm
“The fun part comes in when randsample arbitrarily changes its behavior in response to its first argument having fewer than two elements.”
I actually came here to submit this very wart as a potential blog topic! Earlier this week I had my code break unexpectedly because I wasn’t aware of the changing behavior. I wasn’t using the weighted version so randsample(x,1) didn’t throw an error, but it returned a value in the range 1:x rather than in the singleton “population” [x], which led to an error later in my code. Argh!
Jordi Gutiérrez Hermoso said,
September 28, 2011 at 11:24 pm
Btw, I love the first line here about stimulating a Markov process.
Love the blog, btw. I keep parading it around. Please write more.
crowding said,
February 7, 2012 at 12:11 pm
Would you keep parading it around if I said my opinion of Octave?
MATLAB and Octave are tarred with the same brush, with Octave being rather worse because it’s an intentional replication of all Matlab’s stupidities. I understand how MATLAB got this way, it originated from a community that didn’t communicate much with language researchers and a time before there were good dynamic languages to set the bar. Learn from it and do something better…
I should amend: Octave replicates MATLAB except when it comes to closures, then they break form, declare they don’t like Matlab’s scoping rules for nested functions… and provide no alternatives.
I mean, whatever, if you don’t like setq-type scoping, then implement explicitly-nonlocal scoping or something, but there are a TON of things that are a giant pain to do without at least _some_ kind of lexical closures.
Jordi Gutiérrez Hermoso said,
February 7, 2012 at 12:27 pm
Yes, I know you think Octave is just as bad as Matlab. I think we are only as bads as we have to be. And I try to be as nonbad as possible. I’ve also thought of something like this, but it’s tricky:
http://octave.1599824.n4.nabble.com/A-more-comprehensive-braindead-mode-td4264268.html
I *do* want real closures in Octave too. I just can’t think of a good syntax for them. Nested functions are kind of dumb. I’ve talked to jwe about real closures. The problems are (1) I want real closures and (2) whatever syntax we implement for them, shouldn’t overwrite the existing nested function syntax.
We haven’t implemented nested functions mostly because it’s difficult. Not because we don’t want to.
So what’s the proper syntax for closures that doesn’t break code that already runs in Matlab? The guiding principle should be (even if we’re not there yet): if it works in Matlab, it *should* work in Octave. If it works in Octave, it may break in Matlab. Embrace, extend, extinguish, same dirty tactics they’ve played against us in the past, but this time for the forces of good. >:)
Raster plots and MATLAB « Labrigger said,
December 5, 2011 at 8:45 pm
[...] from an anonymous neuroscientist? Try Abandon MATLAB. This blog’s posts include these gems: Matlab doesn’t know how to draw one ball out of an urn containing one ball The Mathworks don’t even know how to look up functions in their own global namespace And this [...]
Kaushik Ghose said,
March 2, 2012 at 7:00 am
j = nextStates(randsample(length(nextStates), 1, 1, full(weights)));
But, yes, I’m annoyed by Matlab too.
crowding said,
March 2, 2012 at 12:14 pm
Or you could eschew randSample altogether and write (in fewer keystrokes, even)
j = nextStates(find(rand(1) < cumsum(full(weights)), 1, 'last'))
but it's kind of missing the point. randSample is supposed (*) to draw samples from an array you pass it, and perhaps the core axiom of MATLAB's screwed up type system is that everything is an array, even an array of length 1.
Every time you end up building a workaround because the tool you picked out of the toolbox does not behave sensibly, the toolbox loses.
(*) The original documentation reflects the intention of the authors: "Y = RANDSAMPLE(POPULATION,K) returns K values sampled uniformly at
random, without replacement, from the values in the vector POPULATION."
I notice Mathworks reads my blog and has updated the documentation in recent releases to reflect randsample's broken behavior. Guys, next time try fixing the broken behavior instead.
You’re fixing the wrong thing. « Abandon MATLAB said,
March 7, 2012 at 3:03 am
[...] them like normal arrays; you updated the docs on “randsample” to reflect that it draws random samples from arrays only if they are at least 2 elements long; you even updated the docs for “getframe” to clarify that you need to turn off the [...]