Hello, bonjour, zdravstvuyte, mambo,
welcome to the fourth tutorial of Statistical Mechanics: Algorithms and Computations
from the Physics Department of Ecole Normale Superieure.
This week we will study random numbers
and the sampling of discrete and continuous one-dimensional distributions.
You need to understand these relatively simple problems
in order to attack the more complex
sampling problem that we are interested in, for their physical applications.
For three weeks already, we heavily relied on random numbers Upsilon
generated through Upsilon=random.uniform(0,1)
This means that we called on the module 'random' of Python
and use the function 'uniform' between 0 and 1
to generate random numbers between 0 and 1, uniformly distributed with a flat distribution.
We left out many things, for example how these random numbers
are generated
on a computer, the most deterministic of all devices
and whether we can trust them.
To get an idea of what random numbers really are
let us consider a baby version of a random number generator:
it is called the congruential linear random number generator.
Starting from a seed, an initial integer value,
it generates a deterministic
but seemingly erratic sequence of integers,
that are folded back into a given interval through a modulo operation.
It is this modulo operation that gives the name congruential to this and this whole class of generators
Then, these integer numbers are divided by a large number, to give the familiar objects
the random floats between 0 and 1.
Of course, this program is deterministic
and if we run it several times
we will obtain exactly the same output.
Now, we can set the initial seed
equal for example to the system time
and then the output of our random number generator
will be a different sequence at each time that we run it.
This is also done in the random.uniform that you have in Python, unless you take special precautions.
Let us briefly discuss the properties of our naive generator, naive_ran.py
but as it always more fun to find answers by yourself
than to have them told to you
I ask you to download the program from the Coursera website
and then to answer the questions.
You may have to run the program in order to find them.
So, question 1:
the program naive_ran.py generates a sequence of integers idum
and of floats idum / float(m)
what is the range of the integers idum?
Well, the answer is clearly that the integers idum
are contained in an interval including 0 and (m-1).
This is done by the modulo operation: modulo m
Now question 2:
we have understood that the integers idum are between 0 and (m-1)
does the generator generate all the numbers between 0 and (m-1)?
or does it leave out a few of them? Are there holes?
Well, as you may find out for yourself,
for example by sorting the output of the random number generator along the first column,
you see that all the numbers between 0 and (m-1) are represented.
Question 3:
we have seen that the sequence of random integers idum is bounded
so sure enough after certain time we'll hit a number that we have seen before
what happens then? Run the program and see for yourself.
Well, this answer is very simple: the output sequence is periodic
if you hit a number a certain time, you'll generate the same sequence as you did the first time
it is this, plus the fact that there are no holes in the distribution that makes that this distribution is totally flat.
You can learn about the Python generator by reading the documentation
or by running this little program.
In our applications, we need to sample
integers and real, but also physical objects, from a general distribution pi(x)
this is of course what the Markov chain method is all about
but the direct sampling of discrete and one-dimensional continuous systems is an important subject
that Michael will next treat in this tutorial.
Before continuing, please take a moment to download, run and modify the program
naive_ran.py, if you haven't done so already
and also to familiarize yourself with basic usage of random numbers by running this program.
Throughout this course we use the Metropolis algorithm
many times, as Werner explained in the first lecture
the Metropolis acceptance rate is given by pi(a->b)
= minimum of 1 and the ratio of the statistical weights
pi(b) and pi(a) of the states b and a
But what did we actually do in order to decide whether or not a particular move is accepted or rejected?
This can be determined by comparing the acceptance rate
to a random number between 0 and 1.
If the random number is smaller than the acceptance rate the move is accepted, and otherwise it is rejected
Of course we do not even need to make this comparison
if the ratio of the statistical weights is larger than one
cause then - according to Metropolis - the acceptance rate will be 1 and the move is certainly accepted.
We have already used the Metropolis algorithm to sample non-constant probabilities.
For example in last week's homework (homework number 3)
we sampled non-constant distributions of hard-disks
but now let us use the Metropolis algorithm
to sample a simple one-dimensional function, in fact the Gaussian normal distribution
How does this look in Python? We start at x=0
then we propose a random move with a step size between (-delta) and delta
the acceptance rate of this proposed move, according to Metropolis, is shown here
and we have really understood that the minimum does not need to appear in the code.
As we see, the algorithm samples the Gaussian distribution nicely
but it much better to choose another route.
We implement a direct sampling algorithm
in fact a direct sampling rejection algorithm.
Once again, we are now talking about rejections
in a direct sampling algorithm, and not about rejections of the algorithmic transition in a Markov chain algorithm.
So how do we proceed?
Well, we can put our Gaussian into a box between (-x_cut) and x_cut
and if you want, we can do this even on the Monte Carlo beach
in a set-up which is only slightly modified from what we used in lecture 1.
Just like the children threw pebbles randomly into a circle and a square
to estimate the number pi
we can now randomly throw pebbles into our rectangular box
whenever a pebble falls below the Gaussian curve we accept it, otherwise we reject it.
So what is the probability that a pebble that is fallen at a position x
landed below the Gauss function?
Well, clearly this probability is proportional to the function itself
as there's the same number of pebbles landing at each position x
due to the constant box side for all x
then the fraction of the accepted pebbles is proportional to pi(x)
In Python, our new pebble game can be implemented as follows
we first select a random position within the rectangle
then we check whether it is above or below the Gaussian curve
and only if it is below we accept it
This direct sampling rejection algorithm again samples the Gaussian distribution very nicely
but note that we had to introduce this inelegant cutoff (-x_cut), x_cut.
Rejection sampling can be problematic
It is not due to the presence of rejections per se,
but because the rejection rate can become really enormous, even prohibitive.
To see this, let's look at a function that is not quite as friendly as a Gaussian distribution:
the distribution 2 / sqrt(x)
Maybe it troubles you that this distribution diverges for x->0
as sqrt(x)->0
and 1/sqrt(x) becomes infinite.
However, such distributions exist
and Vivien will later show you how to sample them without problems using a simple sample transformation
But for our rejection method, that works so nicely in the simple Gaussian case,
this presents a real headache, and this is because we simply don't know what box size to choose
Should we set y_max equal to 1, 10, 1000, 1000000
Our rejection rate would sky-rocket, and whatever value we choose some part of the curve will still be above our box.
Let's give it a try, we implement the rejection method (just like we did in the Gaussian case) and set y_max=100
although the result does not look so bad, in order to get a feeling for the potential danger of cutting off
the divergent probability distribution at finite y_max, we compare with a Markov chain algorithm.
Note that we do not need to introduce any cut-off as the Markov chain algorithm does not have any problem coping with the infinite probability density.
The programs records the highest values of y and the lowest values of x that it obtains
while it randomly moves along the x-axis between 0 and 1.
This is the output of our program, in a few second of runtime
it goes down to values in x as small as 10^(-7) and up to values in y of a few thousand, way beyond our chosen cutoff of 100.
Should we actually use this Markov chain algorithm?
Alberto in a few moments and Vivien later will show you that this is not at all called for.
So before continuing, please familiarize yourself thoroughly with the way the Metropolis acceptance probability is implemented.
Furthermore take a moment to download, run and modify the programs we just discussed. On the Coursera website you will find markov_gauss.py
that implements the Metropolis choice directly.
There's also the graphic version, markov_gauss_movie.py.
Furthermore you will find the program reject_direct_gauss_cut.py
note that there's the word 'cut' in its name
as we have to cut-off the tails of the distribution at (-x_cut), x_cut.
In this example this was not a problem, as the distribution is suppressed
exponentially for large absolute values of x.
There's also the program reject_direct_inv_sqrt.py
Now we had to introduce a cut-off in y, and it was not exponentially small.
This is really dangerous and we must avoid it at any price.
Before Alberto will next introduce the method of tower sampling
I invite you to take a look at markov_inv_sqrt_movie.py
and take it apart, just in case you still doubt the existence of divergent distributions.
Get it to run for other values of the exponent: positive ones
pi proportional to x, x^2, x^3, for x between 0 and 1
but also for negative ones: pi proportional to x^(-1/2)
x^(-1) and x^(-3/2), again for x between 0 and 1
There are many strange phenomena to observe, and we'll get back to them in week 10, the final week of this course.
Sampling non-uniform finite distributions
has an archetypical example in what we call the Saturday night problem.
As discussed by Werner in lecture 4,
we imagine to have k possible evening activities
and we do not feel equally enthusiastic about them:
go out with friend, visit family, do your homework,
sports, laundry, and so on..
We know all the pi_k probabilities
but we can still have troubles in deciding what to do.
This means that we have troubles in sampling the probability pi_0, pi_1, pi_2 and so on
we can adapt Michael's rejection method, and instead of sampling x between -x_cut and x_cut,
we can draw integer random numbers
0, 1, .., K-1, select an activity
which we'll be accepted if the pebble falls inside the red box.
This occurs with a probability which is proportional to the weight of the activity.
As Michael explained, this algorithm can have problems because of the high rejection rate.
But for the Saturday night problem this is ok, we can program the algorithm and plan our activities week by week.
Bad luck..
Let's look again to the last figure.
You can see that the number of pebbles you should draw before having a single hit
is typically equal to the ratio between pi_max and pi_average.
Where pi_max is the maximal weight and pi_average is the mean value between all the activities.
This number can be very large.
The tower sampling method is much faster and much more elegant.
Instead of placing the boxed one next to the other, we pile up them, one on top of each other,
and the pebble is drawn from the bottom to the top of the tower.
Once again, you choose your activity with a probability which is proportional to the weight of the activity
But the algorithm is much faster, and the output is statistically the same.
in Python, this gives the following rejection-free program
How we can actually find the box j where the pebble falls?
for a small number of activities K
we can just go through the ordered table pi_0, pi_1, ..
until we find a j for which pi_(j-1)<Upsilon<pi_j
For larger K, it is better to use the bisection method
We should first check if Upsilon is smaller or larger than pi_(K/2)
and then cut the range of possible indices j in two for any subsequent iteration.
After log_2(k) operations, we will find the box, and this number gives the complexity of the algorithm.
The tower sampling method is quite amazing because it implements a rejection-free algorithm for sampling, but it is not optimal.
The best method to sample discrete distributions is the Walker method.
Let's consider the colored boxes in the figure.
They represent the discrete distribution that we want to sample.
Let us arrange them in a special way.
We first compute the box height average, as shown in the figure
there are bigger boxes and smaller boxes than average
Take one bigger box and a smaller box,
put on the top of the smaller box the bigger box,
cut it off, return what is left on its original position.
You can continue and repeat this operation for other bigger boxes and other smaller boxes.
Like here, and here.
At the end, you will have a boxing scheme, where for each position you have at most two boxes
which add up to the average. And now sampling is very easy.
You notice that the logarithmic complexity has improved.
We achieved a sampling method that takes a constant work per pebble.
for an arbitrary number of evening activities.
And this ingenious algorithm can be
programmed in Python in a few lines.
Sampling a discrete distribution is so easy that you can think that to sample a continuous one-dimensional distribution you should cut it up in small boxes
and use the tower sampling or the Walker algorithm.
This is not so.
In fact, next in this tutorial, Vivien will show directly the continuous limit of the tower sampling,
And come up with very nice and amazing algorithms.
But before listening to him, please take a moment
to download, run and modify the programs we discussed in this section.
Here you find the tower sampling method,
with the logarithmic search that locates the boxes.
Here we find also the ingenious algorithm by Walker, that produces the nice two-slate parquet.
In the final part of this tutorial, let us now take the continuous limit of the discrete distributions
This corresponds to cases where instead of having a finite or enumerable series of probabilities,
pi_0, pi_1, pi_2...
you have a distribution pi(x) which is a function of the continuous variable x, taking values inside the real line.
To sample such distributions, let's come back to Alberto's tower sampling trick
and remain with a discretized version of the function for a little while, as described here.
To apply Alberto's method we take the different boxes with colors indicated here (red, green, blue, yellow, pink)
and place them on top of each other.
With respect to the Saturday night problem, here
each box corresponds to a position x,
rather than to a specific activity.
To keep track of these x positions, the best way is the one which is displayed here.
Again, if we place a pebble randomly between the bottom and the top of the tower (that we can take equal to one)
we'll fall in the red box with a probability which is proportional to the height of this box.
But this is exactly what we want.
Using this scheme, one can trace from the position of the pebble, to the box and to the position x that we want to sample
according to the distribution (which is discrete) pi(x). This is what is displayed here.
Now let's make the discretization finer and finer.
As you observe, there is an envelope
which develops around the tower construction
and the value of the envelope at a given position x
is in fact simply the integral of the distribution pi(x) up to this point.
In fact, in this example we have taken pi(x) to be a Gaussian normal distribution
so this means that we have a second algorithm, that we call gauss_transform.py, to sample this distribution.
It is a bit complicated, because we first have to compute the integral of this function
Sometimes we can compute analytically the definite integral between minus infinity and x
For the Gaussian distribution - however - Werner's trick
does not work. It only works between minus infinity and plus infinity.
Fortunately for us, the definite integral which is called the error function and which is displayed here
has been studied by mathematicians for centuries.
Although it cannot be reduced to a simpler formula,
we can evaluate it numerically with high efficiency
but now wait a minute, we don't need this function!
So, generically, let's call phi(x) the integral of the distribution pi up to the position x.
We don't need the value phi(x) associated to a value of x, but rather the contrary
we need the position x corresponding to phi, which is the random number between 0 and 1.
All in all, we will need the inverse function of phi(x)
Again, for the Gaussian distribution, the inverse error function cannot be reduced to a simple form
but it can be evaluated numerically as displayed in this program
where its value can be computed by a simple call to the erfinv function
The output of the program is shown here.
To summarize in a mathematical way the mathematical way the discussion that we have had up to now on a specific example,
let us make explicit the correspondence between the discrete and continuous cases.
As displayed here, the first step consists in computing the definite integral of pi between minus infinity and x.
This is the function that we call phi(x).
The second step, which is illustrated there, in the discrete cases
corresponds to determining the position at which the pebble has fallen into.
In the continuous case, it consists in determining the value x which corresponds to the random number Upsilon
distributed between 0 and 1.
Now, let us remark that the procedure that we have described up to now
is in fact the sample transformation described by Werner in this week's lecture, lecture 4.
In any case, the two hurdles that we have at hand
computing the integral and solving the equation phi(x)=Upsilon
can be overcome in general for one-dimensional distribution with various degrees of freedom.
Sometimes the integral can be computed analytically,
sometimes one has to resort to numerical integration
and to interpolations.
Let us also remark that if pi(x) presents divergences or singularities,
however its integral phi(x) up to point x is in fact a much more regular function.
Indeed, since pi is positive and can be normalized,
phi is growing from 0 to 1 in a monotonous manner.
This means - in the end- that the process we have described to sample pi(x) is very stable in general
AS a second example we now turn our attention to the case of a distribution pi(x) equal to a constant
times x to the power of gamma.
where x is a real number between 0 and 1, 0 being excluded.
A few minutes ago, Michael has studied this distribution in the case where gamma = -1/2.
This was a distribution which blew up when x was going to zero.
In fact, we have already considered this distribution for gamma=2
in the program direct_spheres_3d.py.
Let us now apply the sample transformation to sample this distribution.
This is what we use in the following program
with the proper normalization, the function phi(x) can be computed, and it is equal to x^(gamma+1).
The sampling of the distribution is thus given by the following relation
which is x^(1+gamma)=Upsilon
which means that x=Upsilon^(1/(1+gamma))
where as usual Upsilon is a random number between 0 and 1 with a flat distribution.
So, in the end it means that to sample pi you simply have to take this random number
and to take it to the power of 1/(1+gamma).
Isn't this remarkably easy?
It took us a bit of brain work to think about it, but in the end the results justify our effort.
Now it is your turn to download, run and modify the programs we have presented in this section.
Have a look at gauss_transform.py and make sure that you have understood how the sample transformation works.
Further, try out gamma_transform.py
and compare directly the results that you obtain
to the two other algorithms: the Markov chain
algorithm and the direct sampling algorithm that Michael presented earlier.
In conclusion, we have discussed in this tutorial random numbers and probability distribution
and how to reshape them using rejection methods or transformation methods.
We found great rejection-free sampling programs for discrete and continuous one-dimensional distributions,
and solved the direct sampling problem in these cases.
Markov chain methods could generally be avoided.
The situation is much more complicated in high-dimensional situations
for example the hard-spheres simulations that we did last week or the quantum problem that we will next turn our attention to.
In the meantime have fun with Homework session 4.