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.
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.