基于蒙特卡洛的 Pi 计算的 Python 高效矢量化答案

【问题标题】：Python efficient vectorization for Monte Carlo based Pi calculation基于蒙特卡洛的 Pi 计算的 Python 高效矢量化
【发布时间】：2015-03-29 15:52:14
【问题描述】：

为了近似 Pi 的值，请考虑使用随机值填充数组并测试单位圆包含的这种随机方法，

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

非向量化形式的平均速度比向量化形式快 4 倍，例如考虑n = 5000000 和操作系统命令行如下（Python 2.7、Quadcore、8GB RAM、RedHat Linux），

time python pi.py
time python np_pi.py

因此请教如何改进向量化的方法以提高其性能。

【问题讨论】：

请注意，使用 Monte-Carlo 方法计算 pi 是hopelessly inefficient，因此除了演示 Monte Carlo 技术之外，它没有任何用处。如果您想要一种快速计算 pi 的方法，最好使用某种 series approximation。
@BasSwinckels 完全同意，这更像是对 numpy 及其效率的首次探索。

标签： python numpy vectorization montecarlo pi

【解决方案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-array 上调用内置函数 sum 会导致通过遍历数组而不是使用矢量化例程来产生开销。

【讨论】：