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