Re: Issues with rand()

[Daniel J S. <dan...@ie...>]

Thanks Péter.  The rand() function should probably be revisited with your comments in mind.  More below.

Juhász Péter wrote:
> Dear gnuplot developers,
> 
> a recent thread[1] in the newsgroup made me play with the rand()
> function, and I've found some issues with it:
> 
> 1) rand() silently accepts string arguments:
> gnuplot> print rand("foo")
> 0.222457440974512
> 
> While this is not particularly bad in itself, but it's not consistent
> with the behavior of other numerical functions (which reject string
> arguments), and it's undocumented.

Actually, my observation is that functions having real arguments will reject string inputs while functions with optional complex arguments will not reject string inputs.  For example

gnuplot> print gamma("test")

gnuplot> print gamma("test")
               ^
         undefined value

gnuplot> print erfc("test")
1.0

Perhaps the gnuplot code needs to be altered to also reject strings if the argument is possibly complex, unless a string is a valid input argument.  (But I don't know of any, off hand.)


> Its effect is the same as that of
> rand(0) - which have misled the user in the thread[1], because he
> thought that rand("time") sets the random seed to the current time.

Are you proposing there should be an input argument "time"?  The command line documentation says nothing about using the current time as a seed.


> 2) It is possible to lock the PRNG into a state where the returned
> numbers are not random at all:
> gnuplot> print rand(0.5)
> 0.999999999450207
> gnuplot> print rand(0)
> 0.999999999450207
> gnuplot> print rand(0)
> 0.999999999450207
> gnuplot> print rand(0)
> 0.999999999450207
> gnuplot> print rand(0)
> 0.999999999450207
> 
> I'll refrain from linking to the relevant xkcd or Dilbert comics.

Well, here is where gnuplot is deficient in my opinion.  The documentation doesn't indicate what the input value should be, a float or an integer.  Here's the code to illustrate why this is important:

    /* Construct new seed values from input parameter */
    /* FIXME: Ideally we should allow all 64 bits of seed to be set */
    if (real(init) > 0.0) {
	if (real(init) >= (double)(017777777777UL))
	    int_error(NO_CARET,"Illegal seed value");
	if (imag(init) >= (double)(017777777777UL))
	    int_error(NO_CARET,"Illegal seed value");
	seed1 = (int)real(init);
	seed2 = (int)imag(init);
	if (seed2 == 0)
	    seed2 = seed1;
    }
    FPRINTF((stderr,"ranf: seed = %lo %lo        %ld %ld\n", seed1,seed2));

Comparing floating point with an integer cast to a double is slightly dodgy:

real(init) >= (double)(017777777777UL)

Hence the uneasy FIXME statement in the comments.  More than that, casting floating point quantities to integer quantities for the seed is problematic.  This means that any input 0 < x < 1 will put the random number generator in the "zero" state, i.e., stuck.  For example, try

plot rand(0.5)
plot rand(0.34)
plot rand(0.789)

The thing is, reading the documentation, there is no indication that 0 < x < 1 is invalid.  In fact, it's not illogical for the user to assume the value should be in (0,1) since that is what the rand function outputs.

There are a couple ways to fix this.  (Both methods require giving more detailed description of what numeric type the seed can be.)  It seems to be the input of the function can't be either float or integer, it has to be one or the other.  If the input is a float, then perhaps 0 < x < 1 and x should be cast to its equivalence in terms of seed.  (But it probably isn't a unique equivalence, and therein lies a problem with random number generation with this method when the core algorithm is integer based.  Even the existing method of casting loses resolution and has the same problem.)  If the input is an integer, then the values should be translated directly to the seeds of the random number generator.  (My fear is that the gnuplot code, designed to be generic, has no way of preserving full integer resolution because everything is cast to a float from the command line.)


> 3) Consider the following plot:
> 
> set xrange [0:2**22]
> set samples 1000
> plot '+' u 1:(rand($1)) w l
> 
> The rand() function, when called with a nonzero argument, sets the seeds
> based on the argument and returns the first pseudo-random number from
> the sequence associated those seeds. As the plot shows, there is a
> rather obvious dependence between rand's argument and the returned
> number, in fact the dependence is linear if only the lower 21-or-so bits
> of the argument are considered.

This may be true, but I don't know if this is a measure of randomness.


> This is a problem if we want to use the standard trick of initializing
> the PRNG with the current time (as the user wanted to in [1]):
> gnuplot> print rand(real(system("date +%s")))
> 0.596925441003784
> gnuplot> print rand(real(system("date +%s")))
> 0.596924809567054
> gnuplot> print rand(real(system("date +%s")))
> 0.596924178130323
> gnuplot> print rand(real(system("date +%s")))
> 0.596923546693592
> 
> date +%s returns the time in the Unix time_t format (32 bit integer),
> however, the upper bits rarely change in that format.

I'm not sure this is that important.  This is just initializing the random number generator.  So long as the function's output changes for different inputs, it is effectively starting the random number generator in a different state and the very first step should appear to have randomness.  If you are doing an experiment in which trials should have "random"ly independent first values, simply drop the first outcome of the random number generator after the seed (or only seed once, not for each new trial).  (Probably good advice no matter what method is used to seed the random number generator...The random number generator is designed to look independent across outcomes, but starting each trial with a seed constructed from exterior isn't necessarily going to have good random properties.)


> 1) and 2) may be bugs in the implementation that are easy to fix, but 3)
> is a deeper problem. Of course, this kind of behavior can be expected
> from a linear congruence generator - but then maybe it's time to
> consider changing to a different algorithm.

To me, what's problematic is the loss of resolution because of the way gnuplot is programmed.  In other words, we need to go back and address this comment:

    /* FIXME: Ideally we should allow all 64 bits of seed to be set */

which likely means allowing integer inputs.  Do we need some kind of numeric representation which is either a float or an integer, whichever allows the input value to be stored without rounding of some form?  Also, the "rand" documentation should be changed to indicate exactly what numeric type the seed should be.

Dan