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

[Mike Miller]

On Sun, 4 Dec 2005, Paul Kienzle wrote:

> 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

[snip]

> 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).


Thanks, Paul.  Sorry that I needed it to be so explicit.  I like the 
choices you've made.  It is good that we don't have to check that rand is 
zero or one in cases like 1/rand() or 1/(1-rand()).  The range on the 
normals is also very good for most uses.

If we want more of a range on randn or we want numbers closer to zero for 
rand, I suppose we could use some kind of trick.  I designed this to 
produce N random uniforms that hit every value between 0 and 1 at 
increments of 10^-20 starting at 5 x 10^-21 and ending at 1 minus 5 x 
10^-21 :

floor(100*rand(N,10)) * 10.^(-2*[1:10]') + 5*10^-21

But I wonder if that fails in reality because machine precision doesn't 
allow that many digits.

Using a similar idea, it is also possible to get A and B separately for 
uniform (0,1) random numbers of the form A x 10^-B and have any desired 
upper limit for B.  That could allow extension of the range of the normal 
deviates to include much larger possible values.  It would be much slower 
than the rand.oct function, but there might be some uses for a function 
that does that.  Maybe not.

Mike



-------------------------------------------------------------
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
-------------------------------------------------------------