Finding PI digits using Monte Carlo

Question

I have tried many algorithms for finding π using Monte Carlo. One of the solutions (in Python) is this:

def calc_PI():
    n_points = 1000000
    hits = 0

    for i in range(1, n_points):
        x, y = uniform(0.0, 1.0), uniform(0.0, 1.0)

        if (x**2 + y**2) <= 1.0:
            hits += 1

    print "Calc2: PI result", 4.0 * float(hits) / n_points

The sad part is that even with 1000000000 the precision is VERY bad (3.141...).

Is this the maximum precision this method can offer? The reason I choose Monte Carlo was that it's very easy to break it in parallel parts. Is there another algorithm for π that is easy to break into pieces and calculate?

Answer 1

Your fractional error goes by sqrt(N)/N = 1/sqrt(N), So this is a very inefficient way to get a precise estimate. This limit is set by the statistical nature of the measurement and can't be beaten.

You should be able to get about floor(log_10(N))/2-1 digits of good precision for N throws. Maybe -2 just to be safe...

Even at that it assumes that you are using a real RNG or a good enough PRNG.

Answer 2

This is a classic example of Monte Carlo. But if you're trying to break the calculation of pi into parallel parts, why not just use an infinite series and let each core take a range, then sum the results as you go?

http://mathworld.wolfram.com/PiFormulas.html

Answer 3

Use a quasi random number generator (http://www.nag.co.uk/IndustryArticles/introduction_to_quasi_random_numbers.pdf) instead of a standard pseudo RNG. Quasi random numbers cover the integration area (what you're doing is a MC integration) more evenly than pseudo random numbers, giving better convergence.