Hacker News new | past | comments | ask | show | jobs | submit login
Monte Carlo estimation of Pi (github.com/calebmadrigal)
82 points by calebm on March 14, 2017 | hide | past | favorite | 42 comments

BTW, your scatter plots aren't colourblind-friendly - they just look like boxes of indistinguishable dots to me, I had to ask a friend if I was missing something. You seem to have accidentally made an Ishihara test[0], thanks to the particular shades of red and green you've used - red & blue would be better for accessibility.

[0] https://en.m.wikipedia.org/wiki/Ishihara_test

Is there a service that makes it trivially easy to simulate colour blindness on an image or document? Heck, why's it not built into data publishing applications?

Yes, Vischeck: http://www.vischeck.com/vischeck/.

There's a Python library that's usable at the command line, too: https://github.com/joergdietrich/daltonize.

I found a program called simdaltonism on the mac app store that is perfect for this.

Yap, I was thinking "I'm only seeing some random points" too...

Thanks for the heads-up!

One of the first labs I did in college, it brings back pleasant nostalgia.

It was my first independent Monte Carlo calculation -- changed my life.

That project led directly to the 'rand' call in this function [1] 13 years later and fundamentally altered the way I think about data analysis.

[1] https://github.com/4kbt/PlateWash/blob/master/mlib/bootstrap...

Yup - oldie but a goodie. We did it with rice grains and graph paper recently, then a spreadsheet.


Someone here suggested the Steinmetz solid as a follow up which I will try next time I teach that group.

Yep, it was a problem I encountered had back in high school Pascal. Fun times :)

Compare this program, which uses a text object in its own source code to estimate pi: http://ioccc.org/1988/westley.fix.c

Yeah, and you can just make the circle bigger to get a better approximation. :P

Wow. IOCCC spins my head yet again.

I prefer Buffon's needle method as a 'physical' way of calculating pi.

Beer would probably help the calculation.

"Does that look random to you?" https://youtu.be/M34TO71SKGk?t=1m20s


I remember a Mythbusters episode where they were building a buttered-toast-dropper. To test the null hypothesis, they first did 10 runs without buttered toast, just marking one side with an X. They got 7/10 landing the same way up... and declared "It should be 5/10! This is not random enough!"...

One hour Radiolab worth listening to:


the intro to this video explained the calculation very well (up to 1 minute)

Once an Amazon interviewer asked me this (in my second or third year of University) as an interview question.

I didn't get it at the time and I deeply resented the problem. What does this have to do with being a software engineer?

But interesting to see how it's actually done.

It's something I'd expected someone to know in anything related to data science / machine learning / big data. But if it wasn't a position specific to those fields, then yeah that's that a pretty silly riddle of a question.

edit: even within those fields, it's definitely more of a very-specific-domain-knowledge type of question and not something that I think is worth the time to ask it in a job interview.

This is a common thing to learn during your first week or so in an introductory computer science course. They were probably picking something that a second-year student could use to easily demonstrate basic coding ability. Unfortunately, this didn't happen to be part of your intro course.

Monte carlo method can be used to provide software project estimates. The examples I've seen using it are about the same as planning poker

Planning poker has no element of randomness though? I'm confused.

Some addendum:

1. You can do one quadrant rather than four, then multiply that approximation by 4.

2. Montecarlo is good to develop an intuition about Monte carlo approximations but converges really slowly.

Ramanujan-Sato series for pi or Chudnovsky algorithm can be a scalable alternative, gives you a lot of digits.

Now... how many digits you usually need? Very few. Practical applications requiring more than 15 digits are hard to find. e.g: NASA uses 15 digits for interplanetary navigation, even for the most distant spacecraft such as the Voyager 1 which is now 20 billion kilometers from Earth.

Can just use a constant, or in the worst case, if missing, define it as:

    pi = arc tan (1) * 4

Accuracy from one quadrant and four quadrants would be the same, though, and full circle is easier to code (no extra conditions).

To get better accuracy you would need to cut out some parts that are clearly inside or outside the circle, and only throw the points in the areas where the circle can actually be.

Sure, that is equivalent to doing it on one quadrant with 4 times the dots and averaging it.

You can try the same with quasi-random rather than pseudo-random numbers (Sobol, Halton, or even just additive). Should probably get better accuracy with the same number of iterations.

Interestingly, this method is the primary example on the Dart homepage: https://www.dartlang.org/.

BTW there is a nice related course on coursera: "Statistical Mechanics: Algorithms and Computations"[1]. Also notice there is a rosettacode entry for this[2].

1. https://www.coursera.org/learn/statistical-mechanics

2. http://rosettacode.org/wiki/Monte_Carlo_methods

Funny this is here, I did recently a dwitter here to visually do the same: https://www.dwitter.net/d/702

It is also the example I give on a (true) parallel library for front-end javascript with math-intensive code such as calculating PI with Monte Carlo: https://github.com/franciscop/uwork

Today I tried to write this on shell, I couldn't do it in a compact way. Are there some nice PLs for these kind of plays? Cleanest I could do was doing below in a single line via "python -c". I couldn't get it done via bc or awk.

  from random import random
  l = lambda (a, b), c: (a + c, b + 1)
  g = (1 if random()**2 + random()**2 <= 1 else 0 for i in range(10**6))
  r, t = reduce(l, g, (0.0, 0))
  print r / t * 4

  python3 -c 'from random import random; N=1000000;print(4*sum([ int(random()**2 + random()**2 <= 1.0) for i in range(N)])/N)'

F# version that fits in a tweet:

  let r,n=Random(),1e7
  let d()=(r.NextDouble()-0.5)**2.0
  |>Seq.sumBy(function _ when d()+d()<0.25->4.0/n|_->0.0)
  |>printfn "%f"

You could try R. This ends up being a pretty straightforward one-liner:

  sum(runif(10^6, 0, 1)^2 + runif(10^6, 0, 1)^2 < 1) / 10^6 * 4

If for some reason you'd like to actually Run this notebook,I uploaded it to:


You can click on it, sign in to run, edit, etc.

Sign in: any outlook, xbox, hotmail, or msft ID.

https://www.youtube.com/watch?v=I-BC_vI4CAE - I thought this was genius, obtaining a value of Pi using the monte carlo method, via rain drops!

LOL we all have those moments, I thought I was a genius for thinking I could calculate PI by doing a double integral over x,y from 0 to 1 using the Heaviside step function with (x^2+y^2) or something I can't remember... Basically it was that in the integral the step function is on if x^2+y^2 is less than 1 and off if >1, so you would end up with PI/4. Anyways it didn't work out :)

Props to Ken Liu for having the character Luan Zya derive the Monte Carlo estimator in "The Grace of Kings". That was a wonderful read.

Guidelines | FAQ | Lists | API | Security | Legal | Apply to YC | Contact