How to Solve it by Computer #7

Algorithm 3.6: Use the linear congruential method to generate a uniform set of pseudo-random numbers.

This is great and something I never did before – like most of the stuff here. The method were using is the linear congruential method which should create a uniform distribution of numbers.

The basic idea is to use the following formula to create the next pseudo-random number:

for n >= 0

However, there are specific criteria for each parameter.

Parameter :
Parameter m: m >= length sequence required, (ax+b) mod m is an integer
Parameter a: if m is a power of 2 then a mod 8 = 5; if m is a power of 10 then a mod 200 = 21
sqrt(m) < a < m – sqrt(m)
(a-1) should be a multiple of every prime dividing into m
if m is a multiple of 4 then (a-1) should be a multiple of 4
Parameter b: b should be odd, not a multiple of 5

These are a lot of requirements and the easiest thing to do is use existing values. Wikipedia provides common parameters.

The ones provided by the book are:

m = 4096
b = 853
a = 109

Here’s the code:

(defn gen-pseudo-random
  " Generates a sequence of random integers
    less than 4096."
  [seed]
  {:pre [((complement neg?) seed), (< seed 4096)]}
    (let [a 109
          b 853
          m 4096
          random-number (mod (+ (* a seed) b) m)]
  (lazy-seq
      (cons random-number (gen-pseudo-random random-number)))))

Super straight forward. Let’s see how good it works. We want to compare the average.

Problem 3.6.1: Confirm that the algorithm repeats after generating m random numbers. Compute the mean value and variance for the set of m pseudo-numbers.

user=> (last (take 4096 (gen-pseudo-random 40)))
40
user=> (last (take 4096 (gen-pseudo-random 2345)))
2345

Seems to work.

user=> (def pseudo-set (take 4096 (gen-pseudo-random 222)))
#'user/pseudo-set
user=> (average pseudo-set)
4095/2
user=> (variance pseudo-set)
1398101.25

Looks fine. Let’s see how big the variance over the average over the set is.

 

user=> (dotimes [n 20] (println "Seed:" n "Average:" (average (take 4096 (gen-pseudo-random n)))))
Seed: 0 Average: 4095/2
Seed: 1 Average: 4095/2
Seed: 2 Average: 4095/2
Seed: 3 Average: 4095/2
Seed: 4 Average: 4095/2
Seed: 5 Average: 4095/2
Seed: 6 Average: 4095/2
Seed: 7 Average: 4095/2
Seed: 8 Average: 4095/2
Seed: 9 Average: 4095/2
Seed: 10 Average: 4095/2
Seed: 11 Average: 4095/2
Seed: 12 Average: 4095/2
Seed: 13 Average: 4095/2
Seed: 14 Average: 4095/2
Seed: 15 Average: 4095/2
Seed: 16 Average: 4095/2
Seed: 17 Average: 4095/2
Seed: 18 Average: 4095/2
Seed: 19 Average: 4095/2

Yup. By the way, here are the functions for calculating the average and variance:

 

(defn average
  "Calculates the average of a collection"
  [coll]
  (/ (reduce + coll) (count coll)))


(defn variance
  "Calculates the variance of a collection"
  [coll]
  (* (/ 1 (count coll))
     (let [my-average (average coll)]
       (reduce +
               (map #(Math/pow (- % my-average) 2) coll)))))

Pretty standard.