The method you want is very much like throwing darts. You create random points inside a rectangle that you know the area of and count how many of them are underneath the curve that you're trying to integrate.
For example, in Matlab we might do this to get the curve we're trying to integrate:
>> x = 0:0.01:1; % Creates a list [0, 0.01, ... 0.99, 1] >> plot(x,x.^2)
which gives the following plot:

We then need to generate some random points to scatter throughout the plot:
>> N = 100; >> xpts = rand(N,1); >> ypts = rand(N,1); >> hold on % Plot over the top of our last plot >> plot(xpts,ypts,'.')

Now you need to know which of the points fall under the curve (and let's plot them too, for fun)
>> under = ypts < xpts.^2; >> plot(xpts(under),ypts(under),'r.')

The vector under
is now a vector of 1s wherever the point (x,y) is under the curve and 0s when it is above the curve. To approximate the area under the curve we find the average of this vector (with the mean
function) and multiply it by the area of the rectangle (which is 1 in this case, but might not be in general).
>> area = 1 * mean(under); >> disp(area) 0.3800
We know that the exact area is 1/3, so this isn't too bad an approximation.
If you wanted to find out something about the variance of the approximation, you could write a loop that does this 1000 times, to give you 1000 different estimates of the area, and look at some of its statistics:
>> for i = 1:1000 >> X = rand(N,1); >> Y = rand(N,1); >> I(i) = mean(Y < X.^2); >> end
You can look at the mean, variance and standard deviation of I
:
>> mean(I) ans = 0.3321 >> var(I) ans = 0.0022 >> std(I) ans = 0.0469
So the mean is close to 1/3, the variance is close to the theoretical value of (1/3) * (2/3) / 100 = 0.00222... and the standard deviation is around 0.05, which means that your estimate with 100 points will be within 0.23 and 0.43 about 95% of the time. By using more points you could make this much more accurate, although obviously it would be slower.