7
$\begingroup$

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.

3 Answers 3

1

$\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 $ enter image description here

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

 #include  #include  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
  • 0
    @ChrisTaylor I agree. They are different but it serves for the same purpose. Thanks.2014-02-01
18

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:

quadratic

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,'.') 

points

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.') 

pointsred

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.

  • 0
    This answer from over two and a half years ago suddenly got a flurry of upvotes - anyone know why?2014-02-01
9

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

enter image description here

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

  • 2
    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.2011-11-24