Algorithms in Python, #n+1

[Antti J Ylikoski]

I have continued my research in literature algorithms in Python.  The
algorithms in Knuth's volumes 1 -- 3 either have been incorporated
into Python, or they can be easily implemented with Python.  Quite as
John Nagle said here.  However, the Fascicles in Vol. 4 to my opinion
contain nontrivially useful material.

And, Vol. 2 has something interesting, however -- e. g. the statistical
tests on random number generators.

I have included in this entry a working Python program which carries
out the chi squared test on different random number generators.  It
demonstrates the paradigm of building Knuth's equations directly into
one's Python program, instead of the DFA/FSA approach.  Moreover,
someone over there may wish to test his/her random number generators.
Therefore, I feel that it is relevant to show this program (which is
about two pages long) in this context.

I have carried out the chi squared test on Python's own random number
generator, which to my knowledge is based on the Mersenne twister:

http://en.wikipedia.org/wiki/Mersenne_twister

and Knuth's (old) Linear Congruential method:

http://en.wikipedia.org/wiki/Linear_congruential_generator

If you run the program (with Python 3) you will see that the both
generators will pass the chi squared test.

------------------------------------------------------------------------

# Random number generation and analysis after D. E. Knuth.
#
# See The Art of Computer Programming, VOL 2, Seminumerical Algorithms,
# 3rd Edition, pp. 1 --
#
# AJY 05-16-2012.

import math, random

# First do the MIX simulator.
# AJY 05-16-2012.

# ---------------------------------------------------------------------

# Python based MIX simulator.
#
# Written in the object oriented manner.
#
# AJY 05-16-2012.
#
# Do as follows:
#
# MIX_computer = MIXSIM(memory_size, word_size)
# A, X = MIX_computer.MUL(A, X)
#
# etc etc. in the same vein.
# This constitutes just a maximally simple demo.


class MIXSIM(object):
     """ Maximally simple object oriented MIX simulator. """
     def __init__(self, mem_size, word_size):
         self.mem_size  = mem_size
         self.word_size = word_size
         self.w         = 10**word_size # Assume a base 10 computer, 
after Knuth.

     # No other methods in this class.  I said it only constitutes a demo.


# ----------------------------------------------------------------------

class LCRNG(object):
     """ The Linear Congruential Random Number Generator. AJY 
05-16-2012. """

     def __init__(self, modulus, multiplier, increment, seed):
         self.modulus    = modulus    # Modulus == m in Knuth's text.
         self.multiplier = multiplier # Multiplier == a in Knuth's text.
         self.increment  = increment  # Increment == c in Knuth.
                                      # Seed == X0 in Knuth.
         self.reg_X      = seed       # Initial value X0 to Register X.
         self.MIX        = MIXSIM(5000, 12) # Demo sample MIX computer.

     def rand(self):
         # See Knuth VOL 2 on p. 12.

         # Following calculate the next pseudo random integer:
         self.reg_A = (self.multiplier * self.reg_X) % (self.MIX.w + 1)

         # Transfer A to X for the next round.
         self.reg_X = self.reg_A

         # And now return a floating point number x such that 0 <= x < 1.
         return float(self.reg_A) / float (self.modulus)


class chi_squared_test(object):
     """ Instances are chi squared tests of random number generators. """

     def __init__(self, generator_f, buckets, iterations, name):
         self.generator_f = generator_f  # The function to be tested.
         self.buckets     = buckets      # How many discrete buckets are 
used
         self.iterations  = iterations   # The # of random numbers 
generated.
         self.name        = name         # name of this test.
         self.B = [0] * self.buckets     # initially all buckets are empty

     def run_test(self):
         """ Generate random numbers and place them in buckets. """

         for i in range(self.iterations):
             random_nr = self.generator_f() # Random number, in interval 
[0, 1)
             ind_bucket = int(self.buckets * random_nr) # Index, 0 .. 
#buckets-1
             self.B[ind_bucket] += 1

     def chi2test(self):
         """ Carry out the chi squared test of the material generated. """

          # Calculate the chi squared statistics value.  See Knuth VOL 2
          # on Page 43.

         Vsum = 0
         for i in range(self.buckets):
             Vsum += (self.B[i] - self.iterations * (1.0 / 
float(self.buckets)))**2 \
                     / (self.iterations * (1.0 / float(self.buckets)))

         V = Vsum

         print("\n Chi squared test.  Name = ", self.name)
         print("\n Chi squared value V     = ", V)
         print("\n Degrees of freedom:       ", self.buckets - 1)

         chi2_1  = (self.buckets - 1) + math.sqrt(2.0 * self.buckets) \
                   * -2.33 + (2.0/3.0)*(-2.33)**2 - (2.0/3.0)

         chi2_99 = (self.buckets - 1) + math.sqrt(2.0 * self.buckets) \
                   * 2.33 + (2.0/3.0)*(2.33)**2 - (2.0/3.0)

         print("\nValue of chi squared at the 1% level:  ", chi2_1)
         print("\nValue of chi squared at the 99% level: ", chi2_99)



# Main code.
#

if __name__ == "__main__":

     sd = 77755533311 # The seed.

     rand_x = LCRNG(10**12, 200200341, 0, sd)

     chi2_Mersenne = chi_squared_test(random.random, 1000, 1000000, 
"Mersenne")
     chi2_Linear   = chi_squared_test(rand_x.rand,   1000, 1000000, 
"Linear")

     chi2_Mersenne.run_test()
     chi2_Linear.run_test()

     chi2_Mersenne.chi2test()
     chi2_Linear.chi2test()



-----------------------------------------------------------------------------------------------




yours and V/R, Antti J Ylikoski
Helsinki, Finland, the E.U.