Random number generators
Marina
I have a question regarding the subroutine RANDOM_NUMBER,
does it generate a list of uniformly distributed numbers, between 0
and 1, and an average of 0.5?
Which other intrinsic random number generating subroutines are
available for Fortran 90, for other distributions?
(I found that the function rand() doesn't work in my compiler, for
example.)
Thanks,
Marina
steve
>
> I have a question regarding the subroutine RANDOM_NUMBER,
> does it generate a list of uniformly distributed numbers, between 0
> and 1, and an average of 0.5?
The description for RANDOM_NUMBER is given in 13.14.85 of the
Fortran 95 standard. It states
Description. Returns one pseudorandom number or an array of
pseudorandom numbers from the uniform distribution over the
range 0 <= x < 1.
Note, it says nothing about an average of 0.5. If you only draw 2
values from the distribution is highly unlikely that the average is
0.5. OTOH, if you draw hundred or thousands or millions of samples,
then the average will be near 0.5.
> Which other intrinsic random number generating subroutines are
> available for Fortran 90, for other distributions?
The Fortran standard requires that a comforming compiler must
supply RANDOM_NUMBER. Vendors can choose to supply additional
random number generators, so check your documentation.
Note, there are numerous Fortran implementations of PRNG. Try
searching on 'metcalf KISS' and 'ranlux' to name just two. There
are others. As always, you should test whatever PRNG you choose
to determine if it is adequate for your application.
--
steve
frank
>> I have a question regarding the subroutine RANDOM_NUMBER, does it
>> generate a list of uniformly distributed numbers, between 0 and 1, and
>> an average of 0.5?
>
> The description for RANDOM_NUMBER is given in 13.14.85 of the Fortran 95
> standard. It states
>
> Description. Returns one pseudorandom number or an array of
> pseudorandom numbers from the uniform distribution over the range 0
> <= x < 1.
>
> Note, it says nothing about an average of 0.5. If you only draw 2
> values from the distribution is highly unlikely that the average is 0.5.
> OTOH, if you draw hundred or thousands or millions of samples, then
> the average will be near 0.5.
Isn't there an asymmetry in the unit interval though as to which endpoint
is included? So if there's N outcomes on one side of .5 there would be N
+-1 on the other.
--
frank
"Guns: yes, they are harmful."
steve
Can you restate your question without a little more detail? For
example,
for N=4, random_number in gfortran returns 0.998, 0.567, 0.966, 0.748.
If I draw 10000 numbers and count the number in [0,0.5) and [0.5,1),
the
outcome is 4926 and 5074. If I draw 10000000 numbers, I get 4999679
and 50000321. If I reseed the rng with a different seed, these
numbers
change.
--
steve
Arjen Markus
Basically, from a uniform random number you can derive a random number
with any distribution. It is a matter of using the inverse cumulative
distribution function or one of several standard procedures.
You may find the code and documentation by Alan Miller useful -
a zip file miller_all.zip at http://groups.google.com/group/gg95.
Regards,
Arjen
robin
"Marina" <levin....@gmail.com> wrote in message
news:94929b54-dacb-4d20...@m13g2000vbf.googlegroups.com...
Rand is not a standard function.
Gordon Sande
> On Nov 16, 1:58�pm, Marina <levin.mar...@gmail.com> wrote:
>>
>> I have a question regarding the subroutine RANDOM_NUMBER,
>> does it generate a list of uniformly distributed numbers, between 0
>> and 1, and an average of 0.5?
>
> The description for RANDOM_NUMBER is given in 13.14.85 of the
> Fortran 95 standard. It states
>
> Description. Returns one pseudorandom number or an array of
> pseudorandom numbers from the uniform distribution over the
> range 0 <= x < 1.
>
> Note, it says nothing about an average of 0.5. If you only draw 2
> values from the distribution is highly unlikely that the average is
> 0.5. OTOH, if you draw hundred or thousands or millions of samples,
> then the average will be near 0.5.
The standard is even less demanding than that as it says absolutely
nothing about the quality of the random numbers. If it always
returned 0.33 it would be standard conforming and under some definitions
of randomness as surprise it would even be of high quality because
such a sequence would be highly surprising. ;-) It would surely draw
many user complaints about the poor quality of implementation.
There is a common issue about whether the sequence should always
default to the same one for ease of debugging or should it be "random".
steve
> There is a common issue about whether the sequence should always
> default to the same one for ease of debugging or should it be "random".
I'm well aware of this issue because I'm the person that decided
gfortran should default to a given sequence if one does not call
random_seed. I made this decision precisely as you state "for the
ease of debugging."
--
steve
Gordon Sande
The only folks who I have heard question such a choice are beginners
who are surprised when the random sequence is always the same on their
successive runs (that are often for debugging :-( ). It was not really
fair to call it an issue even if it does confuse a few beginners.
nm...@cam.ac.uk
Then try me. I am one of the experts who takes the converse view,
on the grounds of statistical validity. In particular, using the
same numbers to test two related analyses gives erroneous results
if you compare the analyses statistically assuming independence.
A classic error of scientists who ought to know better.
However, I accept that it is a matter of judgement which choice will
lead to more wrong answers, overall.
Regards,
Nick Maclaren.
glen herrmannsfeldt
I think I agree with this in the case of no call to RANDOM_SEED.
It would, however, be nice if a call to RANDOM_SEED() with no
arguments did choose a good seed based on something, such as the
current time/date. (That is, for example, the way that many early
BASIC systems did it, though it is hard to say that there is any
standard for BASIC.)
While one can use DATE_AND_TIME to extract the date/time and then
use that as a seed, there is no way to know what constitutes
a good seed. Duplicating the values as many times as needed
to fill up the seed array may or may not supply a good seed.
It would be nice to have a way to generate a good seed based
on a single default kind integer. Again, duplicating to fill
the seed array is likely not the best choice.
-- glen
Gordon Sande
> In article <2009111711204716807-gsande@worldnetattnet>,
> Gordon Sande <g.s...@worldnet.att.net> wrote:
>> On 2009-11-17 11:13:29 -0400, steve <kar...@comcast.net> said:
>>>
>>>> There is a common issue about whether the sequence should always
>>>> default to the same one for ease of debugging or should it be "random".
>>>
>>> I'm well aware of this issue because I'm the person that decided
>>> gfortran should default to a given sequence if one does not call
>>> random_seed. I made this decision precisely as you state "for the
>>> ease of debugging."
>>
>> The only folks who I have heard question such a choice are beginners
>> who are surprised when the random sequence is always the same on their
>> successive runs (that are often for debugging :-( ). It was not really
>> fair to call it an issue even if it does confuse a few beginners.
>
> Then try me. I am one of the experts who takes the converse view,
> on the grounds of statistical validity. In particular, using the
> same numbers to test two related analyses gives erroneous results
> if you compare the analyses statistically assuming independence.
> A classic error of scientists who ought to know better.
The usual professional advice is that such experiments need to record
both the random number generater used as well as the seed value so
that rerunning would be possible. Implicitly this says that differing
experiments should have different seeds.
> However, I accept that it is a matter of judgement which choice will
> lead to more wrong answers, overall.
Too many folks act like beginners. A kinder way of saying this is that
people are specialized but problems are not so not all problems are solved
as professionally as one might hope.
> Regards,
> Nick Maclaren.
steve
> steve <kar...@comcast.net> wrote:
> > On Nov 17, 6:41?am, Gordon Sande <g.sa...@worldnet.att.net> wrote:
> >> There is a common issue about whether the sequence should always
> >> default to the same one for ease of debugging or should it be "random".
> > I'm well aware of this issue because I'm the person that decided
> > gfortran should default to a given sequence if one does not call
> > random_seed. I made this decision precisely as you state "for the
> > ease of debugging."
>
> I think I agree with this in the case of no call to RANDOM_SEED.
> It would, however, be nice if a call to RANDOM_SEED() with no
> arguments did choose a good seed based on something, such as the
> current time/date. (That is, for example, the way that many early
> BASIC systems did it, though it is hard to say that there is any
> standard for BASIC.)
I respectfully disagree with your second sentence because person
will need to make 2 calls to RANDOM_SEED(); one to reseed the PRNG
and one to record the set of seeds.
> While one can use DATE_AND_TIME to extract the date/time and then
> use that as a seed, there is no way to know what constitutes
> a good seed. Duplicating the values as many times as needed
> to fill up the seed array may or may not supply a good seed.
>
> It would be nice to have a way to generate a good seed based
> on a single default kind integer. Again, duplicating to fill
> the seed array is likely not the best choice.
gfortran's random_seed() in principle has a mechanism in
place to help prevent a user from making a bad choice of
seeds. Whether the mechanism is effective, I do not know.
--
steve
e p chandler
> Hi,
>
> I have a question regarding the subroutine RANDOM_NUMBER,
> does it generate a list of uniformly distributed numbers, between 0
> and 1, and an average of 0.5?
The average is a property of the distribution, not the generator. The
Fortran standard says nothing about its qualidy, just its interface.
> Which other intrinsic random number generating subroutines are
> available for Fortran 90, for other distributions?
Only uniform(0,1) is intrinsic. Others are gernrally built from the
uniform distribution.
There is an ENORMOUS literature on this subject. For example there are
standard texts by Knuth, Gentle, Law and Kelton, and many others. Also
textbooks on mathematical statistics cover this subject
(transfomration of random variables) in detail.
> (I found that the function rand() doesn't work in my compiler, for example.)
It's not part of the standard. Vendors offer many RNGs as extensions.
My general advice to you is:
1. Determine what distribution you want to model.
2. Consult a specialist in random number generation and statistics.
3. Test, test and test.
---- elliot
glen herrmannsfeldt
(snip, I wrote)
>> It would be nice to have a way to generate a good seed based
>> on a single default kind integer. ?Again, duplicating to fill
>> the seed array is likely not the best choice.
> gfortran's random_seed() in principle has a mechanism in
> place to help prevent a user from making a bad choice of
> seeds. Whether the mechanism is effective, I do not know.
How does it know if it is a user supplied seed value, or a
previously stored seed from RANDOM_SEED? I think what I would
most like to see is separate ways for restoring a previous seed
and specifying a new seed. For a new seed, the routine should
generate a good seed based on user supplied data. When restoring
a seed, it should not be changed. With the assumption that all
seed values are possible, there is no way to tell the difference.
I suppose one could put a CRC on the stored seed, and test it
on reloading. For a given CRC size, that makes it very unlikely
to be accidentally specified by the user.
-- glen
Maciej Marek
> The standard is even less demanding than that as it says absolutely
> nothing about the quality of the random numbers. If it always
> returned 0.33 it would be standard conforming and under some definitions
> of randomness as surprise it would even be of high quality because
> such a sequence would be highly surprising. ;-) It would surely draw
> many user complaints about the poor quality of implementation.
That reminds me a Dilbert strip:
- And here's our random number generator.
- 9,9,9,9,9,9,9,9....
- Are you sure that's random?
- That's the problem with randomness. You can never be sure.
Regards,
Maciej Marek
steve
wrote:
I don't understand your question. A call to random_seed()
manipulates the seeds and then initializes the prng. A call
to random_number() simply uses the current values.
glen herrmannsfeldt
(snip)
> I don't understand your question. A call to random_seed()
> manipulates the seeds and then initializes the prng. A call
> to random_number() simply uses the current values.
There are two ways random_seed(put=...) can be used:
1) Initialize the seed based on some user supplied data.
(Constant, time of day, phase of moon, etc.)
2) Restore a seed previously retrieved using random_seed(get=...)
In the first case, the supplied data should be used to generate
a good quality seed. In the second case, the data should be used
as-is, such that a previous seed is restored exactly.
I have done problems before where I needed to restart part way
though, and verify that the results were the same as previous.
(That was in R, which does have separate restore seed and initial
seed options.)
How do you determine if case 1 or case 2 applies?
-- glen
steve
It appears as if you did not read the message in the fortran @ gnu
mailing list archive that I pointed you at with URL. Your question
doesn't make sense (to me). Random_seed() scrambles the seeds
when called with PUT= and it unscrambles the seeds when called
with GET=. If called with no argument, the default set of seeds
isused. You can use GET= and PUT= with the default set of seeds,
but the scrambling and unscrambling of the seeds still occurs.
For more details, http://gcc.gnu.org/bugzilla/show_bug.cgi?id=32812
glen herrmannsfeldt
(snip, I wrote)
>> There are two ways random_seed(put=...) can be used:
>> 1) Initialize the seed based on some user supplied data.
>> ? ?(Constant, time of day, phase of moon, etc.)
>> 2) Restore a seed previously retrieved using random_seed(get=...)
>> In the first case, the supplied data should be used to generate
>> a good quality seed. ?In the second case, the data should be used
>> as-is, such that a previous seed is restored exactly.
(snip)
> It appears as if you did not read the message in the fortran @ gnu
> mailing list archive that I pointed you at with URL. Your question
> doesn't make sense (to me). Random_seed() scrambles the seeds
> when called with PUT= and it unscrambles the seeds when called
> with GET=. If called with no argument, the default set of seeds
> isused. You can use GET= and PUT= with the default set of seeds,
> but the scrambling and unscrambling of the seeds still occurs.
> For more details, http://gcc.gnu.org/bugzilla/show_bug.cgi?id=32812
To quote from the mentioned site:
"This doesn't make miracles (i.e. provide you with a good
seed when you input a particularly poor one), but at least
it makes using the VALUES of DATE_AND_TIME less frustrating
(by generating visibly different streams of PRN)."
It would be nice to have high quality seeds for all possible
values of DATE_AND_TIME when used to generate seeds.
Different generators have different requirements on what
makes a good seed. I suppose I agree that for a sufficient
amount of scrambling it works, assuming the scrambling
can be inverted.
-- glen
steve
> steve <kar...@comcast.net> wrote:
>
> (snip, I wrote)
>
> >> There are two ways random_seed(put=...) can be used:
> >> 1) Initialize the seed based on some user supplied data.
> >> ? ?(Constant, time of day, phase of moon, etc.)
> >> 2) Restore a seed previously retrieved using random_seed(get=...)
> >> In the first case, the supplied data should be used to generate
> >> a good quality seed. ?In the second case, the data should be used
> >> as-is, such that a previous seed is restored exactly.
>
> (snip)
>
> > It appears as if you did not read the message in the fortran @ gnu
> > mailing list archive that I pointed you at with URL. Your question
> > doesn't make sense (to me). Random_seed() scrambles the seeds
> > when called with PUT= and it unscrambles the seeds when called
> > with GET=. If called with no argument, the default set of seeds
> > isused. You can use GET= and PUT= with the default set of seeds,
> > but the scrambling and unscrambling of the seeds still occurs.
>
> To quote from the mentioned site:
>
> "This doesn't make miracles (i.e. provide you with a good
> seed when you input a particularly poor one), but at least
> it makes using the VALUES of DATE_AND_TIME less frustrating
> (by generating visibly different streams of PRN)."
>
> It would be nice to have high quality seeds for all possible
> values of DATE_AND_TIME when used to generate seeds.
So, what bullet-proof, general-purpose, algorithm to you recommend?
> Different generators have different requirements on what
> makes a good seed. I suppose I agree that for a sufficient
> amount of scrambling it works, assuming the scrambling
> can be inverted.
Of course, the scrambling has to be invertible by some
means; otherwise, PUT= followed by an immediate GET=
wouldn't return the desired result.
--
steve
glen herrmannsfeldt
(snip, I quoted and wrote)
>> To quote from the mentioned site:
>> ? "This doesn't make miracles (i.e. provide you with a good
>> ? ?seed when you input a particularly poor one), but at least
>> ? ?it makes using the VALUES of DATE_AND_TIME less frustrating
>> ? ?(by generating visibly different streams of PRN)."
>>
>> It would be nice to have high quality seeds for all possible
>> values of DATE_AND_TIME when used to generate seeds.
> So, what bullet-proof, general-purpose, algorithm to you recommend?
Well, if I wrote the standard I would have had different options
for restoring and initial seeding. A little late for that, though.
Otherwise, my favorite bit scrambler is the CRC32 LFSR. Some
combination of CRC32 across the words, XORing the result with
each word as it goes along. I think it isn't hard to invert
that, and XOR should be fairly fast. Byte at a time CRC32
requires a 256 word table of constants. Bit at a time is
slower, and doesn't need the table.
>> Different generators have different requirements on what
>> makes a good seed. ?I suppose I agree that for a sufficient
>> amount of scrambling it works, assuming the scrambling
>> can be inverted.
> Of course, the scrambling has to be invertible by some
> means; otherwise, PUT= followed by an immediate GET=
> wouldn't return the desired result.
Note that the standard doesn't require PUT followed by GET
to return the same value, but does require the same random
numbers be generated if either value is used with PUT.
(That is, K may be larger than the length actually used
for the seed.)
-- glen
robin
Sent: Wednesday, 18 November 2009 2:07 AM
Subject: Re: Random number generators
>On Nov 17, 12:53 am, "robin" <robi...@bigpond.com> wrote:
>> "Marina" <levin.mar...@gmail.com> wrote in message
>
>> news:94929b54-dacb-4d20...@m13g2000vbf.googlegroups.com...
>> | Hi,
>> |
>> | I have a question regarding the subroutine RANDOM_NUMBER,
>> | does it generate a list of uniformly distributed numbers, between 0
>> | and 1, and an average of 0.5?
>> |
>> | Which other intrinsic random number generating subroutines are
>> | available for Fortran 90, for other distributions?
>> | (I found that the function rand() doesn't work in my compiler, for
>> | example.)
>
>> Rand is not a standard function.
>Sure it is.
> 7.20.2.1 The rand function
> Synopsis
>1 #include <stdlib.h>
> int rand(void);
No it isn''t. That's C, not the Fortran under discussion.
robin
If that behaviour is required, then a call to randon_seed will obtain it.
There's no point in making the default do exactly the same thing.
steve
> steve <kar...@comcast.net> wrote:
>
> (snip, I quoted and wrote)
>
> >> To quote from the mentioned site:
> >> ? "This doesn't make miracles (i.e. provide you with a good
> >> ? ?seed when you input a particularly poor one), but at least
> >> ? ?it makes using the VALUES of DATE_AND_TIME less frustrating
> >> ? ?(by generating visibly different streams of PRN)."
>
> >> It would be nice to have high quality seeds for all possible
> >> values of DATE_AND_TIME when used to generate seeds.
> > So, what bullet-proof, general-purpose, algorithm to you recommend?
>
> Well, if I wrote the standard I would have had different options
> for restoring and initial seeding. A little late for that, though.
You're evading my question. :)
You opine on generating high quality seeds for all possible
values of DATE_AND_TIME, and ask you for your bullet-proof,
generate-purpose algorithm to achieve this aim. Some how you
mutate that request into if you wrote the standard, it would
have been done differently.
> Otherwise, my favorite bit scrambler is the CRC32 LFSR. Some
> combination of CRC32 across the words, XORing the result with
> each word as it goes along. I think it isn't hard to invert
> that, and XOR should be fairly fast. Byte at a time CRC32
> requires a 256 word table of constants. Bit at a time is
> slower, and doesn't need the table.
Well, I suppose this is start. Now, take the values from
DATE_AND_TIME
and provide the algorithm that generates a high quality set of seeds.
--
steve
Richard Maine
[about the Fortran standard's random number generator spec]
> Well, if I wrote the standard I would have had different options
> for restoring and initial seeding.
Of course, that assumes that you sucessfully got your intention into the
written words. :-(
The person who did mostly write the random number part of the standard
claimed that he had intended to specify things about the issue of when
the initial seed was the same for each run and when it was different,
but he accidentally failed to put his intent into the written words. I
take that as a precautionary tale worth noting - something about good
ideas not being worth much if you don't manage to communicate them.
(No, I'm not meaning to imply anything about posts or people in this
thread; it is just a more general precautionary tale that the thread
brought to mind.)
I don't actually recall what the intended specification was because the
answer to an interp request basically said that his unstated (and thus
undiscussed in committee) intent did not matter and that the words of
the standard ruled instead. Those words are the ones that allowed the
choice to be vendor-dependent. Different vendors had already implemented
different choices in released compilers before the interpretation
request came back.
--
Richard Maine | Good judgment comes from experience;
email: last name at domain . net | experience comes from bad judgment.
domain: summertriangle | -- Mark Twain
David Jones
> On Nov 16, 6:04 pm, frank <fr...@example.invalid> wrote:
>> On Mon, 16 Nov 2009 16:05:14 -0800, steve wrote:
>>>> I have a question regarding the subroutine RANDOM_NUMBER, does it
>>>> generate a list of uniformly distributed numbers, between 0 and 1,
>>>> and an average of 0.5?
>>
>>> The description for RANDOM_NUMBER is given in 13.14.85 of the
>>> Fortran 95 standard. It states
>>
>>> Description. Returns one pseudorandom number or an array of
>>> pseudorandom numbers from the uniform distribution over the range 0
>>> <= x < 1.
>>
>>> Note, it says nothing about an average of 0.5. If you only draw 2
>>> values from the distribution is highly unlikely that the average is
>>> 0.5. OTOH, if you draw hundred or thousands or millions of samples,
>>> then
>>> the average will be near 0.5.
>>
>> Isn't there an asymmetry in the unit interval though as to which
>> endpoint is included? So if there's N outcomes on one side of .5
>> there would be N +-1 on the other.
>
> Can you restate your question without a little more detail?
The point centres on the explicit statement of the range in the standard as 0<= x < 1. Thus an exact value of zero is allowed as an output. while an exact value of one is not. So there is asymmetry on this point. For some uses, the question of whether or not exact zeroes and ones are ever returned is important and is sometimes addressed by testing for such values and re-generating another as necessary. Obviously such testing should be avoided if the zero or one case is guaranteed not to occur. Random number generators are such as to naturally include zero as a possible output, and to avoid this would involve slowing the code down.
A description of a RNG should explicity state the range and discretisation of the values that might be reurned, so that it would be possible to see whether just rejecting zero values would leave the set of allowed values symmetric about zero.
David Jones
glen herrmannsfeldt
> steve wrote:
>> On Nov 16, 6:04 pm, frank <fr...@example.invalid> wrote:
>>> Isn't there an asymmetry in the unit interval though as to which
>>> endpoint is included? So if there's N outcomes on one side of .5
>>> there would be N +-1 on the other.
>> Can you restate your question without a little more detail?
> <snip>
> The point centres on the explicit statement of the range in the
> standard as 0<= x < 1. Thus an exact value of zero is allowed
> as an output. while an exact value of one is not. So there is
> asymmetry on this point. For some uses, the question of whether
> or not exact zeroes and ones are ever returned is important and
> is sometimes addressed by testing for such values and re-generating
> another as necessary. Obviously such testing should be avoided if
> the zero or one case is guaranteed not to occur. Random number
> generators are such as to naturally include zero as a possible
> output, and to avoid this would involve slowing the code down.
Consider a true random number generator (not a pseudo-random
number generator, as is usually used) generating an IEEE double
with a 53 bit significand generating uniformly 2**53 possible
values, 0<=x<1. The expected standard deviation of the mean
after averaging N values is 1/sqrt(N). To notice a deviation
of 1 in 2**53, the offset due to the asymmetry, requires N of
about 2**104, or about 1e31. At 1GHz, that takes 1e22 seconds,
or, if I did it right, about 6e14 years. Good PRNG generate
enough additional bits internally that it will be about as
hard to notice.
> A description of a RNG should explicity state the range and
> discretisation of the values that might be reurned, so that
> it would be possible to see whether just rejecting zero values
> would leave the set of allowed values symmetric about zero.
If processor speeds double every year for the next 100 years
then you might have a chance to notice. If by then the generators
switch to a longer return value then you won't.
-- glen
Ron Shepard
<ab2053c6-0fcf-4c46...@i12g2000prg.googlegroups.com>,
steve <kar...@comcast.net> wrote:
> Note, there are numerous Fortran implementations of PRNG. Try
> searching on 'metcalf KISS' and 'ranlux' to name just two. There
> are others. As always, you should test whatever PRNG you choose
> to determine if it is adequate for your application.
Most of my applications involving random numbers do not place stringent
requirements on the PRNG. Almost anything would work in these
situations. If I do something more sophisticated that does require long
cycles, etc., then I would use ranlux or something similar where I knew
the algorithm and could verify that it met whatever requirements where
necessary.
The fortran standard does not specify very much about the intrinsic
PRNG, RANDOM_NUMBER(). But if it did, then what kind of requirements
should it meet? What are the most commonly needed requirements? Should
it be required to have a cycle that is at least as long as the number of
floating point numbers between 0 and 1, or should it be allowed to be
restricted only to the number of values with a fixed exponent? What
kind of restrictions might be imposed on correlations between the
numbers that are generated?
$.02 -Ron Shepard
Ron Shepard
frank <fr...@example.invalid> wrote:
> Isn't there an asymmetry in the unit interval though as to which endpoint
> is included? So if there's N outcomes on one side of .5 there would be N
> +-1 on the other.
There are more floating point values between 0 and .5 than there are
between .5 and 1. It is not just a difference of +-1 value. Should all
such values occur in the pseudorandom sequence, or only a subset of such
values?
$.02 -Ron Shepard
Barry W Brown
Our library, RANDLIB has routines for uniform random numbers and
generators for all of
the common statistical distributions. It is available free for the
source download
from
http://biostatistics.mdanderson.org/SoftwareDownload
glen herrmannsfeldt
> There are more floating point values between 0 and .5 than there are
> between .5 and 1. It is not just a difference of +-1 value. Should all
> such values occur in the pseudorandom sequence, or only a subset of such
> values?
Interesting question. In the days of 32 bit generators, ones I
knew just took 24 bits of the integer, added on the appropriate
exponent, maybe add zero to normalize it, and return. For IEEE
you can't do that, as it has the hidden one. Since most have
more state bits than significand bits it wouldn't be hard to
return the extra bits.
-- glen
David Flower
> does it generate a list of uniformly distributed numbers, between 0
> and 1, and an average of 0.5?
>
> available for Fortran 90, for other distributions?
> (I found that the function rand() doesn't work in my compiler, for
> example.)
>
> Marina
Posters may be interested in the following reference:
A.C.M. Trans. Math, Software, 5, #2, 132 (1979) by Linus Schrage
It is a random number generator that is machine and compiler
independent, useful for testing.
I cannot vouch for its efficiency or randomness.
If anyone is interested I will post a copy of the source code
Dave Flower
robin
>The point centres on the explicit statement of the range in the standard as 0<= x < 1. Thus an exact value of zero is
>allowed as an output. while an exact value of one is not. So there is asymmetry on this point. For some uses, the
>question of whether or not exact zeroes and ones are ever returned is important and is sometimes addressed by testing
>for such values and re-generating another as necessary. Obviously such testing should be avoided if the zero or one
>case is guaranteed not to occur. Random number generators are such as to naturally include zero as a possible output,
>and to avoid this would involve slowing the code down.
But not by much, and probably wouldn't be noticed.
robin
| Ron Shepard <ron-s...@nospam.comcast.net> wrote:
|
| > There are more floating point values between 0 and .5 than there are
| > between .5 and 1. It is not just a difference of +-1 value. Should all
| > such values occur in the pseudorandom sequence, or only a subset of such
| > values?
|
| Interesting question. In the days of 32 bit generators, ones I
| knew just took 24 bits of the integer, added on the appropriate
| exponent, maybe add zero to normalize it, and return. For IEEE
| you can't do that, as it has the hidden one.
Since you'd be using integer arithmetic to prepare the FPN,
the question of IEEE doesn't come into it.
And that applies to both cases.
robin
news:8b69b09f-fd08-44a2...@z41g2000yqz.googlegroups.com...
>Posters may be interested in the following reference:
>A.C.M. Trans. Math, Software, 5, #2, 132 (1979) by Linus Schrage
It's a bit ancient.
Those by George Marsaglia are not only portable,
they also have extremely long periods.
His RNGs include 32-bit generators and 64-bit generators.
nm...@cam.ac.uk
32-bit generators should never be used in any simulation which
uses more than a million numbers in total. To a first approximation,
ALWAYS use 64-bit ones.
Not all of Marsaglia's generators are free from 32-bit defects, even
when they are used in 64-bit forms, though most of them are good or
very good. Their actual code isn't always very portable, but that's
easy to clean up.
Regards,
Nick Maclaren.
glen herrmannsfeldt
(snip)
> 32-bit generators should never be used in any simulation which
> uses more than a million numbers in total. To a first approximation,
> ALWAYS use 64-bit ones.
All the ones I knew from 40 years ago were 32 bit (or 31 bit).
For many problems, that is probably enough. I have noticed that
the games that come with Windows have bad RNGs. Of course closed
source so we really don't know how they do it.
> Not all of Marsaglia's generators are free from 32-bit defects, even
> when they are used in 64-bit forms, though most of them are good or
> very good. Their actual code isn't always very portable, but that's
> easy to clean up.
-- glen
frank
That's on one end of the spectrum. On the other is me, PC Joe in single
precision. I've worked a lot with this material, for example when I was
working up different sorts.
Things look pretty glompy down here:
program rand
implicit integer(a-t)
integer, parameter :: r = 10
real, dimension(1:r) :: myarray
real u
integer m
call init_seed
do i =1,r
call random_number(u)
myarray(i) = u
end do
m = 2**30 * 1.9
print *, "myarray is ", myarray
print *, "m is ", m
contains
subroutine init_seed()
integer :: n, ival(8), v(3), i
integer, allocatable :: seed(:)
call date_and_time(values=ival)
v(1) = ival(8) + 2048*ival(7)
v(2) = ival(6) + 64*ival(5) ! value(4) isn't really 'random'
v(3) = ival(3) + 32*ival(2) + 32*8*ival(1)
call random_seed(size=n)
allocate(seed(n))
call random_seed() ! Give the seed an implementation-dependent kick
call random_seed(get=seed)
do i=1, n
seed(i) = seed(i) + v(mod(i-1, 3) + 1)
enddo
call random_seed(put=seed)
deallocate(seed)
end subroutine
end program rand
! gfortran ran1.f90 -Wall -Wextra -o out
dan@dan-desktop:~/source$ gfortran ran1.f90 -Wall -Wextra -o out
dan@dan-desktop:~/source$ ./out
myarray is 0.92712229 0.85859424 0.79266590
0.57997841 0.31713003 0.38777798 0.23500705
0.40011257 0.90164834 0.72985947
m is 2040109440
dan@dan-desktop:~/source$
Since 2**31 overflowed, I think m is close to maybe 2.1 billion outcomes
to expect on the unit interval. If m =~ 2n +-1, then n is on the order
of one billion, which was only a large number while the great William
Proxmire was alive.
>
>> A description of a RNG should explicity state the range and
>> discretisation of the values that might be reurned, so that it would be
>> possible to see whether just rejecting zero values would leave the set
>> of allowed values symmetric about zero.
>
> If processor speeds double every year for the next 100 years then you
> might have a chance to notice. If by then the generators switch to a
> longer return value then you won't.
Glen, I would almost be scared if our processor speeds trended along the
same curve for the next forty years.
--
frank
"Guns: yes, they are harmful."
frank
>> Isn't there an asymmetry in the unit interval though as to which
>> endpoint is included? So if there's N outcomes on one side of .5 there
>> would be N +-1 on the other.
>
> Can you restate your question without a little more detail? For
> example,
> for N=4, random_number in gfortran returns 0.998, 0.567, 0.966, 0.748.
> If I draw 10000 numbers and count the number in [0,0.5) and [0.5,1), the
> outcome is 4926 and 5074. If I draw 10000000 numbers, I get 4999679 and
> 50000321. If I reseed the rng with a different seed, these numbers
> change.
At least conceptually I think it's important to get the [ ) [ )
correct on this on the test, as you have it. So the test would be if
(x.gt. (.5)).
steve
With your suggested test, the intervals are [0,0.5] and (0.5,1).
See Ron Shepard's post with respect to finite-precision floating point
numbers.
--
steve
Gib Bogle
How many more?
Gib Bogle
> In article <%ShNm.55991$ze1....@news-server.bigpond.net.au>,
> robin <rob...@bigpond.com> wrote:
>> "David Flower" <DavJF...@AOL.COM> wrote in message
>> news:8b69b09f-fd08-44a2...@z41g2000yqz.googlegroups.com...
>>
>>> Posters may be interested in the following reference:
>>> A.C.M. Trans. Math, Software, 5, #2, 132 (1979) by Linus Schrage
>> It's a bit ancient.
>> Those by George Marsaglia are not only portable,
>> they also have extremely long periods.
>>
>> His RNGs include 32-bit generators and 64-bit generators.
>
> 32-bit generators should never be used in any simulation which
> uses more than a million numbers in total.
Never say never. There are many (probably most, in fact) applications
where it will make absolutely no significant difference to the results
if the RNG repeats the cycle.
Ron Shepard
Gib Bogle <bo...@ihug.too.much.spam.co.nz> wrote:
About half of the exponent range.
In binary floating point, are the same number of (equally spaced)
values between .25 and .50 as there are between .50 and 1.0. Then
there are that same number of values between .125 and .25, and that
same number of values again between .0625 and .125, and so on, all
the way down to the minimum exponent value. In most formats, there
are about the same number of exponents less than zero as there are
greater than zero, so about half of the exponent range corresponds
to floating point values that are between 0 and .5, and only a
single exponent value corresponds to the floating point values
between .5 and 1.
For a given exponent, should a PRNG return all of the associated
floating point values with equal probability, or only some of them?
$.02 -Ron Shepard
nm...@cam.ac.uk
Gib Bogle <bo...@ihug.too.much.spam.co.nz> wrote:
>>
>> 32-bit generators should never be used in any simulation which
>> uses more than a million numbers in total.
>
>Never say never. There are many (probably most, in fact) applications
>where it will make absolutely no significant difference to the results
>if the RNG repeats the cycle.
Eh? Why does that contradict my statement? If you can produce an
example which uses more than a million numbers and where an almost
arbitrarily grotty generator will definitely not give misleading
results, I should be interested to hear of it.
Yes, there are a zillion applications where the results are nonsense
anyway, a zillion others which use only a few random numbers, and
so on. But those are out of order.
Regards,
Nick Maclaren.
robin
news:ron-shepard-C17E...@forte.easynews.com...
| In binary floating point, are the same number of (equally spaced)
| values between .25 and .50 as there are between .50 and 1.0. Then
| there are that same number of values between .125 and .25, and that
| same number of values again between .0625 and .125, and so on, all
| the way down to the minimum exponent value. In most formats, there
| are about the same number of exponents less than zero as there are
| greater than zero, so about half of the exponent range corresponds
| to floating point values that are between 0 and .5, and only a
| single exponent value corresponds to the floating point values
| between .5 and 1.
|
| For a given exponent, should a PRNG return all of the associated
| floating point values with equal probability, or only some of them?
Probably a case for integer RNGs.
robin
| In article <%ShNm.55991$ze1....@news-server.bigpond.net.au>,
| robin <rob...@bigpond.com> wrote:
| >"David Flower" <DavJF...@AOL.COM> wrote in message
| >news:8b69b09f-fd08-44a2...@z41g2000yqz.googlegroups.com...
| >
| >>Posters may be interested in the following reference:
| >
| >>A.C.M. Trans. Math, Software, 5, #2, 132 (1979) by Linus Schrage
| >
| >It's a bit ancient.
| >Those by George Marsaglia are not only portable,
| >they also have extremely long periods.
| >
| >His RNGs include 32-bit generators and 64-bit generators.
|
| 32-bit generators should never be used in any simulation which
| uses more than a million numbers in total.
George's 32-bit RNG KISS has a period > 10**35 numbers.
He has put his RNGs through exhaustive tests.
| To a first approximation,
| ALWAYS use 64-bit ones.
|
| Not all of Marsaglia's generators are free from 32-bit defects, even
| when they are used in 64-bit forms, though most of them are good or
| very good. Their actual code isn't always very portable, but that's
| easy to clean up.
Don't know about his early ones, but those of the
past decade or so have been rigorously tested
and have extremely long runs.
The 64-bit RNG KISS has a period greater than 10**75.
Their performances are pretty impressive.
robin
| 32-bit generators should never be used in any simulation which
| uses more than a million numbers in total. To a first approximation,
| ALWAYS use 64-bit ones.
32-bit generators have been used since time immemorial,
and most were OK for pretty well all kinds of work.
George Marsaglia's KISS RNG has a period greater than
10**35, which is far far far greater than the 10**6 figure
that you suggest.
robin
| Not all of Marsaglia's generators are free from 32-bit defects, even
| when they are used in 64-bit forms, though most of them are good or
| very good. Their actual code isn't always very portable, but that's
| easy to clean up.
His RNGs are portable to other languages. He published them
in C and Fortran, and I ported some to PL/I.
They produce the same results regardless of language.
As cross-platform RNGs, they produce the same RNs.
That's about as portable as you can get.
I agree that some of the Fortran versions need tidying up
for portability across Fortran compilers.
nm...@cam.ac.uk
robin <rob...@bigpond.com> wrote:
>
>| 32-bit generators should never be used in any simulation which
>| uses more than a million numbers in total. To a first approximation,
>| ALWAYS use 64-bit ones.
>
>32-bit generators have been used since time immemorial,
>and most were OK for pretty well all kinds of work.
Don't bet on it. I have seen a lot of people get wrong answers
by using them, cured by moving to 64-bit ones.
>George Marsaglia's KISS RNG has a period greater than
>10**35, which is far far far greater than the 10**6 figure
>that you suggest.
Sigh. I am not talking about the period, but the precision, and
the discreteness starts to become a serious problem if any one
simulation uses more than about a million numbers.
Regards,
Nick Maclaren.
frank
I thought I'd try to enumerate the reals in my implementation in the
neighborhood of .5 . I was expecting that many representations would
correspond to a single decimal one, but not so if this achieves its goal:
dan@dan-desktop:~/source$ gfortran real1.f90 -Wall -o out
dan@dan-desktop:~/source$ ./out
mult was 151
mult was 151
mult was 151
mult was 300
mult was 300
mult was 300
mult was 300
mult was 300
mult was 300
mult was 300
mult was 300
mult was 300
mult was 300
mult was 300
myarray is 0.49999991 0.49999994 0.49999997
0.50000000 0.50000006 0.50000012 0.50000018
0.50000024 0.50000030 0.50000036 0.50000042
0.50000048 0.50000054 0.50000060 0.50000066
epsilon is 3.00000000000000039E-008
number of reals in unit interval sp gfortran on ubuntu:
33333333.333333328
dan@dan-desktop:~/source$ cat real1.f90
implicit none
INTEGER, PARAMETER :: SP=SELECTED_REAL_KIND(6,37)
INTEGER, PARAMETER :: DP=SELECTED_REAL_KIND(15,307)
REAL(KIND=SP) :: begin, B, C
REAL(KIND=DP) :: V, accum, summa, epsilon
integer, parameter :: r = 15
real(KIND=SP), dimension(1:r) :: myarray
integer i, mult
v = .0000000001_dp
begin = .4999999_sp
myarray(1) = begin
do i =2,r
! inner control
mult = 0
accum = 0.0_dp
b = myarray ( i - 1)
c = myarray ( i - 1)
summa =0.0_dp
do while ( b == c )
accum = mult * v
!print *, "accum is ", accum
summa = accum + real(b, kind=dp)
c = real(summa, kind=sp)
mult = mult + 1
!print * , mult
end do
print * ,"mult was", mult
myarray(i) = c
end do
print *, "myarray is ", myarray
epsilon = mult * v
print *, "epsilon is ", epsilon
print *, "number of reals in unit interval sp gfortran on ubuntu: ", 1.0/
epsilon
end program
! gfortran real1.f90 -Wall -o out
dan@dan-desktop:~/source$
I'd be very surprised if this were methodologically sound, but it was fun
to make the try after months of not having touched a fortran compiler.
Cheers,
robin
Your claims are vague and are unsubstantiated.
You have not said anything about about Marsaglia's specific
32-bit and 64-bit generators.
Nor have you suppled any information about the 32-bit
generators that you claim to have noticed.
nm...@cam.ac.uk
robin <rob...@bigpond.com> wrote:
>|
>| >George Marsaglia's KISS RNG has a period greater than
>| >10**35, which is far far far greater than the 10**6 figure
>| >that you suggest.
>|
>| Sigh. I am not talking about the period, but the precision, and
>| the discreteness starts to become a serious problem if any one
>| simulation uses more than about a million numbers.
>
>Your claims are vague and are unsubstantiated.
>You have not said anything about about Marsaglia's specific
>32-bit and 64-bit generators.
>Nor have you suppled any information about the 32-bit
>generators that you claim to have noticed.
Do a literature search. I published the analysis.
Regards,
Nick Maclaren.
Harold Stevens
[Snip...]
> Do a literature search. I published the analysis.
I admire your cheerful optimism where "robin" is involved. :)
--
Regards, Weird (Harold Stevens) * IMPORTANT EMAIL INFO FOLLOWS *
Pardon any bogus email addresses (wookie) in place for spambots.
Really, it's (wyrd) at airmail, dotted with net. DO NOT SPAM IT.
I toss GoogleGroup (http://twovoyagers.com/improve-usenet.org/).
steve
> In article <yQuOm.56946$ze1.52...@news-server.bigpond.net.au>,
Is it this paper?
Title: Cryptographic pseudo-random numbers in simulation
Author: Maclaren, N.
Source: Cambridge Security Workshop. Fast Software Encryption.
1994. p. 185-90
It's unavailable from the Univ. of Washington library? Do you
have a pdf that is accessible?
--
steve
nm...@cam.ac.uk
steve <kar...@comcast.net> wrote:
>> >|
>> >| Sigh. =A0I am not talking about the period, but the precision, and
>> >| the discreteness starts to become a serious problem if any one
>> >| simulation uses more than about a million numbers.
>>
>> >Your claims are vague and are unsubstantiated.
>> >You have not said anything about about Marsaglia's specific
>> >32-bit and 64-bit generators.
>> >Nor have you suppled any information about the 32-bit
>> >generators that you claim to have noticed.
>>
>
>Is it this paper?
>
>Title: Cryptographic pseudo-random numbers in simulation
>Author: Maclaren, N.
>Source: Cambridge Security Workshop. Fast Software Encryption.
> 1994. p. 185-90
No. It's this one:
A Limit on the Usable Length of a Pseudorandom Number Sequence (Journal of
Statistical Computation and Simulation (1992), vol.\ 42, pp.\ 47-54).
Regrettably, I have only offprints, and I managed to bury them in
a previous office move. SOMEWHERE I will have the source of the
paper! I will take a look tonight.
Regards,
Nick Maclaren.
Ron Shepard
wrote:
> A Limit on the Usable Length of a Pseudorandom Number Sequence (Journal of
> Statistical Computation and Simulation (1992), vol.\ 42, pp.\ 47-54).
I have no doubt that the original statement in this discussion, that
less than a million numbers should be sampled from a 32-bit PRNG, is
true for many statistical simulations. However, there are indeed
many applications that are just fine with relatively short cycles.
For example, is the simulation depends not only on the current state
but on the history, then this can effectively increases the cycle
length.
Here is an example of this. Suppose that you want to generate a few
random vectors of a million real elements, and supposed your method
is successful as long as the vectors are linearly independent. In
this case, your code will work correctly with short PRNG cycles
(even those shorter than a million in length in some cases) that
repeat during this process *unless* the cycle begins to repeat for
the same vector index for two of the vectors. The odds of that
happening are one out of a million, so effectively, the history of
the sequences, which determines the index at which the cycle begins
to repeat, works to eliminate the problems with the short cycle.
This is not an obscure application of PRNGs, this is one of the more
common ones.
$.02 -Ron Shepard
Craig Powers
If you have access to the UW library, do you also have access to
Interlibrary Loan? These days, they'll deliver electronically via pdf.
frank
> A Limit on the Usable Length of a Pseudorandom Number Sequence (Journal
> of Statistical Computation and Simulation (1992), vol.\ 42, pp.\ 47-54).
>
> Regrettably, I have only offprints, and I managed to bury them in a
> previous office move. SOMEWHERE I will have the source of the paper! I
> will take a look tonight.
>
>
> Regards,
> Nick Maclaren.
I hope it turns up. It sounds pretty interesting. I haven't had good
library access since Twin Cities, where you can be Joe Anybody and they
let you export liberally from the local university libraries. They
shipped in an anthology of Gauss's works from the University of
Wisconsin. No charge, but I always paid my late fees and created them as
a matter of civic pride.
robin
| In article <168fdf4b-0398-451f...@b36g2000prf.googlegroups.com>,
| steve <kar...@comcast.net> wrote:
| >> >|
| >> >| Sigh. =A0I am not talking about the period, but the precision, and
| >> >| the discreteness starts to become a serious problem if any one
| >> >| simulation uses more than about a million numbers.
| >>
| >> >Your claims are vague and are unsubstantiated.
| >> >You have not said anything about about Marsaglia's specific
| >> >32-bit and 64-bit generators.
| >> >Nor have you suppled any information about the 32-bit
| >> >generators that you claim to have noticed.
| >>
| >> Do a literature search. =A0I published the analysis.
| >
| >Is it this paper?
| >
| >Title: Cryptographic pseudo-random numbers in simulation
| >Author: Maclaren, N.
| >Source: Cambridge Security Workshop. Fast Software Encryption.
| > 1994. p. 185-90
|
| No. It's this one:
|
| A Limit on the Usable Length of a Pseudorandom Number Sequence (Journal of
| Statistical Computation and Simulation (1992), vol.\ 42, pp.\ 47-54).
That predates ALL of George Marsaglia's RNGs published this century, right !!!
Including those to which I referred, published in this very ng.
robin
You do it. It's your responsibility.
You are the one making the vague and unsubstantiated claims.
steve
Yes, of course, I have access to the library. Unfortunately, Nick did
not
give a proper citation to his paper. The databases of scientific
papers
I searched found 2 papers by Nick; neither of which is the paper in
question. Now, that Nick gave a citation, I've found that the UW
math library has the journal.
--
steve
nm...@cam.ac.uk
frank <fr...@example.invalid> wrote:
>
>> A Limit on the Usable Length of a Pseudorandom Number Sequence (Journal
>> of Statistical Computation and Simulation (1992), vol.\ 42, pp.\ 47-54).
>>
>> Regrettably, I have only offprints, and I managed to bury them in a
>> previous office move. SOMEWHERE I will have the source of the paper! I
>> will take a look tonight.
I haven't had time - university politics :-(
>I hope it turns up. It sounds pretty interesting. I haven't had good
>library access since Twin Cities, where you can be Joe Anybody and they
>let you export liberally from the local university libraries. They
>shipped in an anthology of Gauss's works from the University of
>Wisconsin. No charge, but I always paid my late fees and created them as
>a matter of civic pride.
I can summarise the results very simply.
For cryptographic use, never use more numbers than the square root
of the period, because the (realistic) birthday test starts to fail.
For statistical use, the excess uniformity will show in the (again
realistic) chi-squared test at period^(2/3).
You need to stay well below the inverse of the precision, because
the (again realistic) spacings test starts to fail. That means
no more than about a million numbers for IEEE 32-bit, or 10^14 for
IEEE 64-bit. But it ALSO means no more than about a hundred
million for 32-bit integers converted to IEEE 64-bit results.
Because those rules are based on statistical analyses of the number
representation, and NOTHING else, they are true for all classes
of generator, past, present and future.
I didn't publish the iterated spacing test results, but they are
the ones that reject most of Marsaglia's generators in the forms
they are published. Many are fixable. I didn't invent even that
test, incidentally, and it dates from 1948.
Regards,
Nick Maclaren.
Gib Bogle
> In article <heeddb$h6f$1...@soup.linux.pwf.cam.ac.uk>, nm...@cam.ac.uk
> wrote:
>
>> A Limit on the Usable Length of a Pseudorandom Number Sequence (Journal of
>> Statistical Computation and Simulation (1992), vol.\ 42, pp.\ 47-54).
>
> I have no doubt that the original statement in this discussion, that
> less than a million numbers should be sampled from a 32-bit PRNG, is
> true for many statistical simulations. However, there are indeed
> many applications that are just fine with relatively short cycles.
> For example, is the simulation depends not only on the current state
> but on the history, then this can effectively increases the cycle
> length.
This is the case with virtually all the uses I make of RNGs. nmm1's
statement is nonsense for my applications.