`scipy.stat.distributions`内置的概率密度函数比用户提供的函数速度慢吗？

Question

`scipy.stat.distributions`内置的概率密度函数比用户提供的函数速度慢吗？

5

假设我有一个数组：adata=array([0.5, 1.,2.,3.,6.,10.])，我想计算这个数组的Weibull分布的对数似然，给定参数[5.,1.5]和[5.1,1.6]。我从未想过需要为此编写自己的Weibull概率密度函数，因为它已经在scipy.stat.distributions中提供了。所以，这应该可以解决问题：

from scipy import stats
from numpy import *
adata=array([0.5, 1.,2.,3.,6.,10.])
def wb2LL(p, x): #log-likelihood of 2 parameter Weibull distribution
    return sum(log(stats.weibull_min.pdf(x, p[1], 0., p[0])), axis=1)

而且：

>>> wb2LL(array([[5.,1.5],[5.1,1.6]]).T[...,newaxis], adata)
array([-14.43743911, -14.68835298])

要么我重新发明轮子并编写一个新的威布尔概率密度函数，例如：

wbp=lambda p, x: p[1]/p[0]*((x/p[0])**(p[1]-1))*exp(-((x/p[0])**p[1]))
def wb2LL1(p, x): #log-likelihood of 2 paramter Weibull
    return sum(log(wbp(p,x)), axis=1)

并且：

>>> wb2LL1(array([[5.,1.5],[5.1,1.6]]).T[...,newaxis], adata)
array([-14.43743911, -14.68835298])

承认，我总是把如果某个东西已经在scipy中了看作理所当然，因为它应该被很好地优化过了，重新发明轮子很少会让它更快。但是出现了一个惊喜：如果我使用timeit，100000次调用wb2LL1(array([[5.,1.5],[5.1,1.6]])[...,newaxis], adata)需要2.156秒，而100000次调用wb2LL(array([[5.,1.5],[5.1,1.6]])[...,newaxis], adata)需要5.219秒，内置的stats.weibull_min.pdf()要慢两倍以上。

检查源代码python_path/sitepackage/scipy/stat/distributions.py并没有提供一个简单的答案，至少对我来说不是。如果有什么，我希望从中得到的是stats.weibull_min.pdf()几乎和wbp()一样快。

相关源代码：2999-3033行。

class frechet_r_gen(rv_continuous):
    """A Frechet right (or Weibull minimum) continuous random variable.

    %(before_notes)s

    See Also
    --------
    weibull_min : The same distribution as `frechet_r`.
    frechet_l, weibull_max

    Notes
    -----
    The probability density function for `frechet_r` is::

        frechet_r.pdf(x, c) = c * x**(c-1) * exp(-x**c)

    for ``x > 0``, ``c > 0``.

    %(example)s

    """
    def _pdf(self, x, c):
        return c*pow(x,c-1)*exp(-pow(x,c))
    def _logpdf(self, x, c):
        return log(c) + (c-1)*log(x) - pow(x,c)
    def _cdf(self, x, c):
        return -expm1(-pow(x,c))
    def _ppf(self, q, c):
        return pow(-log1p(-q),1.0/c)
    def _munp(self, n, c):
        return special.gamma(1.0+n*1.0/c)
    def _entropy(self, c):
        return -_EULER / c - log(c) + _EULER + 1
frechet_r = frechet_r_gen(a=0.0, name='frechet_r', shapes='c')
weibull_min = frechet_r_gen(a=0.0, name='weibull_min', shapes='c')

所以问题是： stats.weibull_min.pdf() 真的比较慢吗？如果是，为什么？

- CT Zhu

1

我可能对那个手写的函数存在一些数值稳定性方面的担忧。为了弥补这一点，可能需要在内置函数中花费一些时间。 - user2357112

1

我有同感，有点对scipy/python失去了信心。这个库函数并不像我预期的那样先进。 - colinfang

1

在 scipy 中，%timeit ss.weibull_min.pdf(0.1*np.arange(1,11),5.,scale=1) 报告了 1000 loops, best of 3: 595 µs per loop 的结果。而在 R 中，benchmark(pweibull((1:10)*0.1, 5 ,1), replications=10^6) 则报告了总计 22.83 秒的运行时间。作为一种编程语言，python 比 R 更加出色是有争议的。但很多时候，我也希望 scipy 库的函数能够更快一些。 - CT Zhu

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- HYRY · Accepted Answer

pdf()方法是在rv_continuous类中定义的，它调用了frechet_r_gen._pdf()。 pdf()的代码如下：

def pdf(self,x,*args,**kwds):
    loc,scale=map(kwds.get,['loc','scale'])
    args, loc, scale = self._fix_loc_scale(args, loc, scale)
    x,loc,scale = map(asarray,(x,loc,scale))
    args = tuple(map(asarray,args))
    x = asarray((x-loc)*1.0/scale)
    cond0 = self._argcheck(*args) & (scale > 0)
    cond1 = (scale > 0) & (x >= self.a) & (x <= self.b)
    cond = cond0 & cond1
    output = zeros(shape(cond),'d')
    putmask(output,(1-cond0)+np.isnan(x),self.badvalue)
    if any(cond):
        goodargs = argsreduce(cond, *((x,)+args+(scale,)))
        scale, goodargs = goodargs[-1], goodargs[:-1]
        place(output,cond,self._pdf(*goodargs) / scale)
    if output.ndim == 0:
        return output[()]
    return output

因此，它有许多参数处理代码，这使得它变慢。