5 ms·

Show HN: Pi Approximation using Monte Carlo

• deleted 12y ago
  [deleted]
• krat0sprakhar 12y ago
  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/cos326/a1/a1.ml#L204-L216 https://github.com/prakhar1989/ocaml-experiments/blob/master...
  • Goopplesoft 12y ago
    Thanks for the share. Tried this in python using some higher iterations and started getting good numbers around 1E6 Link: http://nbviewer.ipython.org/gist/pnegahdar/eb62475b28d5b394d824 http://nbviewer.ipython.org/gist/pnegahdar/eb62475b28d5b394d...
  • fphilipe 12y ago
    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 12y ago
      Wow! Nicely done. The idea of showing relative error is pretty cool :D
  • nemetroid 12y ago
    > 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/ http://jsfiddle.net/qmw4d9p9/1/
• clemsen 12y ago
  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 12y ago
    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 12y ago
      It would be interesting to see how quickly a quasi-random, low-discrepancy sequence converges to PI (e.g. a Sobol Sequence).
  • teisman 12y ago
    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 12y ago
    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.)
• kken 12y ago
  Small remark: The constant you are calculuating is usually denoted as a lower case pi, not upper case.
  • krat0sprakhar 12y ago
    Haha! Thanks for pointing this out.
• skykooler 12y ago
  I notice that it seems to overshoot more than it undershoots. Any idea why this is?
  • lelandbatey 12y ago
    Interestingly, for me it was under-shooting vs. overshooting.
• wfunction 12y ago
  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?
  • robzyb 12y ago
    > There is nothing random about pi. That is exactly what makes Monte Carlo neat! There is nothing random about pi, yet using randomness we can arrive at pi.
    • sherjilozair 12y ago
      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.
      • wfunction 12y ago
        Calling it "pseudo-Monte Carlo" is like calling a landline phone a "wired cell phone". It kind of misses the point.
        • retroencabulato 12y ago
          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.
          • wfunction 12y ago
            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.
• MrQuincle 12y ago
  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 12y ago
    Would you care to estimate the magnitude of this error, and how many steps you would need before it actually became relevant?
    • darkmighty 12y ago
      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 12y ago
        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.
  • darkmighty 12y ago
    The line is a 0-measure space: it shouldn't matter in theory but if your boundary is "squiggly" enough it indeed does, in practice.
  • stephencanon 12y ago
    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.
• deleted 12y ago
  [deleted]
• scott_wilson46 12y ago
  f you are interested. I did something similar for fpga: https://github.com/scottwilson46/FPGARandom https://github.com/scottwilson46/FPGARandom
• killercup 12y ago
  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
• gus_massa 12y ago
  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 12y ago
    Thanks for the great suggestions. Will add these in soon!
• Houshalter 12y ago
  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.