## Random Number Test

### July 7, 2017

First, we need a random number generator to test. Here is a simple Park-Miller random number generator:

`(define seed 20170707)`
`(define (lcg) (set! seed (modulo (* 16807 seed) 2147483647)) seed)`

We also need some auxiliary functions; `zero-pad` uses `digits` from the Standard Prelude:

```(define (zero-pad n len)
(let* ((ds (digits n 2)))
(append (make-list (- len (length ds)) 0) ds)))```
```(define (make-bits n)
(let loop ((n n) (bits (list)))
(if (zero? n) bits
(loop (- n 1) (append (zero-pad (lcg) 31) bits)))))```
`(define (one? n) (= n 1))`

Our first test counts the number of 1-bits in the output; for n 31-bit random numbers, there should be 31 n / 2 1-bits. We look at a range of two standard deviations around the mean:

```(define (count-ones n)
(let* ((bits (make-bits n))
(expected (* n 31 1/2))
(std-dev (sqrt (/ n 4)))
(lo (- expected (* 2 std-dev)))
(hi (+ expected (* 2 std-dev))))
(display "Number of ones: ")
(display (length (filter one? bits)))
(newline)
(display "Expected: ")
(display (inexact->exact (round lo)))
(display " to ")
(display (inexact->exact (round hi)))
(newline)))```
```> (count-ones 100000)
Number of ones: 1549923
Expected: 1549684 to 1550316```

It took me a few tries to get a test in range; most of them were just slightly outside the expected range. Our second test counts runs:

```(define (max-run n)
(let* ((bits (make-bits n))
(expected (- (/ (log (* 0.5 31 n)) (log 0.5))))
(std-dev (/ (log 2)))
(lo (- expected (* 2 std-dev)))
(hi (+ expected (* 2 std-dev))))
(let loop ((bits bits) (current-run 0) (max-run 0))
(if (null? bits)
(begin
(display "Maximum run length: ")
(display max-run)
(newline)
(display "Expected: ")
(display (inexact->exact (round lo)))
(display " to ")
(display (inexact->exact (round hi)))
(newline))
(if (zero? (car bits))
(loop (cdr bits) 0 (max current-run max-run))
(loop (cdr bits) (+ current-run 1) max-run))))))```
```> (max-run 100000)
Maximum run length: 21
Expected: 18 to 23```

Most of the tests were within range, so it was a better result than the count of 1-bits. As Park and Miller explained, their random number generator is a minimum standard, so we don’t expect a perfect result. You can run the program at http://ideone.com/wW5NJO.

By the way, RANDU, the random number generator that Knuth described as “truly horrible,” doesn’t fare badly on our tests:

`(define (randu) (set! seed (modulo (* 65539 seed) 2147483648)) seed)`
```> (count-ones 100000)
Number of ones: 1549930
Expected: 1549684 to 1550316```
```> (max-run 100000)
Maximum run length: 22
Expected: 18 to 23```

That means our test suite is incomplete, not that RANDU is a good random number generator.

Pages: 1 2

### One Response to “Random Number Test”

1. Paul said

I posted on ideone test code in Python. The RANDU generator seems to fail the tests. Other generators like Park-Miller, Java, Python and Ansi C seem to be better.