Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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

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.

share|improve this question
up vote 8 down vote accepted

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

share|improve this answer
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? – user1893354 May 25 '14 at 23:15
$\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 – Donbeo May 25 '14 at 23:18
I agree with that. But when you apply the indicator function I(.), they all get pushed to either 1 or 0. – user1893354 May 26 '14 at 0:48
What do you mean? – Donbeo May 26 '14 at 7:18
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) $ – 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$.

share|improve this answer
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. – user1893354 May 25 '14 at 23:10
@user1893354 Answer revised. Let me know if that helps your intuition. – jsk May 25 '14 at 23:50
Thank you that helps. – user1893354 May 26 '14 at 0:48

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

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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