{"id":2945,"date":"2012-09-26T14:49:04","date_gmt":"2012-09-26T13:49:04","guid":{"rendered":"http:\/\/www.walkingrandomly.com\/?p=2945"},"modified":"2012-10-08T05:04:19","modified_gmt":"2012-10-08T04:04:19","slug":"randomness-in-matlab-are-you-ruining-your-results-due-to-your-choice-of-syntax","status":"publish","type":"post","link":"https:\/\/walkingrandomly.com\/?p=2945","title":{"rendered":"Randomness in MATLAB : Are you ruining your results due to your choice of syntax?"},"content":{"rendered":"

Pop quiz: What does the following line of MATLAB code do?<\/p>\n

rand('state',10)<\/pre>\n
If you said ‘It changes the seed of the random number generator to 10’<\/em> you get half a point.<\/p>\n
‘Only half a point!?’<\/em> I hear you say accusingly ‘but it says so in my book [for example, 1-3], why not a full point?’<\/em><\/p>\n
You only get a full point if you’d said something like ‘It changes the seed of the random number generator to 10 and it also changes the random number generator from the high quality, default Mersenne Twister generator to a lower quality legacy random number generator.<\/strong>‘<\/em><\/p>\n
OK, how about this one?<\/p>\n
rand('seed',10)<\/pre>\n
This behaves in a very similar manner– it changes both the seed and the type of the underlying generator. However, the random number generator it switches to this time is an even older one that was introduced as far back as MATLAB version 4.\u00a0 It is not very good at all by modern standards!<\/p>\n
A closer look<\/strong><\/p>\n
Open up a fresh copy of a recent version of MATLAB and ask it about the random number generator it’s using<\/p>\n
>> RandStream.getGlobalStream\r\nans =\r\nmt19937ar random stream (current global stream)\r\n             Seed: 0\r\n  NormalTransform: Ziggurat<\/pre>\n
mt1993ar refers to a particular variant of the Mersenne Twister algorithm<\/a>— an industry strength random number generator that’s used in many software packages and simulations.\u00a0 It’s been the default generator in MATLAB since 2007a.\u00a0 Change the seed using the modern (since 2011a), recommended syntax and ask again:<\/p>\n
>> rng(10)\r\n>> RandStream.getGlobalStream\r\nans =\r\nmt19937ar random stream (current global stream)\r\n             Seed: 10\r\n  NormalTransform: Ziggurat<\/pre>\n
This is behaving exactly as you’d expect, you ask it to change the seed and it changes the seed…nothing more, nothing less. Now, let’s use the older syntax<\/p>\n
>> rand('state',10)\r\n>> RandStream.getGlobalStream\r\nans =\r\nlegacy random stream (current global stream)\r\n  RAND algorithm: V5 (Subtract-with-Borrow), RANDN algorithm: V5 (Ziggurat)<\/pre>\n
The random number generator has completely changed!\u00a0\u00a0 We are no longer using the Mersenne Twister algorithm, we are now using a ‘subtract with borrow’ [see reference 4 for implementation details] generator which has been shown to have several undesirable issues [5-7].<\/p>\n
Let’s do it again but this time using the even older ‘seed’ version:<\/p>\n
>> rand('seed',10)\r\n>> RandStream.getGlobalStream\r\nans =\r\nlegacy random stream (current global stream)\r\n  RAND algorithm: V4 (Congruential), RANDN algorithm: V5 (Ziggurat)<\/pre>\n
Now, this random number generator is ancient by computing standards.\u00a0 It also has a relatively tiny period of only 2 billion or so.\u00a0 For details see [4]<\/p>\n
Why this matters<\/strong><\/p>\n
Now, all of this is well documented<\/a> so you may wonder why I am making such a big issue out of it.\u00a0 Here are my reasons<\/p>\n