python - 如何加速 Poisson pmf 函数？

Question

我的用例是在小于 10 的所有点上评估 Poisson pmf，我会使用不同的 lambda 多次调用这样的函数。无法提前知道 lambda，因此我无法矢量化 lambda。

我从某个地方听说了一个 secret 技巧，那就是使用_pmf。这样做的缺点是什么？但是，它还是有点慢，有没有什么方法可以改进它而不用从头开始用 C 语言重写 pmf？

%timeit scipy.stats.poisson.pmf(np.arange(0,10),3.3)
%timeit scipy.stats.poisson._pmf(np.arange(0,10),3.3)
a = np.arange(0,10)
%timeit scipy.stats.poisson._pmf(a,3.3)

10000 loops, best of 3: 94.5 µs per loop
100000 loops, best of 3: 15.2 µs per loop
100000 loops, best of 3: 13.7 µs per loop

更新

好吧，只是我懒得用 cython 写。我曾期望有一个更快的解决方案可以针对连续的 x 按顺序(迭代)评估所有离散分布。例如。 P(X=3) = P(X=2) * lambda/3 如果 X ~ Pois(lambda)

相关:Is the build-in probability density functions of `scipy.stat.distributions` slower than a user provided one?

我现在不太相信 Scipy 和 Python。库功能没有我预期的那么先进。

Answer 1

大多数 scipy.stats 发行版都支持向量化评估:

>>> poisson.pmf(1, [5, 6, 7, 8])
array([ 0.03368973,  0.01487251,  0.00638317,  0.0026837 ])

这可能不够快，也可能不够快，但您可以尝试将 pmf 调用移出循环。

关于 pmf 和 _pmf 的区别:真正的工作是在带下划线的函数中完成的 (_pmf, _cdf 等)，而公共(public)函数(pmf、cdf)确保只有有效参数才能进入 _pmf( 的输出如果参数无效，则不能保证 >_pmf 有意义，因此使用风险自负)。

>>> poisson.pmf(1, -1)
nan
>>> poisson._pmf(1, -1)
/home/br/virtualenvs/scipy-dev/local/lib/python2.7/site-packages/scipy/stats/_discrete_distns.py:432: RuntimeWarning: invalid value encountered in log
  Pk = k*log(mu)-gamln(k+1) - mu
nan

更多详情:https://github.com/scipy/scipy/blob/master/scipy/stats/_distn_infrastructure.py#L2721