I am fairly sure that I understand the how Monte Carlo integration works but I am not understanding the formulation of how it is used to estimate Pi. I am going by the procedure outlined in the 5th slide of this presentation http://homepages.inf.ed.ac.uk/imurray2/teaching/09mlss/slides.pdf

I understand the preliminary steps. Pi is equal to 4 times the area of a quarter of the unit circle. And the area of the top-right quarter of the unit circle centered at (0,0) is equivalent to the integral of the curve that is the top-right quarter of the unit circle in $0<x<1$ and $0<y<1$.

What I don't understand is how this integral is

$\iint I((x^2+y^2)<1)P(x,y)dxdy$

where $P(x,y)$ is uniformly distributed in the unit square around the quarter circle (i.e. it is always equal to 1 if $0<x<1$ and $0<y<1$ and 0 otherwise). So this would mean that $I((x^2+y^2)<1)P(x,y)$
is the function that is the top-right quadrant of the unit circle at $0<x<1$ and $0<y<1$ but I do not understand how this is true since the indicator function can only be 1 or 0. I understand that it is probably written in this way to make Monte Carlo sampling easy (i.e. it is an expectation so just sample from $P(x,y)$ and get the average of the samples applied to $I((x^2+y^2)<1)$) but it just does not make intuitive sense to me why that integral represents the area under that curve.

Could someone provide an intuitive explanation of this. Maybe show how that integral was derived in a step-by-step way?


I was able to gain a better understanding by relating the expectation to an area. I will explain it here in case it helps anyone. First start with relating Pi to the area of the top-right quadrant of the unit circle

$\pi=4\times A_{tr}$

Then we place the top-right quadrant into the unit square. And under a uniform distribution over the unit square, the area of the circle quadrant is proportional to the probability of obtaining a sample from it. It follows that the following equality holds


and $A_{square}=1$ so


And substituting into the original equation

$\pi=4\times P(x^2+y^2<1)$

and it is also true that $P(x^2+y^2<1)=E[I(x^2+y^2<1)]$ which is equal to the original double integral.

So I understood it by relating the area to a probability then relating that probability to an expectation that is equivalent to the integral. Let me know if I have made any mistakes.


The area of a circle circle of radius $l$ is equal to $\pi l^2$. It means that a quarter of circle has area $l^2\pi/4$. This means that the square with side the radius of the circle as $area=l^2$.

This means that the ratio between the area of a quarter of circle and the area of the square is $\pi/4$.

A point $(x,y) $ is in the square if $ 0<x<1, 0<y<1$. and it is in the quarter of circle if $ 0<x<1, 0<y<1 ,x^2+y^2<1$.

Your integral is so $∬I((x^2+y^2)<1)P(x,y)= ∬I((x^2+y^2)<1) I(0<x<1)I(0<y<1)$ That is exactly the area described by a quarter of circle

enter image description here

  • $\begingroup$ I guess I'm just having a tough time drawing a connection between the terms inside the integral and the curve itself. If you plotted I(x^2+y^2<1)I(0<x<1)(0<y<1) for different values of x and y, you would not get the curve. Why is that? $\endgroup$ May 25 '14 at 23:15
  • 1
    $\begingroup$ $\left\lbrace (x,y):(x^2+y^2<1), (0<x<1), (0<y<1) \right\rbrace$ are the points on the quarter of circle. I suggest you to try to plot this points $\endgroup$
    – Donbeo
    May 25 '14 at 23:18
  • $\begingroup$ I agree with that. But when you apply the indicator function I(.), they all get pushed to either 1 or 0. $\endgroup$ May 26 '14 at 0:48
  • $\begingroup$ What do you mean? $\endgroup$
    – Donbeo
    May 26 '14 at 7:18
  • 1
    $\begingroup$ The indicator function in an integral is just another way to define the curves where to compute the integral. $\int_{quarter~ of~ circle} = \int 1(x^2+y^2<1)*1(0<x<1)*1(0<y<1) $ $\endgroup$
    – Donbeo
    May 26 '14 at 13:47

The simplest intuitive explanation relies on understanding that $E(I(A)) = P(A)$. Thus, $\int \int I(x^2+y^2 < 1)dxdy = P(x^2 + y^2 < 1)$. Once you realize the double integal is simply a probability, it should make intuitive sense that you could sample $x$ and $y$ from the unit square and compute the proportion of draws for which $x^2 + y^2 <1$.

Perhaps the other piece of intuition missing from your understanding is the connection between area and probability. Since the area of the entire unit square is 1 and points $(x,y)$ are uniformly distributed within the square, the area of any region $A$ within the unit square would correspond to the probability that a randomly chosen point would be within $A$.

  • $\begingroup$ This is the way I understand it as well. But I'm having trouble connecting it to the formulation Pi=4x(area of quarter circle). It does not really make intuitive sense to compare areas to samples. I suppose the connection is that under a uniform distribution, the number of samples is proportional to the area. $\endgroup$ May 25 '14 at 23:10
  • 1
    $\begingroup$ @user1893354 Answer revised. Let me know if that helps your intuition. $\endgroup$
    – jsk
    May 25 '14 at 23:50

I landed on this surfing CV, and I see that the code of the Monte Carlo is in Octave. I happen to have a simulation in R that makes the idea of deriving the number $\pi$ as a bivariate uniform distribution in the $[0,1]$ plane under the constraints of the integrals in the OP very intuitive:

Given that the quarter of a circle is enclosed in a 1-unit square, the area is $\pi/4$. So generating uniformly distributed points in the square $(x,y)$ will end up carpeting the entire square, and calculating the fraction fulfilling $1 < \sqrt{(x^2+y^2)}$ will be tantamount to integrating $∬\textbf{1}((x^2+y^2)<1) \,\textbf{1}(0<x<1)\,\textbf{1}(0<y<1)$ since we are just selecting the fraction of dots within the circle in relation to the unit square:

x <- runif(1e4); y <- runif(1e4)
radius <- sqrt(x^2 + y^2)
# Selecting those values within the circle is obtained with radius[radius < 1]:
(pi = length(radius[radius < 1]) / length(radius)) * 4     =    3.1272

We can plot the values falling within the radius among 10,000 draws:

enter image description here

And we can, naturally, get closer and closer approximation by selecting more points. With 1 million points we get:

(pi = length(radius[radius < 1]) / length(radius)) * 4 [1] 3.141644

a very approximate result. Here is the plot:

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.