166
6 Computing Integrals and Testing Code
f_mean = f_mean + f(x[i], y[j])
f_mean = f_mean/num_inside
area = num_inside/(n**2)*(x1 - x0)*(y1 - y0)
return area*f_mean
(See the file MC_double.py.)
Verification A simple test case is to check the area of a rectangle [0, 2] × [3, 4.5]
embedded in a rectangle [0, 3] × [2, 5]. The right answer is 3, but Monte Carlo
integration is, unfortunately, never exact so it is impossible to predict the output of
the algorithm. All we know is that the estimated integral should approach 3 as the
number of random points goes to infinity. Also, for a fixed number of points, we can
run the algorithm several times and get different numbers that fluctuate around the
exact value, since different sample points are used in different calls to the Monte
Carlo integration algorithm.
The area of the rectangle can be computed by the integral
2
0
4.5
3 dydx, so in
this case we identify f (x, y) = 1, and the g function can be specified as (e.g.)
1 if (x, y) is inside [0, 2] × [3, 4.5] and −1 otherwise. Here is an example, using
samples of different sizes, on how we can utilize the MonteCarlo_double function
to compute the area:
In [1]: from MC_double import MonteCarlo_double
In [2]: def g(x, y):
...:
return (1 if (0 <= x <= 2 and 3 <= y <= 4.5) else -1)
...:
In [3]: MonteCarlo_double(lambda x, y: 1, g, 0, 3, 2, 5, 100)
Out[3]: 2.9484
In [4]: MonteCarlo_double(lambda x, y: 1, g, 0, 3, 2, 5, 1000)
Out[4]: 2.947032
In [5]: MonteCarlo_double(lambda x, y: 1, g, 0, 3, 2, 5, 1000)
Out[5]: 3.0234600000000005
In [6]: MonteCarlo_double(lambda x, y: 1, g, 0, 3, 2, 5, 2000)
Out[6]: 2.9984580000000003
In [7]: MonteCarlo_double(lambda x, y: 1, g, 0, 3, 2, 5, 2000)
Out[7]: 3.1903469999999996
In [8]: MonteCarlo_double(lambda x, y: 1, g, 0, 3, 2, 5, 5000)
Out[8]: 2.986515
Note the compact if-else construction in the definition of g. It is a one-line
alternative to, for example,
if (0 <= x <= 2) and (3 <= y <= 4.5):
return 1
else:
return -1
From the output, we see that the values fluctuate around 3, a fact that supports a
correct implementation, but in principle, bugs could be hidden behind the inaccurate
answers.
6 Computing Integrals and Testing Code
f_mean = f_mean + f(x[i], y[j])
f_mean = f_mean/num_inside
area = num_inside/(n**2)*(x1 - x0)*(y1 - y0)
return area*f_mean
(See the file MC_double.py.)
Verification A simple test case is to check the area of a rectangle [0, 2] × [3, 4.5]
embedded in a rectangle [0, 3] × [2, 5]. The right answer is 3, but Monte Carlo
integration is, unfortunately, never exact so it is impossible to predict the output of
the algorithm. All we know is that the estimated integral should approach 3 as the
number of random points goes to infinity. Also, for a fixed number of points, we can
run the algorithm several times and get different numbers that fluctuate around the
exact value, since different sample points are used in different calls to the Monte
Carlo integration algorithm.
The area of the rectangle can be computed by the integral
2
0
4.5
3 dydx, so in
this case we identify f (x, y) = 1, and the g function can be specified as (e.g.)
1 if (x, y) is inside [0, 2] × [3, 4.5] and −1 otherwise. Here is an example, using
samples of different sizes, on how we can utilize the MonteCarlo_double function
to compute the area:
In [1]: from MC_double import MonteCarlo_double
In [2]: def g(x, y):
...:
return (1 if (0 <= x <= 2 and 3 <= y <= 4.5) else -1)
...:
In [3]: MonteCarlo_double(lambda x, y: 1, g, 0, 3, 2, 5, 100)
Out[3]: 2.9484
In [4]: MonteCarlo_double(lambda x, y: 1, g, 0, 3, 2, 5, 1000)
Out[4]: 2.947032
In [5]: MonteCarlo_double(lambda x, y: 1, g, 0, 3, 2, 5, 1000)
Out[5]: 3.0234600000000005
In [6]: MonteCarlo_double(lambda x, y: 1, g, 0, 3, 2, 5, 2000)
Out[6]: 2.9984580000000003
In [7]: MonteCarlo_double(lambda x, y: 1, g, 0, 3, 2, 5, 2000)
Out[7]: 3.1903469999999996
In [8]: MonteCarlo_double(lambda x, y: 1, g, 0, 3, 2, 5, 5000)
Out[8]: 2.986515
Note the compact if-else construction in the definition of g. It is a one-line
alternative to, for example,
if (0 <= x <= 2) and (3 <= y <= 4.5):
return 1
else:
return -1
From the output, we see that the values fluctuate around 3, a fact that supports a
correct implementation, but in principle, bugs could be hidden behind the inaccurate
answers.