This page was generated from `source/notebooks/L2/3_randomnumbers.ipynb`_.

Random numbers

Random numbers are widely used in science and engineering computations. They can be used to simulate noisy data, or to model physical phenomena like the distribution of velocities of molecules in a gas, or to act like the roll of dice in a game. Monte Carlo simulation techniques, which will be part of a later theory lecture rely heavily and random number. Processes like single photon emission or Brownian motion are stochastic processes, with an intrisic randomness as well. But there are even methods for numerically evaluating multi-dimensional integrals using random numbers.

The basic idea of a random number generator is that it should be able to produce a sequence of numbers that are distributed according to some predetermined distribution function. NumPy provides a number of such random number generators in its library numpy.random.

The random number functions can be imported by

from numpy.random import *

Uniformly distributed random numbers

The rand(num) function creates an array of num floats uniformly distributed on the interval from 0 to 1.

import matplotlib.pyplot as plt
import numpy as np
from numpy.random import *
%config InlineBackend.figure_format = 'retina'
plt.rcParams.update({'font.size': 12,
                     'axes.titlesize': 20,
                     'axes.labelsize': 18,
                     'axes.labelpad': 14,
                     'lines.linewidth': 1,
                     'lines.markersize': 10,
                     'xtick.labelsize' : 16,
                     'ytick.labelsize' : 16,
                     'xtick.top' : True,
                     'xtick.direction' : 'in',
                     'ytick.right' : True,
                     'ytick.direction' : 'in',})

rand()

0.1930417336574951

If you supply the argument num to the rand(num) function you obtain an array of equally distributed numbers with num elements.

rand(5)

array([0.62537446, 0.10574321, 0.72463188, 0.34811094, 0.77941787])

You may also obtain a multi-dimensional array if you give two or more numbers to the the rand function.

rand(5,2)

array([[0.9323968 , 0.20114768],
       [0.20210667, 0.26948309],
       [0.7311368 , 0.09990769],
       [0.64669995, 0.7223557 ],
       [0.7401539 , 0.10669448]])

b=np.linspace(0,1,50)
plt.hist(rand(10000),bins=b,density=1)
plt.show()

Normally distributed random numbers

The function randn(num) produces a normal or Gaussian distribution of num random numbers with a mean of \(\mu=0\) and a standard deviation of \(\sigma=1\). They are distributed according to

\begin{equation} p(x)dx=\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(x-\mu)^2}{2\sigma^2}}dx \end{equation}

Similarly as all the other random number function, you may supply one or multiple arguments to the rand() function. The result is again a multi-dimensional array of random numbers.

b=np.linspace(-5,5,100)
plt.hist(randn(10000),bins=b,density=1)
plt.show()

If you want to create a normal distribution with different standard deviation and mean value, you have to multiply by the new standard deviation and add the mean.

sigma=10
mu=10

d=sigma * np.random.randn(1000) + mu
d.mean()

10.121706338361593

Note: Physics Interlude

Brownian Motion, Random Walk

The following lines create random numbers, which are distributed by a normal distribution. You can use such normally distributed random numbers to generate a random walk. Such a random walk is a simple representation of Brownian motion of colloids suspended in a liquid.

Suppose a colloidal particle in solution is kicked by its surrounding solvent molecules such that it does a small step \((\Delta x_i,\Delta y_i)\) in a random direction. The length of the step \((\Delta x_i,\Delta y_i)\) will be normally distributed. If this motion is in 2 dimensions, then the position in x and y direction after N steps is

\begin{equation} x(N)=\sum_{i=1}^{N}\Delta x_{i} \end{equation}

\begin{equation} y(N)=\sum_{i=1}^{N}\Delta y_{i} \end{equation}

The position is therefore just the cummulative sum of a sequence of normally distributed random numbers, which is easy to realize in Python.

x,y=[randn(10000).cumsum(),randn(10000).cumsum()]

Note that the above example uses the cumsum() function. The cumsum function of an array of number \([x_0,x_1,x_2,..,x_n]\) delivers an array with a progressive sum of elements \([x_0,x_0+x_1,x_0+x_1+x_2,...,x_0+...+x_n]\). The line

x,y=[randn(1000).cumsum(),randn(1000).cumsum()]

is therefore all you need to generate a random walk in two dimensions.

plt.figure(figsize=(6,6))
plt.plot(x,y)
plt.xlabel('x') # set the x-axis label
plt.ylabel('y') # set the y-axis label
plt.show()

Exponentially distributed numbers

A number of processes in physics reveal an exponential statistics. For example the probability to find a molecule under gravity at a certain height \(h\) is distributed by a Boltzmann law for a non-zero temperature.

\begin{equation} p(h)dh = p_{0}\exp(-m\cdot g\cdot h/k_{\rm B}T) dh \end{equation}

On the other hand, the probability to emit a photon spontaneously after a certain time \(t\) follwings the excited state preparation (two level system) is also exponentially distributed.

\begin{equation} p(t)dt=p_{0}\exp(-t/\tau) \end{equation}

Thus after each excitation i.e. by a laser pulse, a molecule emits a single photon with the probability \(p(t)\) and the whole exponential character in the statistics is only appearing after repeating the experiment several times.

The exponential distribution of numpy can be supplied with two numbers.

exponential(b, n)

The parameter b is giving the decay parameter. The number n is optional and giving the number of samples to be provided. The numbers are distributed according to

\begin{equation} \frac{1}{b}\exp(-x/b) \end{equation}

exponential(1)

0.664906606125074

exponential(1,10000)

array([0.89196938, 0.92497306, 0.12805347, ..., 4.39709775, 1.56870637,
       1.35026235])

You may want to test the changes in the exponential distribution with the parameter b.

b=np.linspace(0,10,50)
plt.hist(exponential(2,10000),bins=b,density=1);

Random distribution of integers

The function randint(low, high, num) produces a uniform random distribution of num integers between low (inclusive) and high (exclusive).

randint(1,20,10)

array([ 4,  1, 19, 12,  3, 16, 14, 19,  5,  2])

items = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
shuffle(items)

print(items)

[8, 7, 6, 9, 10, 4, 2, 3, 5, 1]

There are a number of other methods available in the random module of numpy. Please refere to the documentation.