[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

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

From: Mike Miller
Subject: Re: Mersenne Twister range in rand.oct and randn.oct
Date: Mon, 5 Dec 2005 00:17:13 -0600 (CST)

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


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.


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

Octave's home on the web:
How to fund new projects:
Subscription information:

reply via email to

[Prev in Thread] Current Thread [Next in Thread]