integration method

Question

I want to calculate an integral by using the hit and miss method. I can not understand how this method works. I would be grateful if someone could explain me and help me to calculate the value, with a realistic and simple example as
$I=\int_{0}^{1} x^2dx$

or anything you want. Thank you very much for your concern, in advance.

Answer 1

Your integral is the area $A$ under the graph of $y = x^2$ and above $y=0$ for $x$ from $0$ to $1$:

The total area of the square $0 \le x \le 1$, $0 \le y \le 1$ is $1$. If you choose a random point $(X,Y)$ in this square, the probability that this point will be in the red region will be the fraction of the area of the square that is red, namely $A$. Now imagine choosing $n$ random points $(X_j, Y_j)$, independently, in this square. Each has probability $A$ of being in the red region, so if $n$ is large the number $R$ of points in the red region will probably be close to $n A$. Therefore we can use $R/n$, the fraction of our points that fall in the red region, as an estimate for $A$.

Answer 2

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.