I think the reset button should be smaller (you usually don't want to reset the simulation).
I expected that the number in the upper right corner is the number of additional simulations to do. For example, if n=1000 and box=200, the system makes 200 additional simulations and now n=1200.
Also, I'd like to see an advanced mode: (but perhaps I'm not a typical user)
I'd like to see a graph of pi_aproximation vs n, so I can see how the value "converges" to pi.
And I'd like to see an error estimation, for example, the expected variation is ~sqrt(n)/4 (I think that the correct constant is not 4, but it's a similar small number). So you can draw error bars or the pi+-sqrt(n)/4 "bounds".
This is a small hack that I created over the weekend. I learnt about the Monte Carlo method while doing Princeton's Functional Programming course[0] in OCaml[1] as one of the assignment problems. To get the first four digits consistently, I ran around 200,000 - 400,000 simulations - not quite possible in the browser!
It's not the browser (or JS) that's slow, it's all the animations and stuff. If you write Javascript to do only what the OCaml code does, it runs fast enough: http://jsfiddle.net/qmw4d9p9/1/
Nice work. I did a very similar simulation a year ago when I learned about MC in the context of particle filtering in a course about recursive estimation at ETH.
Nicely visualized. One thing however that I realized is that conceptually you actually need a second step if a point falls exactly on the line with distance 1 to the origin. Half of these points should count to the circle, not all of them, or none of them. With 32 bit floating point values there won't be many of these of course...
Every primitive Pythagorean triple has c odd, so strictly speaking (+/-1,0) and (0,+/-1) are the only points with binary FP coordinates that have distance exactly 1 from the origin (decimal FP is more interesting in this regard).
Of course, there are more points with binary FP coordinates whose magnitude computed naively as sqrt(xx + yy) produces a rounded result of 1.0 in binary floating-point; (3/5,4/5) is one such example.
I would address this by using careful computation of the residual of xx + yy - 1.0 to determine which side of the circle the point lies on; something like the following sketch:
float residual = x*x + y*y - 1
if (residual == 0) {
/* x*x + y*y rounded to 1; we need to do a more careful check
to determine whether we're inside or outside. WLOG, assume
that |x| >= |y|. */
if (fabsf(x) < fabsf(y)) { float tmp = x; x = y; y = tmp; }
/* We know that 0.5 <= x*x < 1.0, from which it follows that
fmaf(x,x,-1) is computed without any rounding; we can then
do another fma to isolate the residual (this is also exact
because we know that y*y cannot be *too* tiny unless x is
exactly 1.0, in which case the first fma is exactly zero). */
residual = fmaf(y,y,fmaf(x,x,-1));
}
if (residual < 0) { /* inside circle */ }
else if (residual == 0) { /* on circle */ }
else { /* outside circle */ }
It's worth noting that fma isn't actually necessary here; it suffices to evaluate xx + yy in double precision, except in the |x| = 1 case, which is easy to detect. fma is conceptually cleaner (and faster on modern architectures that support it), however.
I believe the error should be (roughly) ~ size of f.p. interval at 1 * pi ~ pi* 2^(-32/64), i.e. barely representable. Number of steps would be ~ 1/(error^2) ~ 2^(2* 32/64).
Even giving 50% probability at the line the result would be minutely biased because of many other things involving the discreteness of the floating point arithmetic. A less biased version would have to maybe increase the resolution around the boundary as number of steps progresses.
Practically, it might be fun to have a small library that artificially reduces the precision of floating point numbers. I've never heard about such a thing though.
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?
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.
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.)
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.
I will never understand why people think Monte Carlo is such a neat solution. There is nothing random about pi.
Why didn't you just pick the center of the square, use that to partition it into 4 pieces, and recurse for each subsquare?
Or just as easily do this? http://pastebin.com/ezi7S35N (Edit: just noticed a typo -- I should be initializing the numerator/denominator inside the loop, not outside, sorry.)
Monte Carlo methods are never "seriously" used to approximate pi – there are algorithms which converge much, much faster.
This here, however, is not a "serious" attempt at computing pi, but a learning exercise for someone who just found out about Monte Carlo methods.
Regarding the randomness of the problem to solve: A practical application for Monte Carlo methods is integration over a high-dimensional space (with dozens or more degrees of freedom). Traditional deterministic methods have a runtime which is exponential in the number of dimensions, while for Monte Carlo integration, the error of the result decreases as 1/sqrt(N) (where N is the number of samples, which is proportional to the runtime), independent of the dimensionality of the integration space.
I think you missed my point. I know this isn't how pi is calculated, that is completely irrelevant to what I've been trying to say. I also know Monte Carlo has other benefits, and that's entirely my point: the actual benefits aren't illustrated in the example.
what I'm saying is that the entire point of illustrating an algorithm by example is to illustrate the power of the algorithm compared to a naive approach, regardless of whether or not there is a better way to solve the example problem.
i.e., the example should motivate the algorithm.
But computing pi is one of the worst possible illustrations of why anyone would use Monte Carlo, because it's inferior to the naive approach. i.e., it doesn't motivate why anyone would want to use Monte Carlo.
"But computing pi is one of the worst possible illustrations of why anyone would use Monte Carlo, because it's inferior to the naive approach. i.e., it doesn't motivate why anyone would want to use Monte Carlo."
Any suggestions for a simulation that does illustrate the use of Monte Carlo methods while still being explainable with 16+ age range non-specialist mathematics? I'm hacking around with two step simulations like a tree diagram...
My crack at saying how slowly a monte-carlo simulation 'converges' to a value of pi (not really converging, just confidence intervals tightening around the value).
I like computing the volume of the intersection of 3 orthogonal cylinders as a less trivial example problem. It's still easy to compute whether a given point is inside the volume, and there's still an analytic solution to compare to, but the problem is complicated enough that you start to appreciate how much easier MC is than other approaches.
Or compute the area of the mandlebrot set. This one has the advantage of being highly non-convex, which makes it hard to decide how to implement a recursive subdivision routine, since you don't know whether some of the target area is inside a given square just from looking at a few points on its boundary.
Mandlebrot set a bit over the maths level for intended group, but other repeated maps like logistic equation might be OK. I'll go with the cylinders first!
Sounds nice. Easy to visualise and close enough to a crucible design (toroidal volume intersected by a cube off axis and oblique to the central axis of the toroid) I once heard of in the days of Fortran.
The advantage of Monte Carlo isn't that it's the best method, it's that it's so simple. You just sample random points and see how many are in the circle.
Do not underestimate the value of these visualizations. They reach 1000s of kids, students, and who-knows, adults fooling around on Sundays.
So, perhaps you can turn this around, and ask for a special visualization that would be better in your opinion. In that case you might be lucky and have a person like krat0sprakhar nicely visualizing your specific problem.
And just my two cents: most problems I encounter are integration (full Bayesian) and the detection of extremes (MAP). There is nothing random about these problems either, maybe even less than pi (which might be "normal", nobody knows yet)!
I'd be happy to see a visualization of another nice, simple problem, that is more up your alley!
You're right, as a method for computing pi it's silly and inefficient.
As a way of demonstrating the utility of Monte Carlo algorithms in general, though, it's a neat and easily visualisable example, which is why it's used so often in elementary introductions of the concept.
Actual examples of how Monte Carlo algorithms can be useful involve the sampling of many-dimensional spaces via algorithms which sample points you're interested in more often than points you're not interested in -- check out the Metropolis-Hastings algorithm for an example of useful Monte Carlo in action.
We don't necessarily have to use randomness. We can take equally-spaced points in the quarter as well, rather than sampling. This is called pseudo-Monte Carlo.
Could you substantiate this with any references? I thought hard about this, and I find no reason to believe that equally-spaced points would give a biased estimator.
The comment has a point. It is possible that the particular deterministic rule that chooses points could interact with the shape (quarter circle).
In this case, the error would not necessarily go to zero as the number of points increases without bound.
In the particular case of a grid based deterministic probe, and a quarter-circle target, it seems clear that this would not happen.
But consider another example where the underlying target was "all points with rational coordinates". All the probes in a grid sampling scheme would hit the target, but the target has measure zero.
Incidentally, the idea of using a deterministic, low variability sequence for sampling is called quasi Monte Carlo (http://en.wikipedia.org/wiki/Quasi-Monte_Carlo_method). It can give almost order 1/n convergence, much better than the 1/sqrt(n) convergence possible with ordinary Monte Carlo.
> But consider another example where the underlying target was "all points with rational coordinates". All the probes in a grid sampling scheme would hit the target, but the target has measure zero.
And how exactly would the Monte Carlo version differ here? Every single random number generated is going to be rational isn't it?
It's just a smart-aleck comment pointing out that there is a "bias" because we're using a finite number of points and we get the same exact answer every time -- and since that answer isn't exactly equal to pi, the procedure is "biased".
Which no one cares about because (1) this isn't a probabilistic process in the first place, so "bias" itself is meaningless, and (2) there is now an exact bound on the error, unlike the previous case (I guess you can think of this as the "variance" if you want).
Who cares? The error is going to be there either way. And Monte Carlo can't guarantee any hard bound on that error, whereas using equally spaced points gives you a pretty darn good hard bound.
No, you missed the point. In problems with a very large sample space, quasi-MC methods are invaluable. They converge faster and can be shown to approximately have the benefits of traditional MC methods.
I think we're talking past each other here. Talking about "sample space" inherently assumes you're dealing with a random process. My point is pi isn't random, so coming up with a random process just so you can force MC on it and then suddenly deciding you want to make it deterministic just so you can call it "pseudo-MC" is kinda... weird, to say the least.
I really wish people would stop trying to find ways of calling their algorithms "Monte Carlo" just to make their work sound fancy. All it really achieves is it makes them avoid thinking about how to solve the problem. Monte Carlo isn't a universal hammer. First you're supposed to show randomness actually helps you gain something, then you're supposed to start using it.
I think there is some value to this method, but it may not be readily apparent in the present context. If you have a problem that you don't yet have a firm theoretical handle on you can figure out a baseline value for the answer using this technique, and then work backwards to give you the formula that allows you to arrive at the real answer.
I use that trick a lot, for instance with Euler problem 307 I could not find the answer easily in an analytical way, did the simulation, found an approximate answer (to within 6 decimal places or so), used that to figure out what I was doing wrong and then computed the real answer.
Some may see that as cheating, but it works wonders and not just on that particular problem.
I've done that too, and I think there's nothing wrong with that. All I'm saying is people should stop calling it Monte Carlo once they've found the correct solution and realized it has nothing to do with randomness.
I did this a long time ago in Lua thinking it would be really cool to approximate pi. But it never converged to the correct value. No matter how many times I ran it or for how long. I think it was 3.15 or something like that, I don't remember. My best guess it was due to floating point rounding, or because the random function was imperfect.
I think the reset button should be smaller (you usually don't want to reset the simulation).
I expected that the number in the upper right corner is the number of additional simulations to do. For example, if n=1000 and box=200, the system makes 200 additional simulations and now n=1200.
Also, I'd like to see an advanced mode: (but perhaps I'm not a typical user)
I'd like to see a graph of pi_aproximation vs n, so I can see how the value "converges" to pi.
And I'd like to see an error estimation, for example, the expected variation is ~sqrt(n)/4 (I think that the correct constant is not 4, but it's a similar small number). So you can draw error bars or the pi+-sqrt(n)/4 "bounds".