6 ms·

Monte Carlo estimation of Pi

82 points by calebm 10y ago

pimmen89 10y ago
One of the first labs I did in college, it brings back pleasant nostalgia.
- keithpeter 10y ago
  Yup - oldie but a goodie. We did it with rice grains and graph paper recently, then a spreadsheet. http://sohcahtoa.org.uk/pages/maths_montecarlo.html http://sohcahtoa.org.uk/pages/maths_montecarlo.html Someone here suggested the Steinmetz solid as a follow up which I will try next time I teach that group.
- calebm 10y ago
  Yep, it was a problem I encountered had back in high school Pascal. Fun times :)
- ISL 10y ago
  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/bootstrapData.m https://github.com/4kbt/PlateWash/blob/master/mlib/bootstrap...
bingaling 10y ago
With darts: https://www.youtube.com/embed/M34TO71SKGk https://www.youtube.com/embed/M34TO71SKGk
- forinti 10y ago
  Beer would probably help the calculation.
- woof 10y ago
  "Does that look random to you?" https://youtu.be/M34TO71SKGk?t=1m20s https://youtu.be/M34TO71SKGk?t=1m20s Golden!
  - vacri 10y ago
    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!"...
    - woof 10y ago
      One hour Radiolab worth listening to: http://www.radiolab.org/story/91684-stochasticity/ http://www.radiolab.org/story/91684-stochasticity/
- contravariant 10y ago
  I prefer Buffon's needle method as a 'physical' way of calculating pi.
- jps359 10y ago
  the intro to this video explained the calculation very well (up to 1 minute)
charlescearl 10y ago
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.
zht 10y ago
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.
- fiftyacorn 10y ago
  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
  - gkop 10y ago
    Planning poker has no element of randomness though? I'm confused.
- TallGuyShort 10y ago
  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.
- daveFNbuck 10y ago
  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.
mirajshah 10y ago
Interestingly, this method is the primary example on the Dart homepage: https://www.dartlang.org/ https://www.dartlang.org/.
schoen 10y ago
Compare this program, which uses a text object in its own source code to estimate pi: http://ioccc.org/1988/westley.fix.c http://ioccc.org/1988/westley.fix.c
- josephcooney 10y ago
  Wow. IOCCC spins my head yet again.
- pmiller2 10y ago
  Yeah, and you can just make the circle bigger to get a better approximation. :P
grondilu 10y ago
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 https://www.coursera.org/learn/statistical-mechanics 2. http://rosettacode.org/wiki/Monte_Carlo_methods http://rosettacode.org/wiki/Monte_Carlo_methods
kd5bjo 10y ago
See also https://m.youtube.com/watch?v=RZBhSi_PwHU https://m.youtube.com/watch?v=RZBhSi_PwHU
anfractuosity 10y ago
https://www.youtube.com/watch?v=I-BC_vI4CAE 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!
- seiferteric 10y ago
  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 :)
auselen 10y ago
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
- seiferteric 10y ago
  python3 -c 'from random import random; N=1000000;print(4*sum([ int(random()**2 + random()**2 <= 1.0) for i in range(N)])/N)'
- daveFNbuck 10y ago
  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
- endeavour 10y ago
  F# version that fits in a tweet: let r,n=Random(),1e7 let d()=(r.NextDouble()-0.5)**2.0 [0.0..n] |>Seq.sumBy(function _ when d()+d()<0.25->4.0/n|_->0.0) |>printfn "%f"
EvilTerran 10y ago
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 https://en.m.wikipedia.org/wiki/Ishihara_test
- Waterluvian 10y ago
  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?
  - _ihaque 10y ago
    Yes, Vischeck: http://www.vischeck.com/vischeck/ http://www.vischeck.com/vischeck/. There's a Python library that's usable at the command line, too: https://github.com/joergdietrich/daltonize https://github.com/joergdietrich/daltonize.
  - asafira 10y ago
    I found a program called simdaltonism on the mac app store that is perfect for this.
  - scoot 10y ago
    Colorblinding - Chrome extension: https://chrome.google.com/webstore/detail/colorblinding/dgbgleaofjainknadoffbjkclicbbgaa?hl=en https://chrome.google.com/webstore/detail/colorblinding/dgbg...
- paulofalcao 10y ago
  Yap, I was thinking "I'm only seeing some random points" too...
- calebm 10y ago
  Thanks for the heads-up!
franciscop 10y ago
Funny this is here, I did recently a dwitter here to visually do the same: https://www.dwitter.net/d/702 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 https://github.com/franciscop/uwork
partycoder 10y ago
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
- pps43 10y ago
  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.
  - partycoder 10y ago
    Sure, that is equivalent to doing it on one quadrant with 4 times the dots and averaging it.
smortaz 10y ago
If for some reason you'd like to actually Run this notebook,I uploaded it to: https://notebooks.azure.com/smortaz/libraries/MonteCarloPi https://notebooks.azure.com/smortaz/libraries/MonteCarloPi You can click on it, sign in to run, edit, etc. Sign in: any outlook, xbox, hotmail, or msft ID.
pps43 10y ago
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.

82 points by calebm 10y ago

pimmen89 10y ago

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

keithpeter 10y ago
Yup - oldie but a goodie. We did it with rice grains and graph paper recently, then a spreadsheet. http://sohcahtoa.org.uk/pages/maths_montecarlo.html http://sohcahtoa.org.uk/pages/maths_montecarlo.html Someone here suggested the Steinmetz solid as a follow up which I will try next time I teach that group.
calebm 10y ago
Yep, it was a problem I encountered had back in high school Pascal. Fun times :)
ISL 10y ago
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/bootstrapData.m https://github.com/4kbt/PlateWash/blob/master/mlib/bootstrap...

bingaling 10y ago

With darts: https://www.youtube.com/embed/M34TO71SKGk https://www.youtube.com/embed/M34TO71SKGk

forinti 10y ago
Beer would probably help the calculation.
woof 10y ago
"Does that look random to you?" https://youtu.be/M34TO71SKGk?t=1m20s https://youtu.be/M34TO71SKGk?t=1m20s Golden!
- vacri 10y ago
  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!"...
  - woof 10y ago
    One hour Radiolab worth listening to: http://www.radiolab.org/story/91684-stochasticity/ http://www.radiolab.org/story/91684-stochasticity/
contravariant 10y ago
I prefer Buffon's needle method as a 'physical' way of calculating pi.
jps359 10y ago
the intro to this video explained the calculation very well (up to 1 minute)

charlescearl 10y ago

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.

zht 10y ago

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.

fiftyacorn 10y ago
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
- gkop 10y ago
  Planning poker has no element of randomness though? I'm confused.
TallGuyShort 10y ago
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.
daveFNbuck 10y ago
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.

mirajshah 10y ago

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

schoen 10y ago

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

josephcooney 10y ago
Wow. IOCCC spins my head yet again.
pmiller2 10y ago
Yeah, and you can just make the circle bigger to get a better approximation. :P

grondilu 10y ago

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 https://www.coursera.org/learn/statistical-mechanics 2. http://rosettacode.org/wiki/Monte_Carlo_methods http://rosettacode.org/wiki/Monte_Carlo_methods

kd5bjo 10y ago

See also https://m.youtube.com/watch?v=RZBhSi_PwHU https://m.youtube.com/watch?v=RZBhSi_PwHU

anfractuosity 10y ago

https://www.youtube.com/watch?v=I-BC_vI4CAE 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!

seiferteric 10y ago
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 :)

auselen 10y ago

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

seiferteric 10y ago
python3 -c 'from random import random; N=1000000;print(4*sum([ int(random()**2 + random()**2 <= 1.0) for i in range(N)])/N)'
daveFNbuck 10y ago
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
endeavour 10y ago
F# version that fits in a tweet: let r,n=Random(),1e7 let d()=(r.NextDouble()-0.5)**2.0 [0.0..n] |>Seq.sumBy(function _ when d()+d()<0.25->4.0/n|_->0.0) |>printfn "%f"

EvilTerran 10y ago

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 https://en.m.wikipedia.org/wiki/Ishihara_test

Waterluvian 10y ago
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?
- _ihaque 10y ago
  Yes, Vischeck: http://www.vischeck.com/vischeck/ http://www.vischeck.com/vischeck/. There's a Python library that's usable at the command line, too: https://github.com/joergdietrich/daltonize https://github.com/joergdietrich/daltonize.
- asafira 10y ago
  I found a program called simdaltonism on the mac app store that is perfect for this.
- scoot 10y ago
  Colorblinding - Chrome extension: https://chrome.google.com/webstore/detail/colorblinding/dgbgleaofjainknadoffbjkclicbbgaa?hl=en https://chrome.google.com/webstore/detail/colorblinding/dgbg...
paulofalcao 10y ago
Yap, I was thinking "I'm only seeing some random points" too...
calebm 10y ago
Thanks for the heads-up!

franciscop 10y ago

Funny this is here, I did recently a dwitter here to visually do the same: https://www.dwitter.net/d/702 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 https://github.com/franciscop/uwork

partycoder 10y ago

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

pps43 10y ago
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.
- partycoder 10y ago
  Sure, that is equivalent to doing it on one quadrant with 4 times the dots and averaging it.

smortaz 10y ago

If for some reason you'd like to actually Run this notebook,I uploaded it to: https://notebooks.azure.com/smortaz/libraries/MonteCarloPi https://notebooks.azure.com/smortaz/libraries/MonteCarloPi You can click on it, sign in to run, edit, etc. Sign in: any outlook, xbox, hotmail, or msft ID.

pps43 10y ago

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.