Hello, bonjour, sain bainuu, salemetsiz be, welcome to the 10th week of Statistical Mechanics: Algorithms and Computations from the Physics Department of Ecole Normale Superieure. This week's lecture is devoted to a discussion of the alpha and the omega of Monte Carlo. The alpha will be Count Buffon's legendary needle-throwing experiment that has fascinated many people children and adults, learning circles and noble courts ever since its invention in 1777. It will lead us to a discussion of variance reduction methods. The omega will be Levy's stable distributions, that anybody interested in statistics and especially in Monte Carlo algorithms has to be at least aware of. It will lead us to a discussion of the case when the variance of an observable is infinite. So this week's tutorial will contain a short review of all the material covered in these ten weeks, and then we'll have a party. So, let's get started right away with this week of Statistical Mechanics: Algorithms and Computations. The experiment consists in randomly throwing needles of lengths a, onto a wooden floor with cracks that are a distance b apart. Like this.. The positions of the needles are random on an infinite floor the needles do not roll into cracks as they do in real life, they do not interact with each other and the angle phi is randomly distributed between 0 and 2*pi. So, let us introduce variables. x and y (clearly y is irrelevant to the problem) and r_center and phi. All the cracks in the floor are equivalent and there's a symmetry between the tip and the eye of the needle and there's also a symmetry between phi and -phi. So we may look at a reduced problem, where phi goes from 0 to pi/2, that means to 90°, and x_center goes from 0 to b/2. We suppose that the angle phi is uniformly distributed, and also that x_center is uniformly distributed. Now, the x variable of the tip is given by x_center - (a/2) cos(phi) and if x of the tip is smaller than 0, we have a hit. More precisely, the number of hits of a needle centered at x and oriented along phi is given by N_hits(x,phi) which is given by this formula. And the average that we want to compute is given by the integral (as we have many times in this course) over x from 0 to b/2 over phi from 0 to pi/2 N_hits(x,phi) divided by the normalizing integral dx from 0 to b/2 dphi from 0 to pi/2. This gives the following integral. So, we have the solution of our problem under the condition that we know how to integrate dx from 0 to 1 of arccos(x). Well, think about it twice and you'll find out that it would have been wiser to first integrate over x and then integrate over phi. And we find that the expected number of hits is proportional to an integral dphi from 0 to pi/2 of cos(phi), and this integral we know how to do it, it is equal to 1. So we end up with the result that the expected number of hits is (a/b)(2/pi). So in the special case when the length of the needle a is equal to the distance between the floor boards, we find that the number of hits is equal to 2/pi. So, by throwing needles onto a parquet, we can compute the number pi. This is what has fascinated many people since 1777. We can now write a program to do the Buffon experiment ourselves. It suffices to sample x between 0 and b/2 and phi between 0 and pi/2, and then it is necessary to check whether the tip of the needle is on the other side of the crack, as shown here. Well, you see that in this program we have pi input, because we have a random number between 0 and pi/2, and this program is supposed to compute pi so you might be interested in a historically authentic version that does the Buffon experiment without inputting pi. This is show here in the patch. So now we can do Buffon's experiment for as many needles as we want and you see here the output for 2000 needles being thrown on the parquet. Let's look at this program, the output, and let's ask a question: how do the needles hit the cracks? Do they hit the cracks more at the tip? at the eye? or at the center? Can you tell from this picture? So we want to know whether the needles hit the cracks more at the tip, the center or the eye. Mathematically speaking we can introduce an internal variable l, that goes from 0 to a and N_hits can be written as an integral over l from 0 to a, of N_hits(l). Now, we can go into a mathematical argument and a probabilistic reasoning. But it is much easier to look at this question experimentally. So now, in order to find out whether we have more hits at the eye, at the center or at the tip, and in order to avoid getting into probabilistic arguments let's take a second needle. Here, and let's make a little gadget. This will be gadget number 1. All right, you see that the center coordinate of one of the needles is equal to the eye coordinate of the other. Now see the gadget 1 I throw one of the needles randomly onto the floor, but if I don't tell you which one of the needles was thrown randomly, you will not be able to tell. You see that the number of hits at the tip of one of the needles must be equal to the number of hits at the center of the other needle. But because these needles are equivalent, we have proven, more or less, that the probability or the mean number of hits at the center must be equal to the mean number of hits at the tip. Now, we could glue together the needles at different parts.. we could glue them together like this, or like this, or like this.. and we prove like this that the hitting probability N_hits(l) must be independent on l. And by this we can prove that N_hits is equal to a constant times the length of the needle. So the gadget that we had before was of length 3/2 and had 3/2 many hits as the original needle. And this we already saw in our starting calculation, because the number of hits was proportional to a, the length of the needle. So now, can we compute Upsilon without doing a complicated calculation and computing the integral of arccos(x)? Indeed we can. Because the world of gadgets is not limited to straight gadgets: we can glue together the two needles, like this, at an angle. we can do this.. All right, so this is my gadget number 2: it is a needle, a bent needle that also has an internal variable l that goes from 0 to the length of the needle, and I can throw this gadget. Now look at this bent needle gadget number 2, that I throw onto the floor Each point along the interior variable collects hits as if it was by itself So what I will find is that each internal point on this gadget So what I will find is that each internal point on this gadget has the same mean number of hits as it was before. And this remains true even if the needle itself was not just the joint of two straight needles, but even if it was a curved needle. So we find a very strong generalization that the mean number of hits is equal to Upsilon times the length of the needle, independently of the shape of the needle. So now this has a very strong consequence, that was first found out by Barbier in 1860. He considered throwing needles of length pi times b. So let's start with straight needles and let's do this experiment here on the screen so here I throw a first needle, a second needle, a third needle and so on.. You see that the number of hits can be equal to 0, 1, 2 3, or even 4. Now, we can do the same experiment with needles of lengths pi but of different shape, we can take round needles of length pi, and let's throw them: first, second, third, fourth and so on.. So, what Barbier first noticed is that for the round needles we have two hits, two hits, two hits, two hits and two hits. But for the round needles, we also have the law that the mean number of hits, which is equal to two, is equal to Upsilon times the length of the needle. But this length of the needle is pi. So we find that Upsilon must be equal to 2/pi. And you see that Barbier was able to compute for us the outcome of the original needle experiment without ever doing the experiment. He found without any calculation that N_hits must be equal to 2/pi times a/b, and that is the final result. But now, wait a minute. wait a minute, we find that the outcome of the round needle experiment is the same as the outcome of the straight needle experiment? No, of course not. Let us look at the two experiments the straight-needle experiment and the round-needle experiment on their respective landing pads. So for the straight-needle experiment, as a function of x and phi, we see that we can have 0 hits, 1 hit, 2 hits, 3 hits and 4 hits, and we find out that the mean number of hits is equal to 2. This would be quite a complicated calculation, if we did it ourselves. And it would be quite a complicated Monte Carlo simulation, because we would have a lot of fluctuations. Now, let's look at the landing pad of the round-needle experiment. Again we have the variables x and phi and the number of hits is 2 everywhere. We have no fluctuations and the outcome also is 2, but the two experiments are of course different. It is only the observable "mean number of hits" which is the same. So if we are only interested in computing the mean number of hits, we are much better off in changing our experimental set-up and doing the round needles, rather than the straight needles. What I explained here in a very simple setting but the very famous setting of the Buffon-Barbier needle experiment, has many applications and it is called variance reduction. Because you see here, in the round-needle landing pad, we have complete variance reduction, the variance is equal to zero, (of this observable) and in the other the variance is very large. So this variance reduction is a powerful concept, it means that we will re-do a new calculation and think about a new set-up of our algorithm for any observable that we want to compute. So before continuing, please take a moment to download, run and modify the two programs that we discussed in this section. There is direct_needle.py that does the original experiment, and direct_needle_patch.py, that used inspiration from the Monaco games, to avoid computing sines and cosines and using the input of pi to compute pi. And especially try to think over this amazing algorithm by Barbier, telling us that the mean number of hits is proportional to the length of the needle, independent of the shape. This is called Barbier's theory. All during this course, we have insisted on the relationship between sampling and integration. So, to rep it up, for the omega of Monte Carlo, let's do one last integral. We call it the gamma integral: integral from 0 to 1 of x to the gamma. This very simple integral can be done analytically the outcome is 1/(gamma+1) and this integral exists for gamma > -1. So let's do this integral by sampling. Simply take x as a random number between 0 and 1 and compute this random number to the power of gamma. This is done in the program direct_gamma.py. And we are adding a few lines to do the error analysis. We compute the mean value, we compute the mean square value, and afterwards we output the value of the gamma integral +/- the error. Notice that here we have direct sampling algorithm, not a Markov chain sampling algorithm, so that the Gaussian error analysis is ok. So, output of this program direct_gamma.py is shown here for gamma=2 1, 0, -0.2, -0.4, -0.8. You see that most of the times the outcome of the integral +/- the error agrees with the analytical result. But for gamma=-0.8, where the programs runs just fine, we get the result that the gamma integral is equal to 4 +/- 0.1. But this does not agree with the analytical result, which would be the gamma integral is equal to 1/(gamma+1), so it's 1/(1-0.8) it's 1/0.2 = 5. To make progress, let us modify our original program direct_gamma.py and compute the running averages. So, all during a very long simulation, we output what would be the current estimate of the gamma integral, by computing the sum of the observables the random number to the power of gamma, divided by the length of the run. So this gives direct_gamma_running.py Output of direct_gamma_running.py is shown here. We see initially we have very unpredictable results, and then for very long time the simulation zeros in on a result which is different from 5, the exact value. And then, after hundreds of thousands of trials, all of a sudden we hit an enormous value of x to the gamma, notice that gamma is negative, so if x is close to zero, of x to the gamma, notice that gamma is negative, so if x is close to zero, the value becomes very large. So here we see that we approach wrong result and all of a sudden we go through the roof and we continue our calculation. So, clearly because the running average is constant we also have an error estimate the average over dx of x^(2*gamma), which comes out much too small. Then the very large sample hikes up the mean value and also increases the error estimates. Now let us analyze what is the problem here. The problem is that we have an estimation of the variance of our distribution as the sum over the random numbers to the power of (2*gamma), as you see here in the program, and this sum approximates and integral, which is the integral from 0 to 1 in dx x^(2*gamma), so this is the integral from 0 to 1 in dx of x to the -1.6 and this integral is infinite. So the distribution for gamma=-0.8 exists, it is finite, but it has an infinite variance. Clearly the situation is very difficult to analyze in the absence of an analytical solution. Now we continue in two steps to analyze the gamma integral. What we'll do in the first program, direct_gamma_average.py, is to take the average Sigma/N and to compute the probability distribution of Sigma/N as function of the length of the interval over which we average. We take N=1, 10, 100, 1000, 10000.. for N=1, we have the original integral and its probability distribution is given by pi(Sigma/N) =(Sigma/N)^(-1+1/gamma) as we discussed previously in lecture 4. So it goes like this and it decreases like 1/x^2.25, then we take the sum Sigma/N for N=10, we have this distribution then for 100 we get a new distribution, for 1000, ... All the averages are equal to 5 but the most probable values only very slowly approach 5. And this was the reason why we obtained a wrong result because the average turned out to be very different from 5. So finally, let me give you without justification a way to rescale all these curves onto one universal curve. So this rescaling is the following, instead of considering the probability distribution of Sigma/N, we consider Sigma/N minus the mean value divided by N^(-1-gamma). In our case gamma=-0.8, and the mean value is equal to 5. So we have to look at the histogram of the mean value Sigma/N minus 5, divided by N to the -0.2. And you see that the same data that depended so much on N, fall onto the same curve approximately for large N. This is written up in the program direct_gamma_average_rescaled.py which also contains the analytically known limiting distribution in the limit N going to infinity. This distribution that you see approaching here for N larger and larger and that you see written here as analytical formula, is one of the examples of the famous Levy stable distributions. These Levy stable distributions depend on the power gamma of our distribution that has infinite variance and on the precise asymptotic behavior for x going to plus infinity and x going to minus infinity. So in our case the distribution for x going to plus infinity is 1.25 / x^2.25 and it is zero for x going to minus infinity. And what Levy showed in the 1930s is that this very erratic distributions, as we saw in our initial calculations, satisfy very strict laws and obey in the limit N going to infinity if they are properly rescaled, they obey these stable distributions. So before concluding this lecture and this whole lecture series, please take a moment to download, run and modify the programs we discussed in this lecture. There was direct_gamma.py the program that was supposed to compute the integral dx x^gamma between 0 and 1, and that had a lot of problems doing this for gamma between -1 and -0.5, even though the integral existed. Then there was direct_gamma_running.py, the program that computed the running average and had this point going through the roof. Then we had the program direct_gamma_average.py and direct_gamma_average_rescaled.py. So in conclusion we have studied in this lecture two outstanding problems: the original Buffon's needle and the famous Levy stable distributions that brought order and understanding to very erratic random processes. So let me thank you for all of your interest and very active participation during this course.
