Re: Mersenne Twister range in rand.oct and randn.oct

[Paul Kienzle]

On Dec 4, 2005, at 2:28 AM, Mike Miller wrote:

> On Sat, 3 Dec 2005, Paul Kienzle wrote:
>
> > On Dec 1, 2005, at 9:58 AM, Mike Miller wrote:
> >
> > > On Thu, 1 Dec 2005, Mike Miller wrote:
> > > > In one Mersenne Twister FAQ I see this:
> > > >   2002-version mt19937ar.c has versions [0,1],[0,1),[0,1),(0,1) 
> > > > (32-bit
> > > >   precision), and [0,1) with 53-bit precision.
> > > > Which am I getting when I use rand.oct in Octave forge?  What I 
> > > > would really like to know are the largest and smallest values that 
> > > > can be produced both by rand.oct and by randn.oct.
> > > I should say that it's really one .oct file:  rand.oct.
> > > randn is derived from rand.oct.
> > > I took a look at Octave Forge's rand.cc but I can't find in it what 
> > > I need to know (because of my incompetence, obviously).
> >
> > The answer is in Makefile rather than rand.cc.  This sets the 
> > ALL_BITS flag for building randmtzig.c.  Here's the relevant code:
> >
> > /* generates a random number on (0,1) with 53-bit resolution */
> > USE_INLINE double randu53(void)
> > {
> >    const uint32_t a=randi32()>>5;
> >    const uint32_t b=randi32()>>6;
> >    return(a*67108864.0+b+0.4) * (1.0/9007199254740992.0);
> > }
>
>
> OK.  Thanks.  I guess "(0,1)" means that neither 0 nor 1 will ever be 
> generated as a random uniform, but does that last line of the code 
> tell me the largest and smallest possible values for a random uniform? 
>  I can't figure out what the largest and smallest numbers would be.

It does.  Expressed in powers of two the expression is:

   ( [0,2^27-1] * 2^26 + [0,2^26-1] + 0.4 ) / 2^53

So the min is 0.4/2^53 and the max is (2^53-1 + 0.4)/2^53 = 1 - 0.6/2^53

> Is it difficult to determine the range of possible values for the 
> normal deviates (randn with mean 0 and variance 1 in Octave Forge)?

Again, looking at the code:

#define ZIGGURAT_NOR_R 3.6541528853610088
#define ZIGGURAT_NOR_INV_R 0.27366123732975828
...
          /* As stated in Marsaglia and Tsang
           *
           * For the normal tail, the method of Marsaglia[5] provides:
           * generate x = -ln(U_1)/r, y = -ln(U_2), until y+y > x*x,
           * then return r+x. Except that r+x is always in the positive
           * tail!!!! Any thing random might be used to determine the
           * sign, but as we already have r we might as well use it
           *
           * [PAK] but not the bottom 8 bits, since they are all 0 here!
           */
          double xx, yy;
          do
            {
              xx = - ZIGGURAT_NOR_INV_R * log (RANDU);
              yy = - log (RANDU);
            }
          while ( yy+yy <= xx*xx);
          return (rabs&0x100 ? -ZIGGURAT_NOR_R-xx : ZIGGURAT_NOR_R+xx);

Given x^2 < 2y and the largest possible value of 2y is -2 
*log(min(RANDU)), that means |randn| is at most

        sqrt( -2*log(min(RANDU)) ) + 3.6541528853610088

Plugging in 0.4/2^53 for min(RANDU), randn is in (-12.332,12.332).  Not 
quite (-Inf,Inf) like it should be :-(

- Paul



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------