
Show HN: Pi Approximation using Monte Carlo - krat0sprakhar
http://montepie.herokuapp.com
======
gus_massa
Nice visualization! Some UI comments

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".

~~~
krat0sprakhar
Thanks for the great suggestions. Will add these in soon!

------
krat0sprakhar
Hi HN!

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!

[0] -
[http://www.cs.princeton.edu/~dpw/courses/cos326-12/](http://www.cs.princeton.edu/~dpw/courses/cos326-12/)

[1] - [https://github.com/prakhar1989/ocaml-
experiments/blob/master...](https://github.com/prakhar1989/ocaml-
experiments/blob/master/cos326/a1/a1.ml#L204-L216)

~~~
fphilipe
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.

[http://bl.ocks.org/fphilipe/6094255](http://bl.ocks.org/fphilipe/6094255)

~~~
krat0sprakhar
Wow! Nicely done. The idea of showing relative error is pretty cool :D

------
MrQuincle
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...

~~~
afafsd
Would you care to estimate the magnitude of this error, and how many steps you
would need before it actually became relevant?

~~~
darkmighty
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.

~~~
MrQuincle
Thanks!

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.

------
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?

~~~
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).

------
wfunction
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](http://pastebin.com/ezi7S35N) (Edit: just
noticed a typo -- I should be initializing the numerator/denominator inside
the loop, not outside, sorry.)

What's the point of adding randomness?

~~~
adwn
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.

~~~
wfunction
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.

~~~
keithpeter
_" 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...

[http://sohcahtoa.org.uk/pages/maths_montecarlo.html](http://sohcahtoa.org.uk/pages/maths_montecarlo.html)

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).

~~~
jwmerrill
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.

~~~
jwmerrill
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.

~~~
keithpeter
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!

------
Houshalter
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.

------
killercup
Funny, I just recently saw this method as a proposal for an example of Rust's
built-in `rand` library: [https://github.com/rust-
lang/rust/pull/16362](https://github.com/rust-lang/rust/pull/16362)

------
skykooler
I notice that it seems to overshoot more than it undershoots. Any idea why
this is?

~~~
lelandbatey
Interestingly, for me it was under-shooting vs. overshooting.

------
kken
Small remark: The constant you are calculuating is usually denoted as a lower
case pi, not upper case.

~~~
krat0sprakhar
Haha! Thanks for pointing this out.

------
scott_wilson46
f you are interested. I did something similar for fpga:
[https://github.com/scottwilson46/FPGARandom](https://github.com/scottwilson46/FPGARandom)

