Random number generation open/closed interval
19 次查看(过去 30 天)
显示 更早的评论
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).
0 个评论
回答(3 个)
Jan
2011-8-7
编辑:Jan
2018-3-21
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.
5 个评论
Oleg Komarov
2011-8-6
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.
6 个评论
Peter Perkins
2011-8-8
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
2011-8-8
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) .
Charles
2011-8-7
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.
2 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Random Number Generation 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!