Discussion
We could generate real random numbers by accessing, for example, noise on the ethernet network device but that noise might not be uniformly distributed. We typically generate pseudorandom numbers that aren't really random but look like they are. From Ross' Simulation book, we see a very easy recursive mechanism (recurrence relation) that generates values in $[0,1)$ using the modulo (remainder) operation:
$x_{i+1} = a x_i$ modulo $m$
That is recursive (or iterative and not closed form) because $x_i$ is a function of a prior value:
$x_1 = ax_0$ modulo $m$
$x_2 = ax_1$ modulo $m$
$x_3 = ax_2$ modulo $m$
$x_4 = ax_3$ modulo $m$
$...$
Ross indicates that the $x_i$ values are in [0,m-1] but setting any $x_i=0$ renders all subsequent $x_i=0$, so we should avoid that. Practically speaking, then, this method returns values in (0,1).
To get random numbers in [0,1) from $x_i$, we use $x_i / m$.
We must pick a value for $a$ and $m$ that make $x_i$ seem random. Ross suggests choosing a large prime number for $m$ that fits in our integer word size, e.g., $m = 2^{31} - 1$, and $a = 7^5 = 16807$.
Initially we set a value for $x_0$, called the random seed (it is not the first random number). Every seed leads to a different sequence of pseudorandom numbers. (In Python, you can set the seed of the standard library by using random.seed([x]).)
In [20]:
a = 16807 m = pow(2,31)-1 DFLT_SEED = 666 x_i = DFLT_SEED # this is our x_i that changes each runif01() call def runif01(): "Return a random value in U(0,1)" global x_i x_i = a * x_i % m # display(callsviz(varnames=['a','m','x_i'])) return x_i / float(m)
In [51]:
[runif01() for i in range(4)]
Out[51]:
[0.21940766983637944, 0.5847069400291457, 0.1695410698510432, 0.4767609864830789]
In [56]:
def runif(a,b): "Return a random value in U(a,b)" if b<a: # swap t = a a = b b = t return runif01()*(b-a) + a print([runif(0,10) for i in range(3)]) print([runif(5,6) for i in range(3)])
[5.931631259588352, 2.926579901448721, 7.028403648654187] [5.638012293091981, 5.072609996922599, 5.356218278108266]
In [65]:
def setseed(s): "Update the seed global variable but ensure seed > 0" global x_i if s <= 0: s = 666 x_i = s setseed(501) print([runif01() for i in range(3)]) print([runif(5,6) for i in range(3)])
[0.003921010998972231, 0.9004318597262874, 0.5582664197116468] [5.783716093648093, 5.916385943496779, 5.698552350373265]
In [70]:
import matplotlib.pyplot as plt # jupyter notebook command (ignore) %matplotlib inline sample = [runif01() for i in range(5000)] # Get 5000 random variables plt.figure(figsize=(4, 1.5)) plt.hist(sample, bins=10, density=True, alpha=0.3) plt.xlabel('Random value from U(0,1)') plt.ylabel('Probability') plt.show()