Buffon needle pricing experiment approaches pi but doesn't get perfect value

Question

I work on a very typical homework at university. I have to write a C program to evaluate pi through the Buffon Needle using the Monte Carlo method. I think my program works correctly, but I never got pi rights. It is always close to typical 3.14, but sometimes its 3.148910 is sometimes 3.13894. What is the reason for this? Can I do anything to make it more “equal” for exact value? (comments of codes in Polish, comments in () are in English.)

#include <stdio.h> #include <math.h> #include <time.h> #include <stdlib.h> int main(int argc, char **argv) { srand(time(NULL)); int i; //wprowadzenie licznika (counter) double l; // dlugosc igly (lenght of needle) double n;// liczba prob (number of trials) double L; // szerokosc miedzy liniami (width between lines) double x; // kat miedzy igla a normalna do linii (angle) double pi; // szacowana wartosc liczby pi (pi) double P=0; //ilość prób spełniających warunek styku igły z linią (ammount of trials that worked) double d; //odleglosc srodka igly od linii (distance between centre of needle and line) double y; //minimalizacja obliczeń double stosun; //stosunek wykonanych obliczeń (ammount of trials to ammount of succesfull trials) printf("Program szacuje wartosc liczby pi metoda Monte Carlo na podstawie Igly Buffona. Prosze Pamietac o:\n-Podaniu liczby prób jako wartości większej od zera i całkowitej\n-Podaniu długości igły mniejszej niż szerokości między liniami\n"); printf("Prosze podac liczbe prob:\n"); scanf("%lf", &n); printf("Prosze podac dlugosc igly:\n"); scanf( "%lf", &l); printf("Prosze podac szerokość miedzy liniami:\n"); scanf("%lf", &L); if (n <= 0.0 || l > L) { printf("Nastąpił błąd wynikających z podania niepoprawnych danych\n"); return 1; } printf("%f %f %f", n,l,L); //debuging midway for (i=0; i<n; i++) { x = ((double)rand()/RAND_MAX)*3.1415; // losowy kat w radianach (random angle in radians) d = ((double)rand()/RAND_MAX)*L/2; // losowa dlugosc igly mniejsza niz odstep miedzy liniami (random needle lenght) y = (l/2) * sin (x); printf("Próba%d\n", i); if (d<=y) { P++; } } stosun=n/P; printf("Stosun %lf", stosun); pi=stosun*2*l/L; printf("liczba pi=%lf", pi); return 0; } 

Answer 1

The variance of Monte Carlo simulations tends to be very large. (They converge to zero as n increases, but very slowly.) You can read about Buffon Needle, in particular here . I would like to draw your attention to the number on which 5 different tracks of 1000 volumes are plotted. Having examined each of them, it will look as if the value is already converging, but to a completely wrong (and different in each case) result.

Good when you start your analysis. The error estimate carried out strictly in the above source is about 5.6 / n. This is the variance, therefore, to get the standard deviation, we need a square root, which makes it even worse. After 1000 volumes, you can expect π to have an error of around 0.073 (although it can be twice as high or much smaller, just by sheer chance). After a million throws, it should be correct to approximately ± 0.002. After a hundred million, you are likely to get the first three digits after the decimal point on the right and the first error on the fourth or (rarely) fifth. I don't know your n, but the values ​​you get look good in that which is typical. Edit: with n = 1,000,000, they are equal to +3 σ and -1 σ.

Note that this is a common misconception that you can get around this by doing several runs and averaging their results (which should fluctuate above and below the correct value with the same probability, in the end). This, of course, works, but does not provide any benefits, except for the obvious that it can be extended to several machines. The error in executing 100 simulations out of a thousand iterations in each case is exactly the same as when starting one with one hundred thousand iterations.

The obvious way to make the rating better is thus increasing the number of hits. One of the great advantages of Monte Carlo methods is that individual samples are usually completely independent of each other and, thus, provide great opportunities for parallelization (OpenMP, vectorized instructions, graphics processor). But the first step is to really try to make your code faster in order to increase n by an order of magnitude or less. For example, since you know from analysis that accuracy will never be truly overwhelming, you can use float instead of double , which is twice as fast on some platforms. Eliminate the sin call (this is of course more elegant if you don't need the π value to start with). Display a summary only every 1000 steps or so.

Less obvious is the analysis doing the right thing and noticing how you can reduce the numerator in the variance (the denominator will be n what you do). If I'm not mistaken, this experiment works best when the length of the needle is equal to the distance between the lines: less than this is suboptimal, and more than this leads to several intersections that you do not want. There may be other ways to improve it by changing the distribution of the angle and / or position.