Why isn't my python programm converging to pi?

117 views Asked by At

I wrote this method:

def approx_pi(n, p):
    """
    Approximates Pi by putting dots in a square and counting
    the dots in the max possible circle in that square.
    :param n: Number of dots
    :param p: Precision used for urandom
    :return: Approximation of Pi
    """
    in_circle = 0
    k = 100
    max_int = (2**(8*p)-1)
    for i in range(n):
        # get two random Numbers between 0 and 1
        x = int.from_bytes(os.urandom(p), byteorder='big')/max_int
        y = int.from_bytes(os.urandom(p), byteorder='big')/max_int

        if x ** 2 + y ** 2 <= 1:
            in_circle += 1

        # Just for debugging
        if (i+1) % k == 0:
            k = int(k*1.01)
            print(i, '\t',4*in_circle/i)
            sys.stdout.flush()

    return 4*in_circle/n

In most of my Testruns, it is getting stable at 3.141 and then diverges around that value. Is this a weakness of urandom? But if so, why isn't pi shifting in one direction? Or is there something wrong with my code.

2

There are 2 answers

1
ivangtorre On

First of all. Your code is very confusing. It is better if you call the variables by a name that give you an idea about what they mean.

I have simplified the code to the necessary lines. Hence, I use random library for generating the points position.

The problem with this method is that it converges very slowly. With 10**8 points i got 3.14167604.

import random


def approx_pi(points):
    """
    Approximates Pi by putting dots in a square and counting
    the dots in the max possible circle in that square.
    :param points: Number of dots
    :return: Approximation of Pi
    """
    in_circle = 0
    for dummy_i in xrange(points):
        # get two random Numbers between 0 and 1
        x_dot_position = random.random()
        y_dot_position = random.random()

        if x_dot_position ** 2 + y_dot_position ** 2 <= 1:
            in_circle += 1

    return 4.0*in_circle/points

EDITING You were rigth about the function random.random(). So in the next code I have used random.random_integers which is a uniform distribution between [low, high] values.

The code is parallelized and i tried for 10**10 points getting:

PI = 3.14157765 Tiempo de calculo = 3790 seconds

import multiprocessing
import time
import numpy as np

starting_point = time.time()

def approx_pi(point):
    """
    Approximates Pi by putting dots in a square and counting
    the dots in the max possible circle in that square.
    :param points: Number of dots
    :return: Approximation of Pi
    """
    # get two random Numbers between 0 and 1
    x_dot_position = float(np.random.random_integers(0,10**10))/10**10
    y_dot_position = float(np.random.random_integers(0,10**10))/10**10

    if x_dot_position ** 2 + y_dot_position ** 2 <= 1:
        return 1
    else: 
        return 0

###########################################################

total_points     = 1*10**10       
paso    = 1*10**8       

in_circle = 0

for in_este_bucle in xrange(0, total_points, paso):
    print "Procesadores disponibles: " + str(multiprocessing.cpu_count())

    pool = multiprocessing.Pool()
    resultado = pool.map(approx_pi, xrange(in_este_bucle, in_este_bucle + paso))
    pool.close()
    pool.join()
    in_circle += sum(resultado)
    del resultado


print 'Total time: ' + str(time.time()-starting_point) +' seconds'
print
print 4.0*in_circle/total_points
0
ethanpowell On

I'm guessing you're starting to dabble in Monte Carlo methods. When using Monte Carlo methods the random numbers you generate must have a uniform distribution to yield an accurate approximation. Try using:

import random

random.uniform(0,1)

This should yield better results. ivangtorre is absolutely right that a large sample size is needed for an accurate approximation. Here's a link to Monte Carlo methods on wikipedia: https://en.wikipedia.org/wiki/Monte_Carlo_method The introduction section of this article talks about exactly what you're working on now. Cheers.

Edit: You could also increase the range of random.uniform(0,1) and then adjust your algorithm accordingly to get a slight increase in performance.