The contents of this Jupyter Notebook lecture notes are:
As usual I recommend you follow along by typing the code snippets into your own file. Don't forget to call the packages etc. at the start of each code file.
In [2]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
Examples of physical processes that are/can be modelled as random include:
Chaotic systems - again not truely random in the sense of radioactive decay, but can be modelled as random.
Human or animal behaviour can also be modelled as random in some circumstances.
There are many different ways to generate uniform random numbers over a specified range (such as 0-1). Physically, we can for example:
spin a roulette wheel
draw balls from a lottery
throw darts at a board
thow dice
However, when we wish to use the numbers in a computer, we need a way to generate the numbers algorithmically.
Numerically/arithmetically - use a sequential method where each new number is a deterministic function of the previous numbers.
But: this destroys their true randomness and makes them at best, "pseudo-random".
However, in most cases, it is sufficient if the numbers “look” uniformly distributed and have no correlation between them. i.e. they pass statistical tests and obey the central limit theorem.
For example consider the function:
$x' = (ax + c) \mod m$
where $a$, $c$ and $m$ are integer constants, and $x$ is an integer variable. Recall that "$n \mod m$" means you calculate the remainder when $n$ is divided by $m$.
Now we can use this to generate a sequence of numbers by putting the outcome of this equation ($x'$) back in as the new starting value ($x$). These will act like random numbers. Try it.....
Starting from $x = 1$ write a short programme which generates 100 values in this sequence and plots them on a graph. Please use the following inputs:
a = 1664525
c = 1013904223
m = 4294967296
Tip 1: python syntax for "mod m" is:
%m
So your base code will look like:
xp = (a*x+c)%m
Extension problem: this won't work for all values of a, c and m. Can you find some which don't generate pseudo-random numbers?
This is an example of a simple pseudo-random number generator (PRNG). Technically it's a "linear congruential random number generator". Things to note:
For many codes this is sufficient, but you can do better. Fortunately python (Numpy) comes with a number of better versions as in built packages, so we can benefit from the expertise of others in our computational physics codes.
All pseudo-random number generators (PRNG) should possess a few key properties. Namely, they should
be fast and not memory intensive
be able to reproduce a given stream of random numbers (for debugging/verification of computer programs or so we can use identical numbers to compare different systems)
be able to produce several different independent “streams” of random numbers
have a long periodicity, so that they do not wrap around and produce the same numbers again within a reasonably long window.
To obtain a sequence of pseudo-random numbers:
But note:
The length of the sequence of non-repeating numbers is the period of the PRNG.
It is relatively easy to build PRNGs with periods long enough for many practical applications, but one must be cautious in applying PRNG's to problems that require very large quantities of random numbers.
Almost all languages and simulation packages have good built-in generators. In Python, we can use the NumPy random library, which is based on the Mersenne-Twister algorithm developed in 1997.
In [3]:
#Review the documentation for NumPy's random module:
np.random?
Some basic functions to point out (we'll get to others in a bit):
random() - Uniformly distributed floats over [0, 1]. Will include zero, but not one. If you inclue a number, n in the bracket you get n random floats. randint(n,m) - A single random integer from n to m-1
In [21]:
#print 5 uniformly distributed numbers between 0 and 1
print(np.random.random(5))
#print another 5 - should be different
print(np.random.random(5))
In [4]:
#print 5 uniformly distributed integers between 1 and 10
print(np.random.randint(1,11,5))
#print another 5 - should be different
print(np.random.randint(1,11,5))
Notice you have to use 1-11 for the range. Why?
In [22]:
#If you want to save a random number for future use:
z=np.random.random()
print("The number is ",z)
#Rerun random
print(np.random.random())
print("The number is still",z)
Write a programme that generates and prints out two random numbers between 1 and 6. This simulates the rolling of two dice.
Now modify the programme to simulate making 2 million rolls of two dice. What fraction of the time do you get double six?
Extension: Plot a histogram of the frequency of the total of the two dice over the 2 million rolls.
In [3]:
np.random.seed(42)
for i in range(4):
print(np.random.random())
In [4]:
np.random.seed(42)
for i in range(4):
print(np.random.random())
In [6]:
np.random.seed(39)
for i in range(4):
print(np.random.random())
You might want to do this for:
In some circumstances you will want to write code which simulates various events, each of which happen with a probability, $p$.
This can be coded with random numbers. You generate a random number between zero and 1, and allow the event to occur if that number is greater than $p$.
For example, consider a biased coin, which returns a head 20% of the time:
In [10]:
for i in range(10):
if np.random.random()<0.2:
print("Heads")
else:
print("Tails")
Simulate the decay of 1000 thallium atoms over time, using random number generators to mimick the random process of atomic decay.
Thallium-208 decays to stable lead (208) with a half life of 3.053 minutes.
The standard equation of radioactive decay (for the number of atoms in the sample as a function of time) is:
$$N(t) = N(0) 2^{-t/\tau}$$where tau is the half life, N(0) is the number of atoms at t=0. Notice that both t and tau must be in the same units.
The fraction of atoms which have not yet decayed at any time t, is then:
$$\frac{N(t)}{N(0)} = 2^{-t/\tau}$$So then the probability that any given atom has decayed by time t (which is the same as the fraction of atoms that have decayed by that time) is:
$$p(t) = 1 - 2^{-t/\tau}$$Use time steps of 1 second and make a plot of the number of thallium and lead atoms as a function of time until 20 minutes have passed. Overplot the half-life of thallium.
The radioactive decay problem we did in class is an example of a problem which can be coded more efficiently if we could draw random numbers from a distribution other than uniform.
In radioactive decay, the probability of decay in a time step $dt$ is
$$dp = 1 - 2^{-dt/\tau}$$This can be expressed as
$$dp = 1 - \exp(-\frac{dt}{\tau} \ln2)$$and the second term can be expanded using it's Taylor expansion, to give a first order approximation of
$$dp = \frac{\ln2}{\tau} dt$$If we want to know the probability of decay between time t and dt this is then:
$$P(t) dt = 2^{-t/\tau} \frac{\ln2}{\tau} dt$$which is a non-uniform probability (it's larger for small $t$).
We could more quickly calculate the decay of a sample of N atoms by drawing random numbers from the above distribution to give decay times for each atom.
There are at least two methods to generate non-uniform random numbers from a function like the Numpy random function which generates uniform random numbers:
In the rejection method, you generate more random numbers than are needed, and reject them if they are above the value of the probability function you wish to sample from.
This is no faster than the probability method - it's basically equivalent.
The transformation method is faster. In the transformation method, you need to find a function $f(z)$ which converts unfirmly distribtued random numbers $z$ into random numbers with the desired non-uniform sampling.
Most of the Numpy packages for non-uniform random numbers make use of the transformation method. For most physical processes you can make use of a pre-written package which generates random numbers with the desired shape, but in the interests of education, we will work through a couple of examples.
Suppose you have a uniformly distributed random number, $z$, drawn from 0 to 1 (and zero elsewhere). You want a number, $x$ with a probability distribution $p(x)$, and there is going to be some function which relates $x$ and $z$.
It must be the case that
$$p(x) dx = q(z) dz $$(where $q(z)$ is 1 in the interval zero to 1, and zero elsewhere).
If we integrate this equation we get
$$ \int_{-\infty}^{x(z)} p(x) dx = \int_0^z dz $$so then
$$ z = \int_{-\infty}^{x(z)} p(x) dx $$will tell us the function needed to transform the uniform random numbers to random numbers from the distribution p(x).
The exponential probability distribution is
$$p(x) = \mu e^{-\mu x}$$where the $\mu$ is a normalization factor (so that the probability integrates to one from zero to infinity).
This is the radioactive decay probability in the case that $\mu = \ln2/\tau$.
Put this into the transformation equation above we get
$$ z = \int_{-\infty}^{x(z)} \mu e^{-\mu x} dx = 1 - e^{-\mu x} $$which can be used to generate numbers from the exponential distribution.
Redo the radioactive decay problem, but now make use of non-uniform random numbers, to draw N random numbers representing the decay time of N atoms.
You can make the plot much more quickly, with a bit of help from the sort function in numpy.
Extension: Also try the rejection method. Notice how similar the code is to the first way we did the problem.
Simulate the Brownian motion of a particle in two dimension, and make an animation of the output.
Set up the particle to be confined to a square grid of size $L \times L$ (where $L$ is an odd number), and represent it's position using a two-dimensional array of integers.
The particles starts in the middle of the grid. At each time step it moves randomly one lattice point in any direction. This is called a 'random walk'. The particle will "bounce" off the walls of the grid (ie. at the borders it cannot move in the direction which would take it outside of the grid).
Make an animation of the path of the particle over 1 million time steps.
In [ ]: