integration method

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.

-

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$.

-
thank you very much, but my main problem is how i could construct an algorithm in order to calculate an intergal, by using the "hit or miss" method. – johan paul Nov 24 '11 at 19:22
Some detail was given in a comment I made to an earlier version of the same question. You will want a pseudo-random number generator that produces random integers with uniform distribution on $[0,1]$. Some calculators have this built in. Having generated two of them, say $X$ and $Y$, if $Y \le X^3$, call that a hit, otherwise call it a miss. Repeat $n$ times where $n$ is fairly large. (We need a "new" random pair each time.) If the number of hits is $H$, then $H/n$ should give you a reasonable estimate of the area. – André Nicolas Nov 24 '11 at 19:40
@Kostas The variance is given by $\frac{A (1-A)}{n}$ where $A$ is the area of the red portion (= the value of the integral). This is always less than or equal to $\frac{1}{4n}$. – Srivatsan Nov 24 '11 at 19:48
In general, if you use the idea with another function, and a rectangle instead of the square that we got because $1^3=1$, let $p$ be the area of region of interest divided by area of rectangle. Then variance is $p(1-p)/n$. As mentioned by Srivatsan, $p(1-p)\le 1/4$ always. In the case of a rectangle of base $a$ height $b$, generate random numbers $U$, $V$ in $[0,1]$ as usual, let $X=aU$, $Y=bV$, and test whether $(X,Y)$ is in the region whose area you want. If $H$ is the number of times you get a hit, your estimate for the area is $(H/n)ab$. – André Nicolas Nov 24 '11 at 20:16
I would like to thank you all very much, for your help and concern. Everything you told me was quite understandable and useful. Thanks again. – johan paul Nov 24 '11 at 20:31

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.

-
Note that convergence is almost always faster if you use a low-discrepancy sequence instead of random points (the quasi-Monte-Carlo method.) – user7530 Nov 24 '11 at 20:22
I would like to thank you all very much, for your help and concern. Everything you told me was quite understandable and useful. Thanks again. – johan paul Nov 24 '11 at 20:31
Chris, may I ask you something about your answer? I was wondering how did you find the 95% interval confidence, (0.23,0.43). Since the std is close to 0.05, the 95% CI shouldn't be (0.28,0.38), or i do something wrong ? – johan paul Nov 26 '11 at 13:27
A 95% confidence interval is about two standard deviations either side of the mean. Here the standard deviation is 0.05, so twice that is 0.1, which is where I got [0.23, 0.43] as your 95% interval. Note that the 95% rule is only true if you can approximate the thing you're measuring by a normal distribution. To see if you can make the normal approximation, multiply the number of trials $N$ by the smallest of $p$ and $1-p$ and see if that's bigger than 5. Here you have $N=100$, $p=1/3$ and $1-p=2/3$, so you multiply $100\times 1/3 = 33$ which is greater than $5$, so the normal approx is okay. – Chris Taylor Nov 26 '11 at 18:58
Thanks again Chris. – johan paul Nov 29 '11 at 23:11

$\newcommand{\+}{^{\dagger}}% \newcommand{\angles}[1]{\left\langle #1 \right\rangle}% \newcommand{\braces}[1]{\left\lbrace #1 \right\rbrace}% \newcommand{\bracks}[1]{\left\lbrack #1 \right\rbrack}% \newcommand{\ceil}[1]{\,\left\lceil #1 \right\rceil\,}% \newcommand{\dd}{{\rm d}}% \newcommand{\down}{\downarrow}% \newcommand{\ds}[1]{\displaystyle{#1}}% \newcommand{\equalby}[1]{{#1 \atop {= \atop \vphantom{\huge A}}}}% \newcommand{\expo}[1]{\,{\rm e}^{#1}\,}% \newcommand{\fermi}{\,{\rm f}}% \newcommand{\floor}[1]{\,\left\lfloor #1 \right\rfloor\,}% \newcommand{\half}{{1 \over 2}}% \newcommand{\ic}{{\rm i}}% \newcommand{\iff}{\Longleftrightarrow} \newcommand{\imp}{\Longrightarrow}% \newcommand{\isdiv}{\,\left.\right\vert\,}% \newcommand{\ket}[1]{\left\vert #1\right\rangle}% \newcommand{\ol}[1]{\overline{#1}}% \newcommand{\pars}[1]{\left( #1 \right)}% \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\pp}{{\cal P}}% \newcommand{\root}[2][]{\,\sqrt[#1]{\,#2\,}\,}% \newcommand{\sech}{\,{\rm sech}}% \newcommand{\sgn}{\,{\rm sgn}}% \newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}} \newcommand{\ul}[1]{\underline{#1}}% \newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert}$ Indeed, we $\large never$ do Montecarlo with $\ds{\int_{0}^{1}x^{2}\,\dd x}$. We make the following transformation such that the integrand becomes closer to a smooth function: $$\int_{0}^{1}x^{2}\,\dd x =\half\bracks{\int_{0}^{1}x^{2}\,\dd x + \int_{0}^{1}\pars{1 - x}^{2}\,\dd x} =\int_{0}^{1}\bracks{x\pars{x - 1} + \half}\,\dd x$$

$\large\mbox{You can use this little}\quad \verb=C++=\quad \mbox{script}$:

#include <cstdlib>
#include <iostream>
using namespace std;
const double RANDMAX1=double(RAND_MAX) + 1.0;
const unsigned long long ITERATIONS=1000000ULL; // For example, one million.

int main()
{
double result=0,x;

for ( unsigned long long n=0 ; n<ITERATIONS ; ++n ) {
x=rand()/RANDMAX1;
result+=x*(x - 1.0) + 0.5;
}

result/=ITERATIONS;

cout<<"Result = "<<result<<endl;

return 0;
}

-
both are 'smooth' functions. Can you elaborate on the rationale for doing this? – Memming Feb 1 '14 at 16:58
@Memming What I mean is 'closer to a constant function'. This is a usual trick in Montecarlo integration: We try to integrate something which does not change abruptly. Maybe, 'smooth' was not a clear word. English is not my native language. Sorry. – Felix Marin Feb 1 '14 at 18:43
This transformation gives better convergence properties for the algorithm described in your code (which is essentially Gaussian quadrature at random locations) but that's not what I understand by the "hit and miss" method (generating random $(x,y)$ pairs and checking to see if they fall into the region being integrated). – Chris Taylor Feb 1 '14 at 19:12
@ChrisTaylor I agree. They are different but it serves for the same purpose. Thanks. – Felix Marin Feb 1 '14 at 19:27