Monte Carlo simulation of $\pi$

Question

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

Answer 1

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)

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