Pseudo-random number generators¶

The classic and infamous RANDU. This is a simplified version of a version courtesy of ptomato on github.

"...its very name RANDU is enough to bring dismay into the eyes and stomachs of many computer scientists!"

Donald E. Knuth, The Art of Computer Programming

In [1]:
import numpy as np

class RANDU:
"Pseudorandom number generator RANDU, a (flawed) linear congruential PRNG."
def __init__(self, seed=None):
self.seed(seed)

def seed(self, s):
self._state = hash(s) % 2147483648  # == 2^31

def _random(self):
self._state = (65539 * self._state) % 2147483648
return self._state / 2147483648

def random(self, size=None):
"Return one random number or an array of them, with size elements."
if size is None:
return self._random()
else:
r = self._random
return np.array([r() for i in range(size)])

In [2]:
RANDU().random(10)

Out[2]:
array([ 0.68546926,  0.9696463 ,  0.64865447,  0.16511013,  0.15277057,
0.43063228,  0.2088585 ,  0.37746052,  0.38503659,  0.91307486])
In [3]:
%matplotlib notebook
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns

In [4]:
plt.rcParams['figure.figsize'] = (10, 6)
plt.style.use('seaborn')
sns.set_context('talk', font_scale=1.4)

In [5]:
npts = 1_000

pts = RANDU().random(npts*3).reshape(npts, 3)
fig, ax = plt.subplots(subplot_kw=dict(projection='3d'))
ax.scatter(pts[:, 0], pts[:, 1], pts[:, 2], s=5);
plt.title("RANDU");


This is what similar output looks like for Numpy, which uses the Mersenne Twister PRNG, whose period is $2^{19937}-1 (\approx \cal{O}(10^{6001}))$ and which does not suffer from known serial correlations:

In [6]:
pts = np.random.uniform(size=(npts, 3))
fig, ax = plt.subplots(subplot_kw=dict(projection='3d'))
ax.scatter(pts[:, 0], pts[:, 1], pts[:, 2], s=5)
plt.title("Numpy Random");


A terrible PRNG¶

Let's implement a simple linear congruential PRNG of the type

$$X_{k+1} = a X_k \mod M$$

with particularly bad constants $a$, $M$ and seed $X_0$:

In [7]:
a = 3
M = 64
seed = 17

In [8]:
state = seed
numbers = []
for i in range(100):
state = (a * state) % M
numbers.append(state)
numbers = np.array(numbers)
print(numbers)

[51 25 11 33 35 41 59 49 19 57 43  1  3  9 27 17 51 25 11 33 35 41 59 49 19
57 43  1  3  9 27 17 51 25 11 33 35 41 59 49 19 57 43  1  3  9 27 17 51 25
11 33 35 41 59 49 19 57 43  1  3  9 27 17 51 25 11 33 35 41 59 49 19 57 43
1  3  9 27 17 51 25 11 33 35 41 59 49 19 57 43  1  3  9 27 17 51 25 11 33]

In [9]:
sorted(set(numbers))

Out[9]:
[1, 3, 9, 11, 17, 19, 25, 27, 33, 35, 41, 43, 49, 51, 57, 59]
In [10]:
len(_)

Out[10]:
16
In [11]:
np.nonzero(numbers==1)

Out[11]:
(array([11, 27, 43, 59, 75, 91]),)
In [12]:
np.diff(_)

Out[12]:
array([[16, 16, 16, 16, 16]])
In [14]:
%matplotlib inline

In [15]:
randx = numbers[:-1]
randy = numbers[1:]
plt.scatter(randx, randy, facecolors='none', edgecolors='b', linewidth=2);