Thursday, 27th August, 2009
Calculus time!
We all know that Monte Carlo simulations are useful, but what can they really be used for… Aside from designing nukes and H-Bombs, and predicting the stock market, and modeling risk for companies of many types.
Well, as it turns out integration.
Let’s start with an implementation of Riemann Sums.
import simulator
import random
def reimann_integral(f, lower, upper, runs):
randGen = random.Random()
randF = (lambda : randGen.uniform(lower, upper))
randCon = (lambda x: x)
sim = simulator.Simulator()
sim.setFunction(f)
sim.setNumberTests(runs)
sim.setRandomConstraint(randCon)
sim.setRandomSource(randF)
sim.run()
start = (str([lower]),f(lower))
Area = 0
for x in sorted(sim.results()):
x0 = float(start[0].lstrip("[").rstrip("]"))
x1 = float(x[0].lstrip("[").rstrip("]")) # Get x0, x1 out of the result set
xdash = abs(x0-x1) # Width of the strip
A = xdash*x[1] # Find the area of the strip
Area += A
start = x
return Area
Basically, it finds lots of points at random, measures the distance between them on the x-axis and multiplies them by the heights and adds them all up!
Let’s see if it works, let’s find the integral of x^2 + 3x - 4 between 1 and 3.
print riemann_integral((lambda x : (x**2 + 3*x - 4)), 1, 3 10000) 12.6670208327
A quick bit of algebra shows that this integral is 38/3, or 12.66… exactly! We’re pretty close :-D. But what about functions that vary very quickly, we’ll need to account for the rapid change in height of the variable. I think that given enough steps, a trapezoidal integral will help with this, let’s try implementing it. It should have nearly the exact same form as the Riemann sum, we’ll just calculate the area of each trapezoid differently.
def trapezoidal_integral(f,lower,upper,runs):
randGen = random.Random()
randF = (lambda : randGen.uniform(lower, upper))
randCon = (lambda x: x)
sim = simulator.Simulator()
sim.setFunction(f)
sim.setNumberTests(runs)
sim.setRandomConstraint(randCon)
sim.setRandomSource(randF)
sim.run()
start = (str([lower]),f(lower))
Area = 0
for x in sorted(sim.results()):
x0 = float(start[0].lstrip("[").rstrip("]"))
x1 = float(x[0].lstrip("[").rstrip("]")) # Get x0, x1 out of the result set
height = float(abs(x0-x1)) # height of the trapezoid
A = (0.5)*height*(start[1]+x[1]) # Find the area of the trapezoid
Area += A
start = x
return Area
Et Voila! Ca Va, eh? Quick test…
print trapezoidal_integral((lambda x: (x**2+3*x-4)), 1,3,10000) 12.6662471122
More significant digits! The wonders of technology. The trapezoidal rule should almost always be more accurate than the Riemann sums version.
Done! :-D