1 00:00:03,410 --> 00:00:09,490 Hello, bonjour, sain bainuu, salemetsiz be, 2 00:00:10,070 --> 00:00:16,830 welcome to the 10th week of Statistical Mechanics: Algorithms and Computations 3 00:00:16,830 --> 00:00:20,830 from the Physics Department of Ecole Normale Superieure. 4 00:00:21,540 --> 00:00:26,340 This week's lecture is devoted to a discussion 5 00:00:26,340 --> 00:00:30,960 of the alpha and the omega of Monte Carlo. 6 00:00:32,520 --> 00:00:36,720 The alpha will be Count Buffon's 7 00:00:36,720 --> 00:00:40,720 legendary needle-throwing experiment 8 00:00:40,720 --> 00:00:44,720 that has fascinated many people 9 00:00:44,720 --> 00:00:47,150 children and adults, 10 00:00:47,680 --> 00:00:51,800 learning circles and noble courts 11 00:00:51,800 --> 00:00:56,960 ever since its invention in 1777. 12 00:00:58,250 --> 00:01:04,570 It will lead us to a discussion of variance reduction methods. 13 00:01:06,060 --> 00:01:12,070 The omega will be Levy's stable distributions, 14 00:01:12,950 --> 00:01:17,030 that anybody interested in statistics 15 00:01:17,030 --> 00:01:21,360 and especially in Monte Carlo algorithms 16 00:01:21,360 --> 00:01:25,360 has to be at least aware of. 17 00:01:25,360 --> 00:01:29,360 It will lead us to a discussion of the case 18 00:01:29,360 --> 00:01:33,620 when the variance of an observable is infinite. 19 00:01:34,380 --> 00:01:38,420 So this week's tutorial will contain 20 00:01:38,420 --> 00:01:42,100 a short review of all the material covered 21 00:01:42,100 --> 00:01:46,100 in these ten weeks, and then we'll have a party. 22 00:01:46,100 --> 00:01:50,100 So, let's get started right away 23 00:01:50,100 --> 00:01:56,020 with this week of Statistical Mechanics: Algorithms and Computations. 24 00:02:01,890 --> 00:02:07,650 The experiment consists in randomly throwing needles 25 00:02:07,650 --> 00:02:11,650 of lengths a, onto a wooden floor 26 00:02:11,650 --> 00:02:16,390 with cracks that are a distance b apart. 27 00:02:16,910 --> 00:02:21,520 Like this.. 28 00:02:26,670 --> 00:02:31,360 The positions of the needles are random 29 00:02:31,360 --> 00:02:33,760 on an infinite floor 30 00:02:33,760 --> 00:02:37,760 the needles do not roll into cracks 31 00:02:37,760 --> 00:02:40,210 as they do in real life, 32 00:02:40,210 --> 00:02:44,210 they do not interact with each other 33 00:02:45,270 --> 00:02:48,580 and the angle phi 34 00:02:48,580 --> 00:02:53,850 is randomly distributed between 0 and 2*pi. 35 00:03:00,390 --> 00:03:03,160 So, let us introduce variables. 36 00:03:03,500 --> 00:03:09,130 x and y (clearly y is irrelevant to the problem) 37 00:03:09,830 --> 00:03:13,370 and r_center and phi. 38 00:03:14,770 --> 00:03:19,130 All the cracks in the floor are equivalent 39 00:03:20,450 --> 00:03:28,500 and there's a symmetry between the tip and the eye of the needle 40 00:03:28,500 --> 00:03:32,500 and there's also a symmetry between phi and -phi. 41 00:03:33,930 --> 00:03:37,700 So we may look at a reduced problem, 42 00:03:38,320 --> 00:03:45,130 where phi goes from 0 to pi/2, that means to 90°, 43 00:03:46,190 --> 00:03:51,870 and x_center goes from 0 to b/2. 44 00:03:53,630 --> 00:03:59,130 We suppose that the angle phi is uniformly distributed, 45 00:03:59,130 --> 00:04:03,130 and also that x_center is uniformly distributed. 46 00:04:04,620 --> 00:04:15,960 Now, the x variable of the tip is given by x_center - (a/2) cos(phi) 47 00:04:16,880 --> 00:04:24,210 and if x of the tip is smaller than 0, we have a hit. 48 00:04:24,210 --> 00:04:29,370 More precisely, the number of hits of a needle 49 00:04:29,370 --> 00:04:33,370 centered at x and oriented along phi 50 00:04:33,370 --> 00:04:37,370 is given by N_hits(x,phi) 51 00:04:37,370 --> 00:04:40,040 which is given by this formula. 52 00:04:41,360 --> 00:04:44,870 And the average that we want to compute 53 00:04:45,400 --> 00:04:51,400 is given by the integral (as we have many times in this course) 54 00:04:51,400 --> 00:04:55,400 over x from 0 to b/2 55 00:04:55,400 --> 00:04:59,480 over phi from 0 to pi/2 56 00:04:59,480 --> 00:05:02,610 N_hits(x,phi) 57 00:05:02,610 --> 00:05:07,830 divided by the normalizing integral dx from 0 to b/2 58 00:05:07,830 --> 00:05:11,830 dphi from 0 to pi/2. 59 00:05:11,830 --> 00:05:14,940 This gives the following integral. 60 00:05:15,360 --> 00:05:19,240 So, we have the solution of our problem 61 00:05:19,860 --> 00:05:23,640 under the condition that we know how to integrate 62 00:05:23,640 --> 00:05:27,640 dx from 0 to 1 of arccos(x). 63 00:05:27,640 --> 00:05:33,900 Well, think about it twice and you'll find out that it would have been wiser 64 00:05:33,900 --> 00:05:37,900 to first integrate over x 65 00:05:37,900 --> 00:05:41,900 and then integrate over phi. 66 00:05:41,900 --> 00:05:45,900 And we find that the expected number of hits 67 00:05:45,900 --> 00:05:53,420 is proportional to an integral dphi from 0 to pi/2 of cos(phi), 68 00:05:53,420 --> 00:05:57,420 and this integral we know how to do it, it is equal to 1. 69 00:05:57,420 --> 00:06:02,250 So we end up with the result that the expected number of hits 70 00:06:02,250 --> 00:06:07,140 is (a/b)(2/pi). 71 00:06:08,230 --> 00:06:12,410 So in the special case when the length of the needle 72 00:06:13,070 --> 00:06:18,450 a is equal to the distance between the floor boards, we find 73 00:06:18,450 --> 00:06:23,000 that the number of hits is equal to 2/pi. 74 00:06:23,740 --> 00:06:28,180 So, by throwing needles onto a parquet, 75 00:06:28,180 --> 00:06:31,150 we can compute the number pi. 76 00:06:31,150 --> 00:06:37,350 This is what has fascinated many people since 1777. 77 00:06:37,900 --> 00:06:42,990 We can now write a program to do the Buffon experiment ourselves. 78 00:06:44,700 --> 00:06:50,490 It suffices to sample x between 0 and b/2 79 00:06:50,490 --> 00:06:55,890 and phi between 0 and pi/2, 80 00:06:55,890 --> 00:07:04,080 and then it is necessary to check whether the tip of the needle is on the other side of the crack, as shown here. 81 00:07:05,680 --> 00:07:10,840 Well, you see that in this program we have pi 82 00:07:10,840 --> 00:07:15,270 input, because we have a random number between 0 and pi/2, 83 00:07:15,270 --> 00:07:18,580 and this program is supposed to compute pi 84 00:07:18,580 --> 00:07:23,970 so you might be interested in a historically authentic version 85 00:07:23,970 --> 00:07:29,110 that does the Buffon experiment without inputting pi. 86 00:07:29,110 --> 00:07:31,780 This is show here in the patch. 87 00:07:32,420 --> 00:07:36,690 So now we can do Buffon's experiment 88 00:07:36,690 --> 00:07:39,660 for as many needles as we want 89 00:07:40,100 --> 00:07:45,080 and you see here the output for 2000 needles 90 00:07:45,080 --> 00:07:49,080 being thrown on the parquet. 91 00:07:49,080 --> 00:07:56,890 Let's look at this program, the output, and let's ask a question: 92 00:07:56,890 --> 00:08:02,670 how do the needles hit the cracks? 93 00:08:03,870 --> 00:08:06,670 Do they hit the cracks more at the tip? 94 00:08:08,000 --> 00:08:11,680 at the eye? or at the center? 95 00:08:11,680 --> 00:08:15,680 Can you tell from this picture? 96 00:08:15,680 --> 00:08:19,070 So we want to know whether the needles 97 00:08:19,070 --> 00:08:24,510 hit the cracks more at the tip, the center or the eye. 98 00:08:24,510 --> 00:08:31,990 Mathematically speaking we can introduce an internal variable l, 99 00:08:31,990 --> 00:08:35,990 that goes from 0 to a 100 00:08:35,990 --> 00:08:39,990 and N_hits can be written as an integral 101 00:08:39,990 --> 00:08:47,640 over l from 0 to a, of N_hits(l). 102 00:08:49,890 --> 00:08:58,290 Now, we can go into a mathematical argument and a probabilistic reasoning. 103 00:08:58,920 --> 00:09:04,890 But it is much easier to look at this question experimentally. 104 00:09:05,650 --> 00:09:13,640 So now, in order to find out whether we have more hits at the eye, at the center or at the tip, 105 00:09:14,460 --> 00:09:18,540 and in order to avoid getting into probabilistic arguments 106 00:09:18,540 --> 00:09:21,220 let's take a second needle. 107 00:09:22,230 --> 00:09:28,970 Here, and let's make a little gadget. This will be gadget number 1. 108 00:09:39,880 --> 00:09:46,930 All right, you see that the center coordinate of one of the needles 109 00:09:46,930 --> 00:09:51,840 is equal to the eye coordinate of the other. 110 00:09:52,290 --> 00:09:55,840 Now see the gadget 1 111 00:09:55,840 --> 00:09:59,580 I throw one of the needles randomly onto the floor, 112 00:10:06,640 --> 00:10:13,830 but if I don't tell you which one of the needles was thrown randomly, you will not be able to tell. 113 00:10:14,910 --> 00:10:22,530 You see that the number of hits at the tip of one of the needles 114 00:10:22,530 --> 00:10:29,110 must be equal to the number of hits at the center of the other needle. 115 00:10:31,170 --> 00:10:34,450 But because these needles are equivalent, 116 00:10:34,450 --> 00:10:41,870 we have proven, more or less, that the probability 117 00:10:42,490 --> 00:10:47,080 or the mean number of hits at the center 118 00:10:47,080 --> 00:10:52,700 must be equal to the mean number of hits at the tip. 119 00:10:54,670 --> 00:11:03,540 Now, we could glue together the needles at different parts.. 120 00:11:03,540 --> 00:11:08,440 we could glue them together like this, or like this, or like this.. 121 00:11:09,260 --> 00:11:17,310 and we prove like this that the hitting probability N_hits(l) 122 00:11:17,310 --> 00:11:20,250 must be independent on l. 123 00:11:21,820 --> 00:11:29,220 And by this we can prove that N_hits is equal 124 00:11:29,220 --> 00:11:34,840 to a constant times the length of the needle. 125 00:11:35,980 --> 00:11:44,660 So the gadget that we had before was of length 3/2 and had 3/2 many hits as the original needle. 126 00:11:45,920 --> 00:11:55,040 And this we already saw in our starting calculation, because the number of hits was proportional to a, the length of the needle. 127 00:11:56,230 --> 00:12:02,660 So now, can we compute Upsilon without doing a complicated calculation 128 00:12:02,660 --> 00:12:06,660 and computing the integral of arccos(x)? 129 00:12:06,660 --> 00:12:08,360 Indeed we can. 130 00:12:08,800 --> 00:12:15,270 Because the world of gadgets is not limited to straight gadgets: 131 00:12:15,270 --> 00:12:21,500 we can glue together the two needles, like this, at an angle. 132 00:12:23,170 --> 00:12:25,500 we can do this.. 133 00:12:32,210 --> 00:12:42,630 All right, so this is my gadget number 2: it is a needle, a bent needle 134 00:12:42,630 --> 00:12:48,790 that also has an internal variable l that goes 135 00:12:49,580 --> 00:12:55,300 from 0 to the length of the needle, and I can throw this gadget. 136 00:12:57,780 --> 00:13:03,090 Now look at this bent needle gadget number 2, 137 00:13:03,090 --> 00:13:05,720 that I throw onto the floor 138 00:13:11,370 --> 00:13:17,550 Each point along the interior variable collects hits 139 00:13:17,550 --> 00:13:21,550 as if it was by itself 140 00:13:21,550 --> 00:13:25,470 So what I will find is 141 00:13:25,470 --> 00:13:25,550 that each internal point on this gadget So what I will find is 142 00:13:25,550 --> 00:13:30,360 that each internal point on this gadget 143 00:13:31,200 --> 00:13:36,440 has the same mean number of hits as it was before. 144 00:13:37,090 --> 00:13:40,130 And this remains true 145 00:13:40,130 --> 00:13:47,320 even if the needle itself was not just the joint of two straight needles, but even if it was a curved needle. 146 00:13:49,400 --> 00:13:53,690 So we find a very strong generalization 147 00:13:53,690 --> 00:14:00,200 that the mean number of hits is equal to Upsilon 148 00:14:00,200 --> 00:14:04,200 times the length of the needle, 149 00:14:04,200 --> 00:14:08,200 independently of the shape of the needle. 150 00:14:11,420 --> 00:14:15,930 So now this has a very strong consequence, 151 00:14:15,930 --> 00:14:22,220 that was first found out by Barbier in 1860. 152 00:14:24,810 --> 00:14:31,300 He considered throwing needles of length pi times b. 153 00:14:33,180 --> 00:14:35,710 So let's start with straight needles 154 00:14:35,710 --> 00:14:39,710 and let's do this experiment here on the screen 155 00:14:39,710 --> 00:14:43,710 so here I throw a first needle, a second needle, 156 00:14:43,710 --> 00:14:47,710 a third needle and so on.. 157 00:14:47,710 --> 00:14:51,710 You see that the number of hits 158 00:14:51,710 --> 00:14:55,710 can be equal to 0, 1, 2 159 00:14:55,710 --> 00:14:59,710 3, or even 4. 160 00:15:00,830 --> 00:15:04,260 Now, we can do the same experiment 161 00:15:05,320 --> 00:15:08,480 with needles of lengths pi 162 00:15:08,480 --> 00:15:14,260 but of different shape, we can take round needles 163 00:15:15,740 --> 00:15:20,540 of length pi, and let's throw them: first, 164 00:15:20,540 --> 00:15:25,630 second, third, fourth and so on.. 165 00:15:26,820 --> 00:15:33,630 So, what Barbier first noticed is that for the round needles 166 00:15:33,630 --> 00:15:40,170 we have two hits, two hits, two hits, two hits and two hits. 167 00:15:41,600 --> 00:15:45,650 But for the round needles, we also have the law 168 00:15:45,650 --> 00:15:49,650 that the mean number of hits, which is equal to two, 169 00:15:50,540 --> 00:15:55,440 is equal to Upsilon times the length of the needle. 170 00:15:56,800 --> 00:15:59,940 But this length of the needle is pi. 171 00:16:01,440 --> 00:16:07,930 So we find that Upsilon must be equal to 2/pi. 172 00:16:10,070 --> 00:16:15,800 And you see that Barbier was able to compute for us 173 00:16:15,800 --> 00:16:20,440 the outcome of the original needle experiment 174 00:16:20,440 --> 00:16:23,530 without ever doing the experiment. 175 00:16:24,050 --> 00:16:30,420 He found without any calculation that N_hits must be equal to 2/pi 176 00:16:30,420 --> 00:16:34,720 times a/b, and that is the final result. 177 00:16:37,060 --> 00:16:39,280 But now, wait a minute. 178 00:16:40,280 --> 00:16:43,290 wait a minute, we find 179 00:16:43,780 --> 00:16:48,750 that the outcome of the round needle experiment 180 00:16:48,750 --> 00:16:52,750 is the same as the outcome of the straight needle experiment? 181 00:16:54,150 --> 00:16:56,750 No, of course not. 182 00:16:56,750 --> 00:17:00,750 Let us look at the two experiments 183 00:17:00,750 --> 00:17:09,490 the straight-needle experiment and the round-needle experiment on their respective landing pads. 184 00:17:10,960 --> 00:17:17,870 So for the straight-needle experiment, as a function of x and phi, we see that we can have 185 00:17:18,150 --> 00:17:22,200 0 hits, 1 hit, 2 hits, 3 hits and 4 hits, 186 00:17:23,070 --> 00:17:27,740 and we find out that the mean number of hits is equal to 2. 187 00:17:28,700 --> 00:17:33,530 This would be quite a complicated calculation, if we did it ourselves. 188 00:17:34,220 --> 00:17:38,630 And it would be quite a complicated Monte Carlo simulation, 189 00:17:38,630 --> 00:17:42,630 because we would have a lot of fluctuations. 190 00:17:43,100 --> 00:17:46,190 Now, let's look at the landing pad 191 00:17:46,190 --> 00:17:52,300 of the round-needle experiment. Again we have the variables x and phi 192 00:17:53,130 --> 00:17:56,300 and the number of hits is 2 everywhere. 193 00:17:57,740 --> 00:18:00,300 We have no fluctuations 194 00:18:01,120 --> 00:18:06,930 and the outcome also is 2, but the two experiments are of course different. 195 00:18:06,930 --> 00:18:13,690 It is only the observable "mean number of hits" which is the same. 196 00:18:15,230 --> 00:18:19,710 So if we are only interested in computing the mean number of hits, 197 00:18:20,350 --> 00:18:26,430 we are much better off in changing our experimental set-up 198 00:18:26,430 --> 00:18:30,430 and doing the round needles, rather than the straight needles. 199 00:18:31,420 --> 00:18:34,640 What I explained here in a very simple setting 200 00:18:34,640 --> 00:18:39,810 but the very famous setting of the Buffon-Barbier needle experiment, 201 00:18:39,810 --> 00:18:43,030 has many applications 202 00:18:43,030 --> 00:18:47,030 and it is called variance reduction. 203 00:18:47,030 --> 00:18:51,030 Because you see here, in the round-needle 204 00:18:51,030 --> 00:18:57,420 landing pad, we have complete variance reduction, the variance is equal to zero, 205 00:18:58,200 --> 00:19:02,860 (of this observable) and in the other the variance is very large. 206 00:19:03,460 --> 00:19:10,620 So this variance reduction is a powerful concept, it means that we will re-do 207 00:19:10,620 --> 00:19:14,620 a new calculation and think about a new set-up 208 00:19:14,620 --> 00:19:20,260 of our algorithm for any observable that we want to compute. 209 00:19:20,890 --> 00:19:26,660 So before continuing, please take a moment to download, run and modify 210 00:19:26,660 --> 00:19:29,790 the two programs that we discussed in this section. 211 00:19:29,790 --> 00:19:33,790 There is direct_needle.py 212 00:19:34,170 --> 00:19:41,720 that does the original experiment, and direct_needle_patch.py, that used inspiration 213 00:19:41,720 --> 00:19:46,010 from the Monaco games, to avoid computing 214 00:19:46,010 --> 00:19:50,290 sines and cosines and using the input 215 00:19:50,290 --> 00:19:53,080 of pi to compute pi. 216 00:19:53,580 --> 00:19:57,360 And especially try to think over 217 00:19:57,360 --> 00:20:01,960 this amazing algorithm by Barbier, telling us 218 00:20:01,960 --> 00:20:06,240 that the mean number of hits is proportional to the length 219 00:20:06,240 --> 00:20:09,390 of the needle, independent of the shape. 220 00:20:09,390 --> 00:20:13,390 This is called Barbier's theory. 221 00:20:18,260 --> 00:20:21,930 All during this course, we have insisted 222 00:20:21,930 --> 00:20:26,810 on the relationship between sampling and integration. 223 00:20:28,120 --> 00:20:36,080 So, to rep it up, for the omega of Monte Carlo, let's do one last integral. 224 00:20:37,130 --> 00:20:44,180 We call it the gamma integral: integral from 0 to 1 of x to the gamma. 225 00:20:45,660 --> 00:20:50,540 This very simple integral can be done analytically 226 00:20:50,540 --> 00:20:54,540 the outcome is 1/(gamma+1) 227 00:20:54,990 --> 00:21:01,270 and this integral exists for gamma > -1. 228 00:21:03,500 --> 00:21:06,580 So let's do this integral by sampling. 229 00:21:08,930 --> 00:21:15,210 Simply take x as a random number between 0 and 1 230 00:21:16,340 --> 00:21:17,870 and compute 231 00:21:19,080 --> 00:21:22,190 this random number to the power of gamma. 232 00:21:23,780 --> 00:21:29,430 This is done in the program direct_gamma.py. 233 00:21:31,020 --> 00:21:37,680 And we are adding a few lines to do the error analysis. 234 00:21:37,680 --> 00:21:45,350 We compute the mean value, we compute the mean square value, and afterwards we output 235 00:21:45,350 --> 00:21:49,850 the value of the gamma integral +/- the error. 236 00:21:51,240 --> 00:21:55,740 Notice that here we have direct sampling 237 00:21:55,740 --> 00:22:00,830 algorithm, not a Markov chain sampling algorithm, so that the 238 00:22:00,830 --> 00:22:04,830 Gaussian error analysis is ok. 239 00:22:05,750 --> 00:22:11,820 So, output of this program direct_gamma.py 240 00:22:12,460 --> 00:22:16,500 is shown here for gamma=2 241 00:22:16,500 --> 00:22:24,380 1, 0, -0.2, -0.4, -0.8. 242 00:22:26,320 --> 00:22:31,200 You see that most of the times the outcome of the integral 243 00:22:31,200 --> 00:22:35,960 +/- the error agrees with the analytical result. 244 00:22:37,450 --> 00:22:42,020 But for gamma=-0.8, 245 00:22:42,750 --> 00:22:46,560 where the programs runs just fine, we get the result 246 00:22:46,560 --> 00:22:53,690 that the gamma integral is equal to 4 +/- 0.1. 247 00:22:55,440 --> 00:23:00,850 But this does not agree with the analytical result, which would be 248 00:23:00,850 --> 00:23:08,800 the gamma integral is equal to 1/(gamma+1), so it's 1/(1-0.8) 249 00:23:08,800 --> 00:23:13,290 it's 1/0.2 = 5. 250 00:23:14,640 --> 00:23:22,390 To make progress, let us modify our original program direct_gamma.py 251 00:23:23,470 --> 00:23:27,290 and compute the running averages. 252 00:23:27,840 --> 00:23:33,350 So, all during a very long simulation, we output 253 00:23:33,750 --> 00:23:36,640 what would be the current 254 00:23:36,640 --> 00:23:43,090 estimate of the gamma integral, by computing the sum of the observables 255 00:23:43,090 --> 00:23:49,140 the random number to the power of gamma, divided by the length of the run. 256 00:23:49,880 --> 00:23:54,050 So this gives direct_gamma_running.py 257 00:23:55,770 --> 00:24:01,920 Output of direct_gamma_running.py is shown here. 258 00:24:02,540 --> 00:24:11,110 We see initially we have very unpredictable results, and then for very long time 259 00:24:11,620 --> 00:24:15,110 the simulation zeros in 260 00:24:15,110 --> 00:24:21,340 on a result which is different from 5, the exact value. 261 00:24:23,120 --> 00:24:29,320 And then, after hundreds of thousands of trials, 262 00:24:29,670 --> 00:24:34,140 all of a sudden we hit an enormous value 263 00:24:34,140 --> 00:24:40,740 of x to the gamma, notice that gamma is negative, 264 00:24:40,740 --> 00:24:40,750 so if x is close to zero, of x to the gamma, notice that gamma is negative, 265 00:24:40,750 --> 00:24:45,470 so if x is close to zero, 266 00:24:45,470 --> 00:24:47,870 the value becomes very large. 267 00:24:48,780 --> 00:24:58,080 So here we see that we approach wrong result and all of a sudden we go through the roof and we continue our calculation. 268 00:24:59,070 --> 00:25:04,300 So, clearly because the running average is constant 269 00:25:04,300 --> 00:25:08,300 we also have an error estimate 270 00:25:08,300 --> 00:25:14,680 the average over dx of x^(2*gamma), which comes out much too small. 271 00:25:16,320 --> 00:25:22,610 Then the very large sample hikes up the mean value 272 00:25:22,610 --> 00:25:26,610 and also increases the error estimates. 273 00:25:26,610 --> 00:25:30,210 Now let us analyze what is the problem here. 274 00:25:30,210 --> 00:25:34,210 The problem is that we have an estimation 275 00:25:34,210 --> 00:25:37,060 of the variance of our distribution 276 00:25:37,060 --> 00:25:44,480 as the sum over the random numbers to the power of (2*gamma), as you see here in the program, 277 00:25:45,160 --> 00:25:53,390 and this sum approximates and integral, which is the integral from 0 to 1 in dx 278 00:25:53,390 --> 00:26:02,450 x^(2*gamma), so this is the integral from 0 to 1 in dx of x to the -1.6 279 00:26:02,450 --> 00:26:05,320 and this integral is infinite. 280 00:26:06,760 --> 00:26:12,740 So the distribution for gamma=-0.8 exists, 281 00:26:13,900 --> 00:26:17,360 it is finite, but it has an infinite variance. 282 00:26:19,160 --> 00:26:25,070 Clearly the situation is very difficult to analyze 283 00:26:25,070 --> 00:26:28,770 in the absence of an analytical solution. 284 00:26:30,050 --> 00:26:35,270 Now we continue in two steps to analyze the gamma integral. 285 00:26:37,090 --> 00:26:43,690 What we'll do in the first program, direct_gamma_average.py, 286 00:26:44,370 --> 00:26:49,660 is to take the average Sigma/N 287 00:26:50,460 --> 00:26:54,950 and to compute the probability distribution 288 00:26:54,950 --> 00:27:02,410 of Sigma/N as function of the length of the interval over which we average. 289 00:27:04,140 --> 00:27:10,100 We take N=1, 10, 100, 1000, 10000.. 290 00:27:11,240 --> 00:27:15,370 for N=1, we have the original integral 291 00:27:16,430 --> 00:27:23,560 and its probability distribution is given by pi(Sigma/N) 292 00:27:23,560 --> 00:27:32,960 =(Sigma/N)^(-1+1/gamma) 293 00:27:32,960 --> 00:27:37,350 as we discussed previously in lecture 4. 294 00:27:39,370 --> 00:27:47,240 So it goes like this and it decreases like 1/x^2.25, 295 00:27:47,810 --> 00:27:54,820 then we take the sum Sigma/N for N=10, we have this distribution 296 00:27:54,820 --> 00:27:59,850 then for 100 we get a new distribution, for 1000, ... 297 00:27:59,850 --> 00:28:03,850 All the averages are equal to 5 298 00:28:03,850 --> 00:28:10,010 but the most probable values only very slowly approach 5. 299 00:28:11,720 --> 00:28:16,260 And this was the reason why we obtained a wrong result 300 00:28:16,260 --> 00:28:20,820 because the average turned out to be very different from 5. 301 00:28:21,610 --> 00:28:24,820 So finally, let me give you without justification 302 00:28:25,420 --> 00:28:30,480 a way to rescale all these curves onto one universal curve. 303 00:28:31,130 --> 00:28:35,850 So this rescaling is the following, instead of considering 304 00:28:35,850 --> 00:28:39,850 the probability distribution of Sigma/N, 305 00:28:40,530 --> 00:28:45,340 we consider Sigma/N minus the mean value 306 00:28:45,340 --> 00:28:49,730 divided by N^(-1-gamma). 307 00:28:50,640 --> 00:28:58,110 In our case gamma=-0.8, and the mean value is equal to 5. 308 00:28:58,110 --> 00:29:03,330 So we have to look at the histogram of the mean value Sigma/N 309 00:29:03,330 --> 00:29:08,260 minus 5, divided by N to the -0.2. 310 00:29:09,350 --> 00:29:12,260 And you see that the same data 311 00:29:12,260 --> 00:29:21,390 that depended so much on N, fall onto the same curve approximately for large N. 312 00:29:22,220 --> 00:29:25,700 This is written up 313 00:29:25,700 --> 00:29:31,740 in the program direct_gamma_average_rescaled.py 314 00:29:31,740 --> 00:29:36,720 which also contains the analytically known 315 00:29:36,720 --> 00:29:41,520 limiting distribution in the limit N going to infinity. 316 00:29:42,750 --> 00:29:49,070 This distribution that you see approaching here for N larger and larger and that you see written here 317 00:29:49,070 --> 00:29:53,070 as analytical formula, is one of the examples 318 00:29:53,070 --> 00:29:57,070 of the famous Levy stable distributions. 319 00:29:57,990 --> 00:30:04,810 These Levy stable distributions depend on the power gamma 320 00:30:04,810 --> 00:30:08,810 of our distribution that has infinite variance 321 00:30:08,810 --> 00:30:16,440 and on the precise asymptotic behavior for x going to plus infinity and x going to minus infinity. 322 00:30:17,240 --> 00:30:21,650 So in our case the distribution for x going to plus infinity 323 00:30:21,650 --> 00:30:32,280 is 1.25 / x^2.25 and it is zero for x going to minus infinity. 324 00:30:33,560 --> 00:30:37,350 And what Levy showed in the 1930s 325 00:30:37,960 --> 00:30:45,330 is that this very erratic distributions, as we saw in our initial calculations, 326 00:30:45,870 --> 00:30:48,780 satisfy very strict 327 00:30:48,780 --> 00:30:53,860 laws and obey in the limit N going to infinity 328 00:30:53,860 --> 00:30:59,960 if they are properly rescaled, they obey these stable distributions. 329 00:30:59,960 --> 00:31:02,790 So before concluding this lecture 330 00:31:02,790 --> 00:31:09,440 and this whole lecture series, please take a moment to download, run and modify 331 00:31:09,440 --> 00:31:12,600 the programs we discussed in this lecture. 332 00:31:12,600 --> 00:31:16,600 There was direct_gamma.py 333 00:31:17,250 --> 00:31:24,490 the program that was supposed to compute the integral dx x^gamma between 0 and 1, 334 00:31:25,220 --> 00:31:34,040 and that had a lot of problems doing this for gamma between -1 and -0.5, even though the integral existed. 335 00:31:35,320 --> 00:31:40,550 Then there was direct_gamma_running.py, 336 00:31:40,770 --> 00:31:46,830 the program that computed the running average and had this point going through the roof. 337 00:31:48,080 --> 00:31:53,010 Then we had the program direct_gamma_average.py 338 00:31:53,010 --> 00:31:57,360 and direct_gamma_average_rescaled.py. 339 00:31:57,360 --> 00:32:02,950 So in conclusion we have studied in this lecture two outstanding 340 00:32:02,950 --> 00:32:06,950 problems: the original Buffon's needle 341 00:32:06,950 --> 00:32:10,950 and the famous Levy stable distributions 342 00:32:10,950 --> 00:32:16,640 that brought order and understanding to very erratic random processes. 343 00:32:17,650 --> 00:32:26,340 So let me thank you for all of your interest and very active participation during this course.