3
$\begingroup$

I am trying to find the value of $\pi$ using Monte Carlo simulation. However, I don't want to generate two random numbers as coordinates. Instead, I want to select a point on the edge of the square first and select another point on the line that connects the origin and this previously selected point.

In other words, the point of this is to find the area of the circle without using x and y coordinates.

The idea is that if $a$ is a random number between 1 and $\sqrt{2}$ and then $b$ is a random number between 0 and $a$ In theory all values of $b$ that are less than 1 should be inside the circle, so I should have the number of points in the circle and the number of total points, that way I can calculate the fraction that is proportional to $\pi$.

I hope this image can explain better than my words

I wrote the following Python code to achieve that but I get results close to $3.4$ instead of the real value of $\pi$. I asked the question on Stack Overflow and they said the issue is not about the program itself but the way I am trying to find the value. Is there a flaw in my reasoning?

Here's the code I used:

$\pi$ Monte Carlo

import random
import pylab


def MonteCarloPI(numtries):
    circle = 0
    for i in range (numtries):
        a = (float(random.random())+1)**0.5
        b = float(random.random())*a
        if b <= 1:
            circle += 1
        rapportoAree = (circle/numtries)
    return rapportoAree*4

print(MonteCarloPI(1000))
$\endgroup$
9
  • $\begingroup$ Could you post this to StackOverflow, as it's much more on-topic there: stackoverflow.com $\endgroup$
    – rolando2
    Jun 6, 2018 at 19:20
  • 1
    $\begingroup$ @rolando2 I did, there I got advised of posting this here. I think there is something wrong with my statistics, not with my code, take a look at the picture, are my logics wrong there? $\endgroup$
    – user210809
    Jun 6, 2018 at 19:22
  • 1
    $\begingroup$ I think the lengths of those lines that connect the origin and the chosen point on the edge of the square are not uniformly distributed but your algorithm assumes they are. If I select 1 million points (with step size 1/1m, not random) on an edge and calculate the length I see a highly right skewed distribution. $\endgroup$
    – ayhan
    Jun 6, 2018 at 19:38
  • 1
    $\begingroup$ Your method estimates $8 (\sqrt{2} - 1)\approx3.31371$ rather than $\pi.$ $\endgroup$
    – whuber
    Jun 6, 2018 at 20:34
  • 1
    $\begingroup$ In an unrelated question, I think the line "rapportoAree = (circle/numtries)" should be outside the for loop. This doesn't affect the result, but it would avoid dividing a thousand times. $\endgroup$
    – Pere
    Jun 6, 2018 at 20:39

1 Answer 1

4
$\begingroup$

As user2285236 comments, the problem is that the points aren't uniformly distributed in the eighth of a square as they are supposed to be. I assume your first step intends to choose a random point in the edge (although I think you missed an exponent), and that works. However, in the second step you choose a random point in the line to the center and then you get more points in the regions where those lines are closer.

Please forgive me for answering that Python question with R code:

> library(ggplot2)
> 
> n <- 1000
> r1 <- runif(n)
> a <- sqrt(r1^2+1)
> b <- runif(n)*a
> inside <- b<1
> sum(inside)/n*4
[1] 3.556
> y <- b/a
> x <- r1*b/a
> qplot(x,y)

Your random points

Those are the random points you get with your code. Please notice that points aren't evenly distributed and they get closer in the lower part of the triangle, that is, inside the circle. Therefore you get a value larger than pi as points are more likely to fall inside the circle.

Suppressing the "2" exponent in a doesn't change much.

On the other hand, if points are actually uniformly distributed, you get a good approximation of pi:

> n <- 1000
> sum(runif(n)^2+runif(n)^2<1)/n*4
[1] 3.148

Interestingly, you could also compute pi with your algorithm, but you would need to weight the points to account for the non uniform distribution. The weight is just b/a:

> n <- 1000
> r1 <- runif(n)
> a <- sqrt(r1^2+1)
> b <- runif(n)*a
> inside <- (b<1)*b/a
> total <- b/a
> sum(inside)/sum(total)*4
[1] 3.155418
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.