Random number generation open/closed interval

The function rand generates random numbers using open interval (0,1). How can I generate random values for partly closed or fully closed interval like (0,1] or [0,1] or [0,1).

Answer 1

[Jan]

Most RNG create chunks of 32 random bits. They are converted to 0 <= x < 1 DOUBLEs by e.g.:

a = bitshift(randi([0, 2^32-1]), -5);
b = bitshift(randi([0, 2^32-1]), -6);
d = (a * 67108864.0 + b) / 9007199254740992.0

For 0 <= x <= 1:

d = (a * 67108864.0 + b) / 9007199254740991.0

For 0 < x <= 1:

d = (1 + a * 67108864.0 + b) / 9007199254740992.0

A simple rejection mechanism is correct: Creating 0 <= x <= 1 and rejecting 0 and/or 1 conserves the normal distribution.

Another frequently used method:

a = randi([0, 2^32-1]);
b = bitshift(randi([0, 2^32-1]), -11);
d = (2097152.0 * a + b) / 9007199254740992.0;

This is slower if implemented in C with unsigned integers, because UINT32->DOUBLE is implemented in software, while INT32->DOUBLE is a processor command.

[Sakshi]

When I tried above commands, it gave maximum variable size allowed by the program is exceeded. Also how can I verify whether the result includes one of the boundary value since the process is random, it may not do so for a large number of random trials..

[Jan]

@Sakshi: Please post the command, which causes the error.

You can simply check, if the boundaries are reached by using 2^8 instead of 2^32...

[Mark Kittisopikul]

The arguments for randi are incorrect here. For scalar inputs, the first argument is imax and the second input, n, designates the size of the output matrix n x n. imin is implied to be zero. Thus we correct the above example as follows:

a = bitshift(randi(2^32-1), -5);
b = bitshift(randi(2^32-1), -6);

If you give a vector as the first argument then it specifies imin and imax. We can thus alternatively add brackets to fix the example.

a = bitshift(randi([0 2^32-1]), -5);
b = bitshift(randi([0 2^32-1]), -6);

[Walter Roberson]

randi with scalar first argument has an implied imin of 1, not of 0.

[Jan]

I've fixed the wrong input for randi now.

Answer 2

[Oleg Komarov]

rand generates on the closed interval [...]

What comes to my mind:

• to generate on (...)

rand(...)*(1-2*eps) + eps

• to generate on (...]

rand(...)*(1-eps) + eps

• to generate on [...)

rand(...)*(1-eps)

EDIT

It seems the example is wronlgy indicating the closed interval and some tests confirm the that is actually (0,1):

nnz(rand(1e8,1) == 0)
nnz(rand(1e8,1) == 1)

No matches.

[Sakshi]

Thank you for your answer. I appreciate the methods suggested by you. But I just want to confirm whether rand produces number in the open or closed interval because in the documentation,it is given as open interval. I quote the documentation statement below :

"r = rand(n) returns an n-by-n matrix containing pseudorandom values drawn from the standard uniform distribution on the open interval (0,1)."

[Oleg Komarov]

I guess we have to submit a documentation clarification:

"EXAMPLE 1:

Generate values from the uniform distribution on the interval [a, b]:

r = a + (b-a).*rand(100,1);"

[Walter Roberson]

It is more difficult to navigate the documentation for random numbers these days to find out what the boundaries of each one are.

I am sure that the default generator these days cannot generate 0 or 1 exactly. Or rather it is not possible internally for it to generate 1, and the rand routine deliberately rejects 0.

My tests in R2008b with about 2e9 random numbers is no 0's and no 1's... but of course I just might not have seen them.

[Jan]

@Walter, @Oleg: 2e9 and 1e8 random values are a really tiny subset to look for a specific number: If RAND replies a closed interval, the chance to get 0 or 1 is 2^-53.

[Peter Perkins]

Oleg, the "on" in that example is intentional. It is surprising but true that even though rand returns values strictly within the open interval (0,1), when you transform that interval to have arbitrary endpoints, rounding in double precision variously leads to closed, half-closed, or open intervals, depending on the length of the interval relative to the magnitudes of a and b.

Walter, you are correct that for the current generator choices, rand internally cannot create a 1, and for the generators that can (internally) create a 0, those values are (very infrequently) rejected. In the current release, anyways, the range of rand is given in the first sentence of its reference page.

[Walter Roberson]

In some previous releases it was documented that the range for some generators was [2^(-53), 1-2^(-53)] rather than documenting it as (0,1) . This was in the 2008a documentation for example.

Peter, you yourself indicated earlier that the default twister algorithm is NOT (0,1) .

http://www.mathworks.com/matlabcentral/answers/879-does-matlab-have-a-birthday-problem#answer_1262

in which you wrote,

Walter, the User Guide includes this: "mt19937ar: The Mersenne Twister, as described in [9], with Mersenne prime 2^19937-1. This is the generator documented at http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html. It has a period of 2^19937-1. Each U(0,1) value is created using two 32-bit integers from the generator; the possible values are all multiples of 2^-53 strictly within the interval (0,1)."

Notice in particular that this cannot generate anything on the interval [0, 2^(-53))

When I looked through the documentation of the specific ranges for each generator (this from before Mathworks removed the information from the reference pages where it belongs IMHO), I was unable to find any uniform distribution generator that was in fact dense over (0,1) .

Answer 3

[Charles]

I think that the easy solution would be to generate the regular random numbers, then set all values <= 0.5 to 0 and values > 0.5 to 1. That should be pretty random.

[Walter Roberson]

That would also not fit the requirement to generate over [0,1] or [0,1)

[Charles]

Oops, didn't really notice the restriction initially.