[R] Needed: Understading runif() output :-)

[Ivan Frohne]

At 06:38 25-05-00, Kjetil Kjernsmo wrote:
>Dear all,
>
>I have been trying to understand what runif() is telling me.
>I am generating lots of numbers (billions and billions (wow, I've dreamed
>about saying that for many years... :-) ), for a distribution that has the
>following quantile function:
>     1 / (2 * sqrt(1 - p))
>(that is, the distribution has a lower cutoff)
>As you can imagine, this has rather heavy upper tail. I was looking at the
>largest values, and it looked as if the largest values appeared again and
>again. Now, it wasn't in itself that large values were strange, since I'm
>generating so many numbers, but that the largest were very much larger
>than the second largest numbers, and that exactly the same number appeared
>again and again. First I thought it was a bug, and I'm sorry to have
>wasted r-devels time with a bug report.
>
>I started running the same simulation with different RNGs and they all
>seem to generate numbers in "quantized states". Then, I started to look
>into what runif() gives, and let it print 13 digits.
>In the below output, I use the "Mersenne-Twister" RNG and I have generated
>1e+10 numbers (100000 at a time) and I print a line if it the number is
>above 10000 (my dist, the left coloumn), the right coloumn are runif()
>the corresponding outputs.
>[1] 3.276800000000e+04 9.999999997672e-01
>[1] 13377.479981919865     0.999999998603
>[1] 1.158523750296e+04 9.999999981374e-01
>[1] 1.036215143684e+04 9.999999976717e-01
>[1] 1.158523750296e+04 9.999999981374e-01
>[1] 1.036215143684e+04 9.999999976717e-01
>[1] 1.036215143684e+04 9.999999976717e-01
>[1] 13377.479981919865     0.999999998603
>
>So, it seems that the runif() outputs are "quantized" too. The question
>is: What is the reason for this?


Mersenne Twister pseudo-random numbers are based on 32-bit
unsigned integers.  R uses 64-bit doubles, and if you want all
bits random, you need to combine two uniforms.  Something like

x <- (1.0 - 2^-32) * first.uniform + 2 ^ -32 * second.uniform

should work.  For efficiency, use the actual values of the two
constants, and make first.uniform and second.uniform vectors of
1000 or more uniforms.

--Ivan Frohne

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._