为什么我的 python 程序没有收敛到 pi？答案

【问题标题】：Why isn't my python programm converging to pi?为什么我的 python 程序没有收敛到 pi？
【发布时间】：2015-08-31 10:35:34
【问题描述】：

我写了这个方法：

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

在我的大多数测试运行中，它在 3.141 处稳定下来，然后在该值附近发散。这是urandom的弱点吗？但如果是这样，为什么 pi 不朝一个方向移动？还是我的代码有问题。

【问题讨论】：

如果您使用x = int.from_bytes(os.urandom(p), byteorder='big')（y 模拟）和if x ** 2 + y ** 2 <= max_int**2，它是否有效？只是为了消除不精确的浮点运算。
n 和 p 的价值观是什么？
@Stefan，我选择 n := 1.000.000.000 和 p \in {1, 2, 3, 4}。
@koffein 这真的是一个好点！不幸的是，到目前为止它似乎没有任何效果，但我正在增加数字并进行更多测试。
@Finn 还有一件事：选择n>max_int**2 没有任何好处，因为在这种情况下，您还可以采用所有可能的(x,y)-组合，您会得到更好的结果。

标签： python random statistics pi

【解决方案1】：

首先。您的代码非常混乱。最好用一个能让您了解它们的含义的名称来调用变量。

我已将代码简化为必要的行。因此，我使用随机库来生成点位置。

这种方法的问题是收敛速度非常慢。 10**8 分我得到 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

编辑你对函数 random.random() 是正确的。所以在接下来的代码中，我使用了 random.random_integers，它是 [low, high] 值之间的均匀分布。

代码是并行化的，我尝试获得 10**10 分：

PI = 3.14157765 计算时间 = 3790 秒

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

【讨论】：

谢谢。我认为我们不应该使用random.random()，主要是因为返回的值在[0, 1)，所以我们会得到一个 pi' > pi。根据 python 文档，它的保护程序可以使用os.urandom()。

【解决方案2】：

我猜您已经开始涉足蒙特卡洛方法了。使用蒙特卡罗方法时，您生成的随机数必须具有均匀分布才能产生准确的近似值。尝试使用：

import random

random.uniform(0,1)

这应该会产生更好的结果。 ivangtorre 是绝对正确的，即需要大样本量才能进行准确的近似。这是 wikipedia 上 Monte Carlo 方法的链接：https://en.wikipedia.org/wiki/Monte_Carlo_method 本文的介绍部分准确地讨论了您现在正在做什么。干杯。

编辑：您还可以增加 random.uniform(0,1) 的范围，然后相应地调整您的算法以稍微提高性能。

【讨论】：