「scipy.stat.distributions」的內建概率密度函數是否比用戶提供的慢？

Question

假設我有一個數組：adata=array([0.5, 1.,2.,3.,6.,10.])，並且我想計算該陣列的威布爾分佈的對數似然性，給定參數[5.,1.5]和[5.1,1.6]。我從未想過我需要爲此任務編寫自己的威布爾概率密度函數，因爲它已在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]) 

或者我另起爐竈譜寫新的韋伯PDF功能，如：

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,0100000呼叫wb2LL1(array([[5.,1.5],[5.1,1.6]])[...,newaxis], adata)需要2.156s，而100000呼叫wb2LL(array([[5.,1.5],[5.1,1.6]])[...,newaxis], adata)需要5.219s，內置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()真的慢？如果是這樣，怎麼回事？

Answer 1

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 

所以，它有很多參數處理代碼，這使得它緩慢。