[clemsen]

Nice visualization. After playing with it for a while, the number never really converged to Pi. The value even reached 3.135 after 10000 iterations. Do you have any measure how "stable" this method is? Or how many iterations are needed to be relatively certain that the first n digits won't change anymore?

[teisman]

The simulation actually follows a binomial distribution, meaning the standard deviation can be calculated as sqrt(npq), where n is the number of trials, p is the probability of success (pi/4), and q is 1-p. For a simulation of 10.000 iterations, the standard deviation is approximately 41. You say you observed 3.135 after 10.000 iterations, meaning you observed about 7837 hits. The expectation was to observe 7854 hits, meaning there was only a difference of 17 hits - well within 1 standard deviation.

edit: naturally you can also use this information to calculate the number of trials needed get a desired probability that the estimation falls within certain bounds from pi.

[Karellen]

I got 3.1544 after 10000 iterations, which also didn't match my intuition of how quickly it "should" converge. So I decided to knock up a quick deterministic version in C that just iterated over all points in [0, 1), [0, 1) at various resolutions.

For 10x10 (100 iterations), I got 3.44

For 100x100 (10,000 iterations), I got 3.1796

For 1000x1000 (1,000,000 iters), I got 3.14552

For 10,000x10,000 (100,000,000 iters), I got 3.141990

For 100,000x100,000 (10,000,000,000 iters, 2m06s runtime), I got 3.141633.

So the deterministic version takes a long time to converge on the right answer too. Therefore, it seems that the Monte Carlo version appears less accurate than you'd think, not because of the stability of the randomness, but because your intuition (and mine) is way off about how good any kind of sampling is at converging on the right answer.

See @teisman's great reply for figuring out how many iterations you'd need for a given level of accuracy.

Edit: Looking back at the numbers, you need 100x as many samples to roughly get an extra digit of precision. That's 10x the iterations on each axis, which makes sense that that would give you an extra digit of precision in your answer. Therefore, after 10000 (100x100) iterations, you wouldn't expect better accuracy than 2 digits, or 3.1. Which is what we see.

Edit2: I'm an idiot. After converting my test from (sqrt((x * x + y * y)/(samples * samples)) < 1) to ((x * x + y * y) < (samples * samples)) and making the maths integer only, I got the 100,000x100,000 runtime down to 10.8s. (Same answer though.)

[krat0sprakhar]

I tried the OCaml version which was much much faster than JS. Had to run around 100,000+ simulations to get 3.14xxx consistently. This is a really simple algorithm so I'd suggest giving it a shot in a your favorite language to see how fast it converges.

[brotchie]

It would be interesting to see how quickly a quasi-random, low-discrepancy sequence converges to PI (e.g. a Sobol Sequence).