使用蒙特卡罗方法求解圆周率的小数位数

Question

我尝试过很多使用蒙特卡罗方法来寻找π的算法。 其中一种解决方案（使用Python）如下：

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

可悲的是，即使使用了十亿的精度，结果也非常不准确（3.141...）。

这种方法的最大精度是多少呢？ 我选择蒙特卡罗法的原因是它很容易被分成几个并行部分进行计算。 是否有另一种能够轻松拆分成小块进行计算得出π值的算法呢？

- Jon Romero

Answer 1

这是蒙特卡罗方法的一个经典示例。但如果你试图将π的计算分解成并行部分，为什么不使用无穷级数，让每个核心负责一定范围的计算，然后随着计算进行逐步累加结果呢？

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

- Nosredna

---

那是我的第一种方法。但我想尝试一下蒙特卡罗方法，因为它可以应用于许多领域。 - Jon Romero

当很难找到公式时，请使用蒙特卡罗方法。当很容易找到公式时，请使用公式。 - Nosredna

Answer 2

你的分数误差按照 sqrt(N)/N = 1/sqrt(N) 来计算，因此这是一种获得精确估计的非常低效的方法。这个限制是由测量的统计特性决定的，无法被打破。

你应该能够获得大约 floor(log_10(N))/2-1 位好的精度，对于 N 次抛掷。也许为了安全起见，可以再减去 -2。

即使如此，它也假设你使用的是真正的随机数生成器或足够好的伪随机数生成器。

- dmckee --- ex-moderator kitten

Answer 3

使用准随机数生成器（http://www.nag.co.uk/IndustryArticles/introduction_to_quasi_random_numbers.pdf）代替标准伪随机数生成器。准随机数更均匀地覆盖积分区域（您正在进行蒙特卡罗积分），比伪随机数更好地收敛。

- quant_dev

---

我的天真猜测是，虽然这个方法会更快地收敛，但估计置信区间可能更困难。你知道有关此事的文献吗？ - dmckee --- ex-moderator kitten

不过，有一个C库http://www.feynarts.de/cuba/实现了包括置信区间在内的MC积分（它返回绝对误差估计和Chi^2概率，即估计错误的概率）。您可以下载代码并查看实现，或给作者发电子邮件询问他用来编写代码的文献。 - quant_dev

啊！作者在文中提供了一篇论文链接。该论文发表在《计算物理学通讯》上（并且在arXiv上也有http://arxiv.org/abs/hep-ph/0404043）。令我不安的是，他和我从事同样的工作，而这是我第一次听说这个。哎呀！ - dmckee --- ex-moderator kitten

我妻子在她的研究中广泛使用这个库。根据她的经验，在低维情况下，随机蒙特卡罗积分方法的性能被确定性方法（古巴库中的"Cuhre"例程）超越，即使对于带有不连续性的“困难”被积函数也是如此。 - quant_dev