{"id":159,"date":"2008-11-05T15:04:17","date_gmt":"2008-11-05T20:04:17","guid":{"rendered":"https:\/\/blogs.mathworks.com\/loren\/2008\/11\/05\/new-ways-with-random-numbers-part-i\/"},"modified":"2016-08-02T14:40:53","modified_gmt":"2016-08-02T19:40:53","slug":"new-ways-with-random-numbers-part-i","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/loren\/2008\/11\/05\/new-ways-with-random-numbers-part-i\/","title":{"rendered":"New Ways With Random Numbers, Part I"},"content":{"rendered":"<div xmlns:mwsh=\"https:\/\/www.mathworks.com\/namespace\/mcode\/v1\/syntaxhighlight.dtd\" class=\"content\">\r\n   <introduction>\r\n      <p>Today I'd like to introduce a guest blogger, Peter Perkins, who is a statistical software developer here at The MathWorks.  In his life prior to The MathWorks, he worked as a statistician\r\n         doing marine mammal research for the U.S. National Marine Fisheries Service, and one of his interests in MATLAB is random\r\n         number generation and Monte-Carlo simulation.\r\n      <\/p>\r\n   <\/introduction>\r\n   <h3>Contents<\/h3>\r\n   <div>\r\n      <ul>\r\n         <li><a href=\"#3\">Random Streams in MATLAB<\/a><\/li>\r\n         <li><a href=\"#6\">Resetting a Random Stream<\/a><\/li>\r\n         <li><a href=\"#9\">Seeding a Stream<\/a><\/li>\r\n         <li><a href=\"#15\">Using Substreams for Reproducibility<\/a><\/li>\r\n         <li><a href=\"#21\">Saving and Restoring State<\/a><\/li>\r\n         <li><a href=\"#32\">Where Not To Do Any Of This<\/a><\/li>\r\n         <li><a href=\"#33\">Next Time<\/a><\/li>\r\n         <li><a href=\"#34\">Other Situations?<\/a><\/li>\r\n      <\/ul>\r\n   <\/div>\r\n   <p>Almost everyone who's used MATLAB has used the <a href=\"https:\/\/www.mathworks.com\/help\/releases\/R2008b\/techdoc\/ref\/rand.html\"><tt>rand<\/tt><\/a> or <a href=\"https:\/\/www.mathworks.com\/help\/releases\/R2008b\/techdoc\/ref\/randn.html\"><tt>randn<\/tt><\/a> functions at one time or another.  The way people use them varies widely, though.  Sometimes, you might just want to fill\r\n      in a vector or matrix with values other than 0 or 1, in order to try out some new lines of code that you've written.  At the\r\n      other end of the spectrum, you might be running a Monte-Carlo simulation where the statistical properties of the pseudo-random\r\n      values matter.\r\n   <\/p>\r\n   <p>But regardless of the context, one of the questions I consistently hear is, \"How and when should I initialize or reinitialize\r\n      or reset MATLAB's random number generator?\" The answer I usually give is, annoyingly, another question: \"Why do you want to\r\n      do that?\"  That's because messing with the state of a random number generator can affect the statistical properties of the\r\n      values it generates, and the right answer to the original question depends on the context - quite often, the answer is \"Don't.\"\r\n       There's a lot that could be said about the <i>when<\/i> and the <i>why<\/i>, and I'll touch on that. But mostly, I want to show some of the <i>how<\/i> that's new in MATLAB R2008b.\r\n   <\/p>\r\n   <p>If you want to run the code that I used to create this blog entry, you'll need R2008b or later, and you'll want to start with\r\n      a fresh MATLAB session, or else execute these two commands (don't worry, I'll explain what they do).\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">    mtstream = RandStream(<span style=\"color: #A020F0\">'mt19937ar'<\/span>);\r\n    RandStream.setDefaultStream(mtstream);<\/pre><h3>Random Streams in MATLAB<a name=\"3\"><\/a><\/h3>\r\n   <p>New in MATLAB R2008b is the addition of random number streams.  A random number stream is just what it sounds like:  a source\r\n      for a sequence of random values.  Usually, people talk about streams of values from a standard uniform distribution, and that's\r\n      the sense that MATLAB uses too. Actually, MATLAB has always had random number streams in the abstract, underneath <tt>rand<\/tt> and <tt>randn<\/tt>.  But in R2008b, streams are something you can interact with directly, and you can even create your own streams. There's\r\n      lots to talk about there, but this article is about how to initialize or reset a stream, so what follows is just a brief summary.\r\n       There's <a href=\"https:\/\/www.mathworks.com\/help\/releases\/R2008b\/techdoc\/math\/brnuahp.html\">lots more<\/a> in MATLAB's documentation.\r\n   <\/p>\r\n   <p>Quite often, you don't need to worry about streams at all; MATLAB creates one for you when it starts up.  <tt>rand<\/tt>, <tt>randn<\/tt>, and the new <a href=\"https:\/\/www.mathworks.com\/help\/releases\/R2008b\/techdoc\/ref\/randi.html\"><tt>randi<\/tt><\/a> function draw values from that stream, and you can generate random numbers without ever knowing or caring what stream they\r\n      come from.  And there's nothing wrong with that - if the random number generators in MATLAB are doing their job correctly\r\n      (and in R2008b, you have a choice among three generators based on state-of-the-art algorithms), you should be able to just\r\n      generate random numbers using <tt>rand<\/tt>, <tt>randn<\/tt> or <tt>randi<\/tt>, and treat everything they return as independent random values.\r\n   <\/p>\r\n   <p>But how do you go back and reproduce the random values you generate or otherwise control how they're generated?  There's a\r\n      new way to do that: you interact with a new kind of object in MATLAB, one that represents a random stream.  Let's get the\r\n      stream that MATLAB creates at startup.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">defaultStream = RandStream.getDefaultStream()<\/pre><pre style=\"font-style:oblique\">defaultStream = \r\nmt19937ar random stream (current default)\r\n             Seed: 0\r\n         RandnAlg: Ziggurat\r\n<\/pre><p>The \"mt19937ar\" indicates that this stream generates values using the most common of the <a href=\"http:\/\/www.math.sci.hiroshima-u.ac.jp\/~m-mat\/MT\/emt.html\">Mersenne Twister<\/a> algorithms, the same algorithm that's been underneath <tt>rand<\/tt> for several releases.  There are other generator choices; more about that in a moment. MATLAB created the stream using a\r\n      seed of 0; more about <i>that<\/i> in a moment as well.\r\n   <\/p>\r\n   <p>What's <tt>RandnAlg<\/tt>? One of the things that's new in R2008b is that <tt>randn<\/tt> draws values from the same stream as <tt>rand<\/tt>, although of course the uniform values have to be transformed to a standard normal distribution. The stream shown here makes\r\n      that transformation using the <a href=\"http:\/\/www.jstatsoft.org\/v05\/i08\">Ziggurat<\/a> algorithm, although you can change it (sorry, no more about that; you'll have to read <a href=\"https:\/\/www.mathworks.com\/help\/releases\/R2008b\/techdoc\/math\/brnuahp.html\">the doc<\/a>).\r\n   <\/p>\r\n   <p>Finally, the display for this stream shows that it is the <i>default stream<\/i>. That just means that it's the stream that <tt>rand<\/tt>, <tt>randn<\/tt>, and <tt>randi<\/tt> draw their values from.  Below, I'll show how you can create your own stream, and make <i>it<\/i> the default.  There are a couple of reasons why you might want to do that.  But first, let's see how to use the default stream\r\n      to reproduce results.\r\n   <\/p>\r\n   <h3>Resetting a Random Stream<a name=\"6\"><\/a><\/h3>\r\n   <p>Sometimes it's useful to be able to reproduce the output from MATLAB's random number generators.  For example, you might want\r\n      to rerun a Monte-Carlo simulation repeatedly, using exactly the same random values because you are debugging the code. In\r\n      fact, MATLAB initializes the default random stream to the same internal state each time it starts up, so the same calls to\r\n      <tt>rand<\/tt>, <tt>randn<\/tt>, and <tt>randi<\/tt> will return the same \"random\" values every you start MATLAB (unless you take steps - more about that in a minute).  But it\r\n      would be tedious to have to restart MATLAB each time you wanted to rerun a simulation.\r\n   <\/p>\r\n   <p>The simplest way to reproduce random number generator output is to reset the default stream, which puts that stream back to\r\n      the state it had when initially created.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">reset(defaultStream)\r\nz1 = randn(1,5)<\/pre><pre style=\"font-style:oblique\">z1 =\r\n      0.53767       1.8339      -2.2588      0.86217      0.31877\r\n<\/pre><p>If I do it again, I'll get the same values again.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">reset(defaultStream)\r\nz2 = randn(1,5)<\/pre><pre style=\"font-style:oblique\">z2 =\r\n      0.53767       1.8339      -2.2588      0.86217      0.31877\r\n<\/pre><p>While simple, <a href=\"https:\/\/www.mathworks.com\/help\/releases\/R2008b\/techdoc\/ref\/randstream.reset.html\"><tt>reset<\/tt><\/a> is kind of a blunt instrument, and provides reproducibility only on a very gross scale - it's as if you restarted MATLAB.\r\n       It's not something you should do unless you really want to reuse the same random sequence over again.\r\n   <\/p>\r\n   <h3>Seeding a Stream<a name=\"9\"><\/a><\/h3>\r\n   <p>All pseudo-random number generators have the notion of a seed.  Think of a stream as a very long sequence of values, wrapped\r\n      into a circle, because it repeats when exhausted.  Think of a seed as a starting position on the circle that you specify when\r\n      you create a stream.  There are (usually) many more random values than there are seeds, so it's not possible to start just\r\n      anywhere by specifying a seed.  And there's usually no simple description of where, say, the seed 0 lies in relation to, say,\r\n      the seeds 1 or 2 other than knowing that they represent different starting positions. But if you create a stream with the\r\n      same seed, you'll get the same sequence of values every time.\r\n   <\/p>\r\n   <p>When MATLAB starts, it creates the default stream with commands equivalent to<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">stream0 = RandStream(<span style=\"color: #A020F0\">'mt19937ar'<\/span>,<span style=\"color: #A020F0\">'Seed'<\/span>,0)\r\nRandStream.setDefaultStream(stream0);<\/pre><pre style=\"font-style:oblique\">stream0 = \r\nmt19937ar random stream\r\n             Seed: 0\r\n         RandnAlg: Ziggurat\r\n<\/pre><p>We saw that earlier.  If I create a new stream with a different seed, such as<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">stream1 = RandStream(<span style=\"color: #A020F0\">'mt19937ar'<\/span>,<span style=\"color: #A020F0\">'Seed'<\/span>,1)<\/pre><pre style=\"font-style:oblique\">stream1 = \r\nmt19937ar random stream\r\n             Seed: 1\r\n         RandnAlg: Ziggurat\r\n<\/pre><p>and make <i>it<\/i> the default stream,\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">RandStream.setDefaultStream(stream1);<\/pre><p>I'll get a <i>different<\/i> sequence of random values.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">randn(1,5)<\/pre><pre style=\"font-style:oblique\">ans =\r\n     -0.64901       1.1812     -0.75845      -1.1096     -0.84555\r\n<\/pre><p>But if I always use the same seed, I'll always get the same random number stream.  This is a lot like resetting the current\r\n      default stream - in fact, <tt>reset<\/tt> is nothing more than reseeding.  The main difference between how I used <tt>reset<\/tt> above and here is that here, I've created my own stream, so I know <i>exactly<\/i> what I'm getting.  Before, I'd get repeatability, but the specific random values depended on what generator algorithm was\r\n      underneath the current default stream.\r\n   <\/p>\r\n   <p>Reseeding is kind of a big hammer - you would not want to create new streams every time you needed a set of random values.\r\n      It's most useful if you use it as an initialization step, perhaps at MATLAB startup, perhaps before running a simulation.\r\n   <\/p>\r\n   <h3>Using Substreams for Reproducibility<a name=\"15\"><\/a><\/h3>\r\n   <p>We've seen that a drawback of using <tt>reset<\/tt> or <tt>'Seed'<\/tt> to reproduce results is that you essentially have to start over from the beginning.  A new feature in R2008b, substreams,\r\n      provides a new way around that drawback.\r\n   <\/p>\r\n   <p>Remember the circle?  Substreams are nothing more than fixed \"checkpoints\" in a random number stream, defined by their position\r\n      around the circle, spaced evenly at very large intervals.  The substreams themselves are indexed beginning at 1, and you can\r\n      jump to the beginning of a substream by setting a stream's <tt>Substream<\/tt> property.  Values generated from different substreams are statistically independent in the same sense as consecutive values\r\n      in a stream are.\r\n   <\/p>\r\n   <p>Two new generator algorithms in R2008b support substreams: the Combined Multiple Recursive (<tt>'mrg32k3a'<\/tt>) and the Multiplicative Lagged Fibonacci (<tt>'mlfg6331_64'<\/tt>) generators (the doc has <a href=\"https:\/\/www.mathworks.com\/help\/releases\/R2008b\/techdoc\/math\/brrztpq.html\">references<\/a>). I'll create one of the former, and make it the default stream.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">stream = RandStream(<span style=\"color: #A020F0\">'mrg32k3a'<\/span>);\r\nRandStream.setDefaultStream(stream);<\/pre><p>To reposition a stream to the beginning of a particular substream, I just set its <tt>Substream<\/tt> property.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">stream.Substream = 2;<\/pre><p>Much like seeds, substreams are nothing more than aliases for particular points along a sequence of random values.  So what's\r\n      different about them? There are some differences that I'll get into next time, but for now, the most obvious difference is\r\n      that you can jump between substreams without having to create a new stream.  It's a fairly lightweight operation.\r\n   <\/p>\r\n   <p>Let's see what happens if I position a random number stream to the beginning of a substream before each iteration of a loop.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\"><span style=\"color: #0000FF\">for<\/span> i = 1:5\r\n    stream.Substream = i;\r\n    z = rand(1,i)\r\n<span style=\"color: #0000FF\">end<\/span><\/pre><pre style=\"font-style:oblique\">z =\r\n      0.72701\r\nz =\r\n       0.5582      0.85273\r\nz =\r\n      0.16662      0.29236      0.77278\r\nz =\r\n      0.34773      0.38864      0.80161      0.95025\r\nz =\r\n      0.56709      0.68377      0.29879      0.59318      0.69247\r\n<\/pre><p>Now I can reproduce results from the 3rd iteration simply by repositioning the stream to the corresponding substream.  I didn't\r\n      have to rerun the whole loop, and I didn't even have to know beforehand that I wanted to look at this iteration.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">i = 3;\r\nstream.Substream = i;\r\nz = rand(1,i)<\/pre><pre style=\"font-style:oblique\">z =\r\n      0.16662      0.29236      0.77278\r\n<\/pre><p>If that loop were a time-consuming Monte-Carlo simulation, it would be a real boon to be able to go recreate and investigate\r\n      any specific iteration without rerunning the whole thing.\r\n   <\/p>\r\n   <p>Keep in mind that you don't need to overuse substreams.  They're there to use, you don't have to worry about \"using up\" all\r\n      the values in each substream before moving to the next one, but it would be pointless to take the other extreme and jump to\r\n      a different substream every time you generate a new value. Substreams don't add randomness, they just make it easier to reproduce\r\n      values.\r\n   <\/p>\r\n   <h3>Saving and Restoring State<a name=\"21\"><\/a><\/h3>\r\n   <p><tt>rand<\/tt> and <tt>randn<\/tt> are <a href=\"http:\/\/en.wikipedia.org\/wiki\/Pseudorandom_number_generator\">pseudorandom number generators<\/a>, and all pseudorandom generator algorithms have an internal state that determines what the next value will be.  It's a deterministic\r\n      algorithm. Remember the circle?  Each point along the circle corresponds to a different state.  But for almost all modern\r\n      algorithms, the state is more complicated than just the last number generated or a simple index of the current position in\r\n      the stream.\r\n   <\/p>\r\n   <p>Let's get the state of the default generator.  To do that, I'll read the <tt>State<\/tt> property.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">defaultStream = RandStream.getDefaultStream()\r\nsavedState = defaultStream.State;\r\nwhos <span style=\"color: #A020F0\">savedState<\/span><\/pre><pre style=\"font-style:oblique\">defaultStream = \r\nmrg32k3a random stream (current default)\r\n             Seed: 0\r\n         RandnAlg: Ziggurat\r\n  Name             Size            Bytes  Class     Attributes\r\n\r\n  savedState      12x1                48  uint32              \r\n\r\n<\/pre><p>If I call <tt>rand<\/tt> to generate a few random values,\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">u1 = rand(1,5)<\/pre><pre style=\"font-style:oblique\">u1 =\r\n      0.83911      0.51069      0.84469      0.68938      0.23777\r\n<\/pre><p>and then I restore the state that I saved,<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">defaultStream.State = savedState;<\/pre><p><tt>rand<\/tt> will generate those exact same values the next time I call it.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">u2 = rand(1,5)<\/pre><pre style=\"font-style:oblique\">u2 =\r\n      0.83911      0.51069      0.84469      0.68938      0.23777\r\n<\/pre><p>The same is also true for <tt>randn<\/tt> and for <tt>randi<\/tt>.  And because all three functions draw values from the same default stream, there's only one state to worry about - the default\r\n      stream's state.  I can use the state I saved earlier to reproduce results from <tt>randn<\/tt>, for example.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">defaultStream.State = savedState;\r\nz1 = randn(1,5)<\/pre><pre style=\"font-style:oblique\">z1 =\r\n     0.048564      0.87031      0.52365     0.095609     -0.79233\r\n<\/pre><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">defaultStream.State = savedState;\r\nz2 = randn(1,5)<\/pre><pre style=\"font-style:oblique\">z2 =\r\n     0.048564      0.87031      0.52365     0.095609     -0.79233\r\n<\/pre><p>It's important to remember that because there's a single default stream, calling one function affects what specific values\r\n      the others return.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">savedState = defaultStream.State;\r\nu1 = rand(1,5)\r\nu2 = rand(1,5)<\/pre><pre style=\"font-style:oblique\">u1 =\r\n     0.045212      0.48758      0.75138      0.36639       0.5676\r\nu2 =\r\n      0.10969      0.27046      0.10359      0.87516      0.28295\r\n<\/pre><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">defaultStream.State = savedState;\r\n<span style=\"color: #228B22\">% This returns the same vector as before<\/span>\r\nu1 = rand(1,5)\r\n<span style=\"color: #228B22\">% Call randn in between the two calls to rand<\/span>\r\nz = randn(1,5);\r\n<span style=\"color: #228B22\">% This returns different values than before<\/span>\r\nu2 = rand(1,5)<\/pre><pre style=\"font-style:oblique\">u1 =\r\n     0.045212      0.48758      0.75138      0.36639       0.5676\r\nu2 =\r\n      0.20425     0.043083      0.54273      0.10648      0.20139\r\n<\/pre><p>The functions do not act independently with respect to the values they return, however, the values generated by the three\r\n      functions can still be treated as <i>statistically<\/i> independent regardless of what order you generate them in.\r\n   <\/p>\r\n   <p>One practical application of saving and restoring the state might be to reproduce an entire Monte-Carlo simulation run, by\r\n      saving the state before running it. Another might be to recreate one particularly interesting iteration in a Monte-Carlo simulation,\r\n      having saved the state before that iteration. Another might be to recreate just one call to <tt>rand<\/tt> that happened deep in your code, so that you can recreate that call elsewhere, perhaps for debugging.  As long as you have\r\n      saved the state at the appropriate point, you can \"jump to\" <i>anywhere<\/i> in the sequence of random values. That's the advantage over resetting or seeding a stream, but it's the catch too: you have\r\n      to have saved the right state beforehand.\r\n   <\/p>\r\n   <p>As an aside, the internal states of the various generators available in MATLAB have different sizes and formats, and you can\r\n      see what they are by looking at the <tt>State<\/tt> property of different streams.  But I recommend not getting too familiar with, or modifying, the inner workings of these\r\n      state vectors.  Read 'n' write 'em, yes, perhaps even rope 'n' brand 'em, but don't try and understand 'em.\r\n   <\/p>\r\n   <h3>Where Not To Do Any Of This<a name=\"32\"><\/a><\/h3>\r\n   <p>Sometimes, discretion is the better part of valor.  Be careful to understand why and when you do any of this.  Remember, resetting\r\n      or replacing the default stream, or writing to its state, or jumping to a substream affects all subsequent calls to <tt>rand<\/tt>, <tt>randn<\/tt>, and <tt>randi<\/tt>, and is probably not something you want to do deep in your code - too easy to forget that it is in there.  Replacing the\r\n      default stream, especially, is probably not something you need to do a lot of.\r\n   <\/p>\r\n   <h3>Next Time<a name=\"33\"><\/a><\/h3>\r\n   <p>We've seen four ways to reproduce results from MATLAB's random number generators.<\/p>\r\n   <div>\r\n      <ul>\r\n         <li>Resetting a stream<\/li>\r\n         <li>Creating a new stream with a fixed seed<\/li>\r\n         <li>Substreams<\/li>\r\n         <li>Explicitly saving and restoring the generator state<\/li>\r\n      <\/ul>\r\n   <\/div>\r\n   <p>But sometimes what's needed is the \"opposite of reproducibility\" - you want to make sure that different runs of a simulation\r\n      use <i>different<\/i> random inputs, because you want to be able to assume that the different runs are statistically independent.  Until <a href=\"https:\/\/blogs.mathworks.com\/loren\/?p=161\">next time ...<\/a><\/p>\r\n   <h3>Other Situations?<a name=\"34\"><\/a><\/h3>\r\n   <p>Can you think of any cases where you've needed to reproduce random numbers that can't be handled by the examples I've shown?\r\n       Let me know <a href=\"https:\/\/blogs.mathworks.com\/loren\/?p=159#respond\">here<\/a>.\r\n   <\/p><script language=\"JavaScript\">\r\n<!--\r\n\r\n    function grabCode_609fa39fce4d41e0a1068235dcf65f36() {\r\n        \/\/ Remember the title so we can use it in the new page\r\n        title = document.title;\r\n\r\n        \/\/ Break up these strings so that their presence\r\n        \/\/ in the Javascript doesn't mess up the search for\r\n        \/\/ the MATLAB code.\r\n        t1='609fa39fce4d41e0a1068235dcf65f36 ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' 609fa39fce4d41e0a1068235dcf65f36';\r\n    \r\n        b=document.getElementsByTagName('body')[0];\r\n        i1=b.innerHTML.indexOf(t1)+t1.length;\r\n        i2=b.innerHTML.indexOf(t2);\r\n \r\n        code_string = b.innerHTML.substring(i1, i2);\r\n        code_string = code_string.replace(\/REPLACE_WITH_DASH_DASH\/g,'--');\r\n\r\n        \/\/ Use \/x3C\/g instead of the less-than character to avoid errors \r\n        \/\/ in the XML parser.\r\n        \/\/ Use '\\x26#60;' instead of '<' so that the XML parser\r\n        \/\/ doesn't go ahead and substitute the less-than character. \r\n        code_string = code_string.replace(\/\\x3C\/g, '\\x26#60;');\r\n\r\n        author = 'Loren Shure';\r\n        copyright = 'Copyright 2008 The MathWorks, Inc.';\r\n\r\n        w = window.open();\r\n        d = w.document;\r\n        d.write('<pre>\\n');\r\n        d.write(code_string);\r\n\r\n        \/\/ Add author and copyright lines at the bottom if specified.\r\n        if ((author.length > 0) || (copyright.length > 0)) {\r\n            d.writeln('');\r\n            d.writeln('%%');\r\n            if (author.length > 0) {\r\n                d.writeln('% _' + author + '_');\r\n            }\r\n            if (copyright.length > 0) {\r\n                d.writeln('% _' + copyright + '_');\r\n            }\r\n        }\r\n\r\n        d.write('<\/pre>\\n');\r\n      \r\n      d.title = title + ' (MATLAB code)';\r\n      d.close();\r\n      }   \r\n      \r\n-->\r\n<\/script><p style=\"text-align: right; font-size: xx-small; font-weight:lighter;   font-style: italic; color: gray\"><br><a href=\"javascript:grabCode_609fa39fce4d41e0a1068235dcf65f36()\"><span style=\"font-size: x-small;        font-style: italic;\">Get \r\n            the MATLAB code \r\n            <noscript>(requires JavaScript)<\/noscript><\/span><\/a><br><br>\r\n      Published with MATLAB&reg; 7.7<br><\/p>\r\n<\/div>\r\n<!--\r\n609fa39fce4d41e0a1068235dcf65f36 ##### SOURCE BEGIN #####\r\n%% New Ways With Random Numbers, Part I\r\n% Today I'd like to introduce a guest blogger,\r\n% <https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/loadAuthor.do?objectId=1093801&objectType=author Peter Perkins>,\r\n% who is a statistical software developer here at The MathWorks.  In his life\r\n% prior to The MathWorks, he worked as a statistician doing marine mammal\r\n% research for the U.S. National Marine Fisheries Service, and one of his\r\n% interests in MATLAB is random number generation and Monte-Carlo simulation.\r\n\r\n%%\r\n% Almost everyone who's used MATLAB has used the\r\n% <https:\/\/www.mathworks.com\/help\/releases\/R2008b\/techdoc\/ref\/rand.html |rand|>\r\n% or\r\n% <https:\/\/www.mathworks.com\/help\/releases\/R2008b\/techdoc\/ref\/randn.html |randn|>\r\n% functions at one time or another.  The way people use them varies\r\n% widely, though.  Sometimes, you might just want to fill in a vector or\r\n% matrix with values other than 0 or 1, in order to try out some new lines of\r\n% code that you've written.  At the other end of the spectrum, you might be\r\n% running a Monte-Carlo simulation where the statistical properties of the\r\n% pseudo-random values matter.\r\n%\r\n% But regardless of the context, one of the questions I consistently hear is,\r\n% \"How and when should I initialize or reinitialize or reset MATLAB's random\r\n% number generator?\" The answer I usually give is, annoyingly, another\r\n% question: \"Why do you want to do that?\"  That's because messing with the\r\n% state of a random number generator can affect the statistical properties of\r\n% the values it generates, and the right answer to the original question\r\n% depends on the context - quite often, the answer is \"Don't.\"  There's a lot\r\n% that could be said about the _when_ and the _why_, and I'll touch on that. But\r\n% mostly, I want to show some of the _how_ that's new in MATLAB R2008b.\r\n\r\n%%\r\n% If you want to run the code that I used to create this blog entry, you'll\r\n% need R2008b or later, and you'll want to start with a fresh MATLAB session,\r\n% or else execute these two commands (don't worry, I'll explain what they\r\n% do).\r\n    mtstream = RandStream('mt19937ar');\r\n    RandStream.setDefaultStream(mtstream);\r\n\r\n\r\n%% Random Streams in MATLAB\r\n% New in MATLAB R2008b is the addition of random number streams.  A random\r\n% number stream is just what it sounds like:  a source for a sequence of\r\n% random values.  Usually, people talk about streams of values from a standard\r\n% uniform distribution, and that's the sense that MATLAB uses too. Actually,\r\n% MATLAB has always had random number streams in the abstract, underneath\r\n% |rand| and |randn|.  But in R2008b, streams are something you can interact\r\n% with directly, and you can even create your own streams. There's lots to\r\n% talk about there, but this article is about how to initialize or reset\r\n% a stream, so what follows is just a brief summary.  There's\r\n% <https:\/\/www.mathworks.com\/help\/releases\/R2008b\/techdoc\/math\/brnuahp.html lots more>\r\n% in MATLAB's documentation.\r\n%\r\n% Quite often, you don't need to worry about streams at all; MATLAB creates\r\n% one for you when it starts up.  |rand|, |randn|, and the new\r\n% <https:\/\/www.mathworks.com\/help\/releases\/R2008b\/techdoc\/ref\/randi.html |randi|>\r\n% function draw values from that stream, and you can generate random numbers\r\n% without ever knowing or caring what stream they come from.  And there's\r\n% nothing wrong with that - if the random number generators in MATLAB are\r\n% doing their job correctly (and in R2008b, you have a choice among three generators\r\n% based on state-of-the-art algorithms), you should be able to just generate\r\n% random numbers using |rand|, |randn| or |randi|, and treat everything they\r\n% return as independent random values.\r\n%\r\n% But how do you go back and reproduce the random values you generate or\r\n% otherwise control how they're generated?  There's a new way to do that:\r\n% you interact with a new kind of object in MATLAB, one that represents a\r\n% random stream.  Let's get the stream that MATLAB creates at startup.\r\ndefaultStream = RandStream.getDefaultStream()\r\n%%\r\n% The \"mt19937ar\" indicates that this stream generates values using the most\r\n% common of the\r\n% <http:\/\/www.math.sci.hiroshima-u.ac.jp\/~m-mat\/MT\/emt.html Mersenne Twister>\r\n% algorithms, the same algorithm that's been underneath |rand| for several\r\n% releases.  There are other generator choices; more about that in a moment.\r\n% MATLAB created the stream using a seed of 0; more about _that_ in a moment\r\n% as well.\r\n%%\r\n% What's |RandnAlg|? One of the things that's new in R2008b is that |randn|\r\n% draws values from the same stream as |rand|, although of course the uniform\r\n% values have to be transformed to a standard normal distribution. The stream\r\n% shown here makes that transformation using the\r\n% <http:\/\/www.jstatsoft.org\/v05\/i08 Ziggurat>\r\n% algorithm, although you can change it (sorry, no more about that; you'll\r\n% have to read\r\n% <https:\/\/www.mathworks.com\/help\/releases\/R2008b\/techdoc\/math\/brnuahp.html the doc>).\r\n%\r\n% Finally, the display for this stream shows that it is the _default stream_.\r\n% That just means that it's the stream that |rand|, |randn|, and |randi| draw\r\n% their values from.  Below, I'll show how you can create your own stream, and\r\n% make _it_ the default.  There are a couple of reasons why you might want to do\r\n% that.  But first, let's see how to use the default stream to reproduce results.\r\n\r\n\r\n%% Resetting a Random Stream\r\n% Sometimes it's useful to be able to reproduce the output from MATLAB's\r\n% random number generators.  For example, you might want to rerun a\r\n% Monte-Carlo simulation repeatedly, using exactly the same random values\r\n% because you are debugging the code. In fact, MATLAB initializes the default\r\n% random stream to the same internal state each time it starts up, so the same\r\n% calls to |rand|, |randn|, and |randi| will return the same \"random\" values\r\n% every you start MATLAB (unless you take steps - more about that in a\r\n% minute).  But it would be tedious to have to restart MATLAB each time you\r\n% wanted to rerun a simulation.\r\n%\r\n% The simplest way to reproduce random number generator output is to reset\r\n% the default stream, which puts that stream back to the state it had when\r\n% initially created.\r\nreset(defaultStream)\r\nz1 = randn(1,5)\r\n%%\r\n% If I do it again, I'll get the same values again.\r\nreset(defaultStream)\r\nz2 = randn(1,5)\r\n%%\r\n% While simple,\r\n% <https:\/\/www.mathworks.com\/help\/releases\/R2008b\/techdoc\/ref\/randstream.reset.html |reset|>\r\n% is kind of a blunt instrument, and provides\r\n% reproducibility only on a very gross scale - it's as if you restarted\r\n% MATLAB.  It's not something you should do unless you really want to reuse\r\n% the same random sequence over again.\r\n\r\n\r\n%% Seeding a Stream\r\n% All pseudo-random number generators have the notion of a seed.  Think of a\r\n% stream as a very long sequence of values, wrapped into a circle, because it\r\n% repeats when exhausted.  Think of a seed as a starting position on the\r\n% circle that you specify when you create a stream.  There are (usually) many more\r\n% random values than there are seeds, so it's not possible to start just\r\n% anywhere by specifying a seed.  And there's usually no simple description of\r\n% where, say, the seed 0 lies in relation to, say, the seeds 1 or 2 other than\r\n% knowing that they represent different starting positions. But if you create\r\n% a stream with the same seed, you'll get the same sequence of values every\r\n% time.\r\n\r\n%%\r\n% When MATLAB starts, it creates the default stream with commands\r\n% equivalent to\r\nstream0 = RandStream('mt19937ar','Seed',0)\r\nRandStream.setDefaultStream(stream0);\r\n%%\r\n% We saw that earlier.  If I create a new stream with a different seed, such\r\n% as\r\nstream1 = RandStream('mt19937ar','Seed',1)\r\n%%\r\n% and make _it_ the default stream,\r\nRandStream.setDefaultStream(stream1);\r\n%%\r\n% I'll get a _different_ sequence of random values.\r\nrandn(1,5)\r\n%%\r\n% But if I always use the same seed, I'll always get the same random number\r\n% stream.  This is a lot like resetting the current default stream - in fact,\r\n% |reset| is nothing more than reseeding.  The main difference between how I used\r\n% |reset| above and here is that here, I've created my own stream, so I know _exactly_\r\n% what I'm getting.  Before, I'd get repeatability, but the specific random\r\n% values depended on what generator algorithm was underneath the current\r\n% default stream.\r\n%\r\n% Reseeding is kind of a big hammer - you would not want to create new\r\n% streams every time you needed a set of random values. It's most useful if\r\n% you use it as an initialization step, perhaps at MATLAB startup, perhaps\r\n% before running a simulation.\r\n\r\n\r\n%% Using Substreams for Reproducibility\r\n% We've seen that a drawback of using |reset| or |'Seed'| to reproduce results\r\n% is that you essentially have to start over from the beginning.  A new\r\n% feature in R2008b, substreams, provides a new way around that drawback.\r\n%\r\n% Remember the circle?  Substreams are nothing more than fixed \"checkpoints\"\r\n% in a random number stream, defined by their position around the circle,\r\n% spaced evenly at very large intervals.  The substreams themselves are\r\n% indexed beginning at 1, and you can jump to the beginning of a substream by\r\n% setting a stream's |Substream| property.  Values generated from different\r\n% substreams are statistically independent in the same sense as consecutive\r\n% values in a stream are.\r\n%\r\n% Two new generator algorithms in R2008b support substreams: the Combined\r\n% Multiple Recursive (|'mrg32k3a'|) and the Multiplicative Lagged\r\n% Fibonacci (|'mlfg6331_64'|) generators (the doc has\r\n% <https:\/\/www.mathworks.com\/help\/releases\/R2008b\/techdoc\/math\/brrztpq.html references>).\r\n% I'll create one of the former, and make it the default stream.\r\nstream = RandStream('mrg32k3a');\r\nRandStream.setDefaultStream(stream);\r\n\r\n%%\r\n% To reposition a stream to the beginning of a particular substream, I just\r\n% set its |Substream| property.\r\nstream.Substream = 2;\r\n\r\n%%\r\n% Much like seeds, substreams are nothing more than aliases for particular\r\n% points along a sequence of random values.  So what's different about them?\r\n% There are some differences that I'll get into next time, but for now, the\r\n% most obvious difference is that you can jump between substreams without\r\n% having to create a new stream.  It's a fairly lightweight operation.\r\n%\r\n% Let's see what happens if I position a random number stream to the beginning\r\n% of a substream before each iteration of a loop.\r\nfor i = 1:5\r\n    stream.Substream = i;\r\n    z = rand(1,i)\r\nend\r\n%%\r\n% Now I can reproduce results from the 3rd iteration simply by repositioning\r\n% the stream to the corresponding substream.  I didn't have to rerun the whole\r\n% loop, and I didn't even have to know beforehand that I wanted to look at\r\n% this iteration.\r\ni = 3;\r\nstream.Substream = i;\r\nz = rand(1,i)\r\n%%\r\n% If that loop were a time-consuming Monte-Carlo simulation, it would be a\r\n% real boon to be able to go recreate and investigate any specific iteration\r\n% without rerunning the whole thing.\r\n\r\n%%\r\n% Keep in mind that you don't need to overuse substreams.  They're there to\r\n% use, you don't have to worry about \"using up\" all the values in each\r\n% substream before moving to the next one, but it would be pointless to take\r\n% the other extreme and jump to a different substream every time you generate\r\n% a new value. Substreams don't add randomness, they just make it easier to\r\n% reproduce values.\r\n\r\n\r\n%% Saving and Restoring State\r\n% |rand| and |randn| are\r\n% <http:\/\/en.wikipedia.org\/wiki\/Pseudorandom_number_generator pseudorandom number generators>,\r\n% and all pseudorandom generator algorithms have an internal state that\r\n% determines what the next value will be.  It's a deterministic algorithm.\r\n% Remember the circle?  Each point along the circle corresponds to a different\r\n% state.  But for almost all modern algorithms, the state is more complicated\r\n% than just the last number generated or a simple index of the current\r\n% position in the stream.\r\n%\r\n% Let's get the state of the default generator.  To do that, I'll read the\r\n% |State| property.\r\ndefaultStream = RandStream.getDefaultStream()\r\nsavedState = defaultStream.State;\r\nwhos savedState\r\n%%\r\n% If I call |rand| to generate a few random values,\r\nu1 = rand(1,5)\r\n%%\r\n% and then I restore the state that I saved,\r\ndefaultStream.State = savedState;\r\n%%\r\n% |rand| will generate those exact same values the next time I call it.\r\nu2 = rand(1,5)\r\n%%\r\n% The same is also true for |randn| and for |randi|.  And because all three\r\n% functions draw values from the same default stream, there's only one state\r\n% to worry about - the default stream's state.  I can use the state I saved\r\n% earlier to reproduce results from |randn|, for example.\r\ndefaultStream.State = savedState;\r\nz1 = randn(1,5)\r\n%%\r\ndefaultStream.State = savedState;\r\nz2 = randn(1,5)\r\n%%\r\n% It's important to remember that because there's a single default stream,\r\n% calling one function affects what specific values the others return.\r\nsavedState = defaultStream.State;\r\nu1 = rand(1,5)\r\nu2 = rand(1,5)\r\n%%\r\ndefaultStream.State = savedState;\r\n% This returns the same vector as before\r\nu1 = rand(1,5)  \r\n% Call randn in between the two calls to rand\r\nz = randn(1,5); \r\n% This returns different values than before\r\nu2 = rand(1,5)  \r\n%%\r\n% The functions do not act independently with respect to the values they\r\n% return, however, the values generated by the three functions can still be\r\n% treated as _statistically_ independent regardless of what order you generate\r\n% them in.\r\n\r\n%%\r\n% One practical application of saving and restoring the state might be to\r\n% reproduce an entire Monte-Carlo simulation run, by saving the state before\r\n% running it. Another might be to recreate one particularly interesting\r\n% iteration in a Monte-Carlo simulation, having saved the state before that\r\n% iteration. Another might be to recreate just one call to |rand| that\r\n% happened deep in your code, so that you can recreate that call elsewhere,\r\n% perhaps for debugging.  As long as you have saved the state at the\r\n% appropriate point, you can \"jump to\" _anywhere_ in the sequence of random\r\n% values. That's the advantage over resetting or seeding a stream, but it's\r\n% the catch too: you have to have saved the right state beforehand.\r\n\r\n%%\r\n% As an aside, the internal states of the various generators available in\r\n% MATLAB have different sizes and formats, and you can see what they are by\r\n% looking at the |State| property of different streams.  But I recommend not\r\n% getting too familiar with, or modifying, the inner workings of these state\r\n% vectors.  Read 'n' write 'em, yes, perhaps even rope 'n' brand 'em, but\r\n% don't try and understand 'em.\r\n\r\n\r\n%% Where _Not_ To Do Any Of This\r\n% Sometimes, discretion is the better part of valor.  Be careful to understand\r\n% why and when you do any of this.  Remember, resetting or replacing the\r\n% default stream, or writing to its state, or jumping to a substream affects\r\n% all subsequent calls to |rand|, |randn|, and |randi|, and is probably not\r\n% something you want to do deep in your code - too easy to forget that it is\r\n% in there.  Replacing the default stream, especially, is probably not\r\n% something you need to do a lot of.\r\n\r\n\r\n%% Next Time\r\n% We've seen four ways to reproduce results from MATLAB's random number\r\n% generators. \r\n%\r\n% * Resetting a stream\r\n% * Creating a new stream with a fixed seed\r\n% * Substreams\r\n% * Explicitly saving and restoring the generator state\r\n%\r\n% But sometimes what's needed is the \"opposite of\r\n% reproducibility\" - you want to make sure that different runs of a simulation\r\n% use _different_ random inputs, because you want to be able to assume that\r\n% the different runs are statistically independent.  Until \r\n% <https:\/\/blogs.mathworks.com\/loren\/?p=161 next time ...>\r\n\r\n\r\n%% Other Situations?\r\n% Can you think of any cases where you've needed to reproduce random numbers\r\n% that can't be handled by the examples I've shown?  Let me know \r\n% <https:\/\/blogs.mathworks.com\/loren\/?p=159#respond here>.\r\n\r\n\r\n##### SOURCE END ##### 609fa39fce4d41e0a1068235dcf65f36\r\n-->","protected":false},"excerpt":{"rendered":"<p>\r\n   \r\n      Today I'd like to introduce a guest blogger, Peter Perkins, who is a statistical software developer here at The MathWorks.  In his life prior to The MathWorks, he worked as a... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/loren\/2008\/11\/05\/new-ways-with-random-numbers-part-i\/\">read more >><\/a><\/p>","protected":false},"author":39,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[16,6,28],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/159"}],"collection":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/users\/39"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/comments?post=159"}],"version-history":[{"count":1,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/159\/revisions"}],"predecessor-version":[{"id":1893,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/159\/revisions\/1893"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/media?parent=159"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/categories?post=159"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/tags?post=159"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}