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