I saw something like this in C++ as an introduction to Monte Carlo, so I decided to make something similar in Python. My original code used for loops, but I vectorized it with no small amount of effort, and it now runs orders of magnitude faster. For example, I can calculate pi to .002% error with 100,000,000 randomized coordinates in approximately 15 seconds. Careful to start small, because memory fills up quickly, and the computer will run slow if you overstep your RAM. I'm able to go up to a bit more than 150 million without compromising speed and functionality in other tasks.
For those who are curious, vectorization is a technique whereby numpy (or similar) arrays replace things like loops. They're a bit tricky to write (at least for me), but they work beautifully.
It might be useful for visualization to plot the occurrence of data points, and observe the randomness
Download
Copy to clipboard1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 | from random import *
from numpy import *
darts=input("How many darts? ")
x=random.random(darts)
y=random.random(darts)
one=ones(darts)
dist=sqrt(x*x+y*y)
print "dist ",dist
hit=(less_equal(dist,one)==True)
print "hit ",hit
hits=sort(hit)
print "hits ",hits
hitsnum=hits.searchsorted(True, side="right") - hits.searchsorted(True, side="left")
print "hitsnum ",hitsnum
pic=4.0*hitsnum/darts
print "Pi calculated as ",pic
error=100*(pic-3.1415926535897932)/(3.1415926535897932)
print "Error is ", error, "%"
|
This is a good intro to monte carlo methods, and it helped me grasp some of the concepts.
As a next step, one might try to compute the definite integral of a function experimentally. Create an array of random x-values, plug them into a weird function, and figure out what percentage fall inside the area of the curve. Multiply by the total area of the plot (the domain of the random values) and you've got your area.
