4

I was referring to this lecture http://videolectures.net/mlss09uk_murray_mcmc/. However, I didn't get how the pi value was calculated. Here is a screenshot

enter image description here

As far as I know the integral gives the area under the curve. So I didn't get how the double integral is giving the value of pi using P(x,y).

Also how did it convert to mean in the octave. I am just a beginner. So please some guidance is required

user34790
  • 6,757
  • 10
  • 46
  • 69
  • 2
    This is two questions: (1) Why does a double integral express an area? (2) How can an area be estimated by an average? The second question is essentially the same as your immediately previous one at http://stats.stackexchange.com/questions/39128. The first is a matter of mathematics, where you might hope for good explanations on the math site--but it really comes down to that's the definition of area. – whuber Oct 10 '12 at 19:56
  • Note that double integral is volume under a surface, so in that formula, they compute volume of cylinder with unit height. – sitems Oct 10 '12 at 20:17
  • 1
    Forget about probability. The red area inside the square is $1/4$ the area of a circle of radius $1$ (the side of the square). What is the value of the red area? – Zen Oct 10 '12 at 20:40

1 Answers1

4

With respect to the double integral $4 \iint I(x^2 + y^2 < 1) P(x, y) dx~dy$, we can easily see that this evaluates to $\pi$ if we actually perform the calculation. First notice that the function P is 1 only on the inside of the unit square in the upper right quadrant. Consequently, we can rewrite the integral as $4 \int_0^1 \int_0^1 I(x^2 + y^2 < 1) dx~dy$.Now we notice that the indicator function is 1 when ever the coordinate (x, y) is inside the unit circle. This means that multiplying the two values will yield 1 when ever (x, y) is in the upper right quadrant of the unit circle, and will be 0 everywhere else. We can therefore transform our double integral into $4 \int_0^1 \int_0^{\sqrt{1 - y^2}} 1~dx~dy = 4 \int_0^1 \sqrt{1-y^2} dy = 4 \frac{\pi}{4} = \pi$.

If we think carefully about what the indicator and P are doing, intuitively it should be obvious what is going on. These two functions will be 1 in the red area of your diagram. If we think about what the area of a quarter of a unit circle is, it should be obvious that it is $\pi/4$ since $\pi~r^2 = \pi~1^2 = \pi$. This integral is just a round about way of calculating the area of a quarter of the unit circle.

With respect to the octave code, if we go through it piece by piece everything should be clear.

First we have a=rand(S, 2). All this does is generate S samples from the unit square in the upper right quadrant. In other words, it generates samples from the square shown in his diagram. Note that P will evaluate to 1 for all the S samples in a.

Next we have this mysterious expression a.*a. All this does is multiply each entry in the matrix by itself. The sum with the 2 says to add each row in a, so in effect, sum(a.*a, 2) says compute $x^2 + y^2$ for each of my S samples. The < operator will return a 1 if true or a 0 if false, so it basically serves like our indicator function in the integral. In other words, the result of sum(a.*a, 2) < 1 will be a bunch of 1's and 0's, where a 1 indicates that a sample was inside the upper right quadrant of the unit circle and a 0 indicates that the sample was outside the unit circle (but still inside the unit square of course).

The key observation is that the ratio of 1's to 0's in this sampling scheme will be identical to the integral $\int_0^1 \int_0^1 I(x^2 + y^2 < 1) dx~dy$. Thus all we have to do is take the mean of the 1's and 0's, which is exactly the ratio we want. Finally we multiply by 4, and we have $\pi$. Thus the code 4*mean(sum(a.*a, 2)<1) really does approximate the integral we care about, which as we just showed, is exactly the value of $\pi$.

jlund3
  • 1,158
  • If it would help to see this exact same algorithm in a less obtuse language like Python, I'm more than willing to add to my answer... – jlund3 Oct 10 '12 at 21:19