<jupyter_start><jupyter_text>Pseudo-random number generators**<NAME>** [<EMAIL>](mailto:<EMAIL>)  Properties of PRNGs+ dimension of output    - commonly 32 bits, but some have more    + number of states    - dimension of state space in bits    - sometimes state = output, but better generators generally have output = f(state)+ period    - maximum over initial states of the number of states visited before repeating    - period &le; number of states    - if state has $s$ bits, period $\le 2^s$    - for some PRNGs, period is much less than number of states    - for some seeds for some PRNGs, number of states visited is much less than period    + $k$-distribution    - suppose $\{X_i\}$ is sequence of $P$ $w$-bit integers    - define $t_v(X_i)$ to be the first $v$ bits of $X_i$    - $\{X_i\}$ is $k$-distributed to $v$-bit accuracy if each of the $2^{kv}-1$ possible nonzero $kv$-bit vectors occurs equally often among the $P$ $kv$-bit vectors$$ (t_v(X_i),\,t_v(X_{i+1}), \ldots ,t_v(X_{i+k-1}))\quad (0\le i<P),$$and the zero vector occurs once less often.    - amounts to a form of uniformity in $k$-dimensional space, over an entire cycle    - does not measure dependence or other "serial" properties+ sensitivity to initial state; burn-in    - many PRNGs don't do well if the seed has too many zeros    - some require many iterations before output behaves well    - for some seeds, some PRNGs repeat very quickly  Some PRNGs  Middle SquareDates to Franciscan friar ca. 1240 (per [Wikipedia](https://en.wikipedia.org/wiki/Middle-square_method)); reinvented by von Neumann ca. 1949.Take $n$-digit number, square it, use middle $n$ digits as the "random" and the new seed.E.g., for $n=4$, take $X_0 = 1234$.$1234^2 = 1522756$, so $X_1 = 2275$.$2275^2 = 5175625$, so $X_2 = 7562$.+ $10^n$ possible states, but not all attainable from a given seed+ period at most $8^n$, but can be very short. E.g., for $n=4$,    - 0000, 0100, 2500, 3792, & 7600 repeat forever    - 0540 &rightarrow; 2916 &rightarrow; 5030 &rightarrow; 3009 &rightarrow; 0540  Linear Congruential Generators (LCGs) $$ X_{n+1} = (aX_n +c)\mod m.$$LCG period is at most $m$. Hull-Dobell Theorem: the period of an LCG is $m$ for all seeds $X_0$ iff+ $m$ and $c$ are relatively prime+ $a-1$  is divisible by all prime factors of $m$+ $a-1$  is divisible by 4 if $m$ is divisible by 4  Marsaglia (PNAS, 1968): Random Numbers Fall Mainly in the PlanesMultiplicative congruential generators ($c=0$), Lehmer (1949).![Marsaglia68](./Figs/marsaglia68.png)  RANDURANDU is a particularly bad [linear congruential generator](https://en.wikipedia.org/wiki/Linear_congruential_generator) promulgated in the 1960s and widely copied.RANDU is given by the recursion$$ X_{j+1} = 65539 X_j\mod 2^{31}.$$Period is ($2^{29}$); all outputs are odd integers. Triples of values from RANDU fall on 15 planes in 3-dimensional space, as shown below.<jupyter_code>from __future__ import division
%matplotlib inline
import math
import numpy as np
import scipy as sp
from scipy.misc import comb, factorial
from scipy.optimize import brentq
from scipy.stats import chisquare, norm
import scipy.integrate
from random import Random
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# LCG; defaults to RANDU, a particularly bad choice
class lcgRandom: # defaults to RANDU: BEWARE!
    def __init__(self, seed=1234567890, A=0, B=65539, M = 2**31): 
        self.state = seed
        self.A = A
        self.B = B
        self.M = M
        
    def getState(self):
        return self.state, self.A, self.B, self.M
    
    def setState(self,seed=1234567890, A=0, B=65539, M = 2**31):
        self.state = seed
        self.A = A
        self.B = B
        self.M = M

    def nextRandom(self):
        self.state = (self.A + self.B * self.state) % self.M
        return self.state/self.M

    def random(self, size=None):  # vector of rands
        if size==None:
            return self.nextRandom()
        else: 
            return np.reshape(np.array([self.nextRandom() for i in np.arange(np.prod(size))]), size)
    
    def randint(self, low=0, high=None, size=None):  # integer between low (inclusive) and high (exclusive)
        if high==None:  # numpy.random.randint()-like behavior
            high, low = low, 0
        if size==None:
            return low + np.floor(self.nextRandom()*(high-low)) # NOT AN ACCURATE ALGORITHM! See below.
        else:
            return low + np.floor(self.random(size=size)*(high-low))
        
# generate triples using RANDU
reps = int(10**5)
randu = lcgRandom(12345)
xs = np.transpose(randu.random(size=(reps,3)))
# plot the triples as points in R^3
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(xs[0],xs[1], xs[2])

plt.rcParams['figure.figsize'] = (18.0, 18.0) 

ax.view_init(-100,110)

plt.show()<jupyter_output><empty_output><jupyter_text>Wichmann-Hill (1982)Sum of 3 LCGs. Period is 6,953,607,871,644.     def WH(s1, s2, s3):        s1 = (171 * s1) % 30269        s2 = (172 * s2) % 30307        s3 = (170 * s3) % 30323        r = (s1/30269 + s2/30307 + s3/30323) %  1        return [r, s1, s2, s3] The right way, the wrong way, and the Microsoft way.WH generally not considered adequate for statistics, but was (nominally) the PRNG in Excel for severalgenerations. Excel did not allow the seed to be set, so analyses were not reproducible.  ![mcCullough](./Figs/notWichmannHill08.png)<NAME>., 2008. Microsoft Excel's 'Not The Wichmann–Hill' random number generators_Computational Statistics & Data Analysis_, _52_, 4587&ndash;4593doi:10.1016/j.csda.2008.03.006  Mersenne Twister (MT)  Matsumoto & Nishimura (1997)+ example of "twisted generalized feedback shift register"+ period $2^{19937}-1$, a Mersenne Prime+ $k$-distributed to 32-bit accuracy for all $k \in \{1, \ldots, 623\}$. + passes DIEHARD and most of TestU01 (see below)+ standard in many packages:    - GNU Octave, Maple, MATLAB, Mathematica, Python, R, Stata    - Apache, CMU Common Lisp, Embeddable Common Lisp, Free Pascal, GLib, PHP, GAUSS, IDL, Julia, Ruby, SageMath, Steel Bank Common Lisp, Scilab, Stata, GNU Scientific Library, GNU Multiple Precision Arithmetic Library, Microsoft Visual C++.    - SPSS and SAS offer MT, as does C++ (v11 and up)+ generally considered adequate for statistics (but not for cryptography); however, will trouble that in this work, esp. for "big data"+ usual implementation has 624-dimensional state space, but TinyMT uses only 127 bits+ seeding complicated, since state is an array+ can take a while to "burn in," especially for seeds with many zeros+ output for close seed states can be close+ 2002 update improves seeding+ completely predictable from 624 successive outputs+ problems discovered in 2007 (see TestU01, below)<jupyter_code># Python implementation of MT19937 from Wikipedia 
# https://en.wikipedia.org/wiki/Mersenne_Twister#Python_implementation

def _int32(x):
    # Get the 32 least significant bits.
    return int(0xFFFFFFFF & x)

class MT19937:

    def __init__(self, seed):
        # Initialize the index to 0
        self.index = 624
        self.mt = [0] * 624
        self.mt[0] = seed  # Initialize the initial state to the seed
        for i in range(1, 624):
            self.mt[i] = _int32(
                1812433253 * (self.mt[i - 1] ^ self.mt[i - 1] >> 30) + i)

    def extract_number(self):
        if self.index >= 624:
            self.twist()

        y = self.mt[self.index]

        # Right shift by 11 bits
        y = y ^ y >> 11
        # Shift y left by 7 and take the bitwise and of 2636928640
        y = y ^ y << 7 & 2636928640
        # Shift y left by 15 and take the bitwise and of y and 4022730752
        y = y ^ y << 15 & 4022730752
        # Right shift by 18 bits
        y = y ^ y >> 18

        self.index = self.index + 1

        return _int32(y)

    def twist(self):
        for i in range(624):
            # Get the most significant bit and add it to the less significant
            # bits of the next number
            y = _int32((self.mt[i] & 0x80000000) +
                       (self.mt[(i + 1) % 624] & 0x7fffffff))
            self.mt[i] = self.mt[(i + 397) % 624] ^ y >> 1

            if y % 2 != 0:
                self.mt[i] = self.mt[i] ^ 0x9908b0df
        self.index = 0<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup