Python高效向量化的蒙特卡罗方法计算Pi值

Question

对于近似计算圆周率的方法，考虑使用随机方法，通过生成随机数并测试其是否在单位圆内来进行计算。

import random as rd
import numpy as np

def r(_): return rd.random()

def np_pi(n):
    v_r = np.vectorize(r)
    x = v_r(np.zeros(n))
    y = v_r(np.zeros(n))

    return sum (x*x + y*y <= 1) * 4. / n

请注意，随机数生成依赖于Python标准库；但是考虑使用numpy进行随机数生成。

def np_pi(n):
   x = np.random.random(n)
   y = np.random.random(n)

    return sum (x*x + y*y <= 1) * 4. / n

考虑非向量化方法，

import random as rd

def dart_board():
    x,y = rd.random(), rd.random()
    return (x*x + y*y <= 1)

def pi(n):
    s = sum([dart_board() for _ in range(n)])
    return s * 4. / n

非矢量化形式的平均速度比矢量化形式快四倍，例如考虑 n = 5000000 和操作系统命令行如下所示（Python 2.7，四核，8GB RAM，RedHat Linux）。

time python pi.py
time python np_pi.py

因此，如何改进矢量化方法以提高其性能是一个需要探讨的问题。

- elm

---

请注意，使用蒙特卡罗方法来计算π是完全低效的（https://dev59.com/D3XYa4cB1Zd3GeqP6nqk#18141358），因此它除了演示蒙特卡罗技术外没有其他用途。如果您想要一种快速计算π的方法，最好使用某种级数逼近方法（http://mathworld.wolfram.com/PiFormulas.html）。 - Bas Swinckels

@BasSwinckels 完全同意，这更多是对numpy及其效率的初步探索。 - elm

Answer 1

你正在调用Python内置的sum函数，而不是NumPy的向量化方法sum：

import numpy as np
import random as rd

def np_pi(n):
    x = np.random.random(n)
    y = np.random.random(n)

    return (x*x + y*y <= 1).sum()

def dart_board():
    x,y = rd.random(), rd.random()
    return (x*x + y*y <= 1)

def pi(n):
    s = sum([dart_board() for _ in range(n)])

现在的计时结果有很大不同：

In [12]: %timeit np_pi(10000)
1000 loops, best of 3: 250 us per loop

In [13]: %timeit pi(10000)
100 loops, best of 3: 3.54 ms per loop

我猜测，在numpy数组上调用内置的sum会导致开销，因为它会通过迭代数组来计算，而不是使用矢量化算法。

- Oliver W.