title: Is the build-in probability density functions of `scipy.stat.distributions` slower than a user provided one?
body: <p>Suppose I have an array: <code>adata=array([0.5, 1.,2.,3.,6.,10.])</code> and I want to calculate log likelihood of Weibull distribution of this array, given the parameters <code>[5.,1.5]</code> and <code>[5.1,1.6]</code>. I have never thought I need to write my own Weibull probability density functions for this task, as it is already provided in <code>scipy.stat.distributions</code>. So, this ought to do it:</p>
  
  <pre><code>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)
  </code></pre>
  
  <p>And: </p>
  
  <pre><code>&gt;&gt;&gt; wb2LL(array([[5.,1.5],[5.1,1.6]]).T[...,newaxis], adata)
  array([-14.43743911, -14.68835298])
  </code></pre>
  
  <p>Or I reinvent the wheel and write a new Weibull pdf function, such as:</p>
  
  <pre><code>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)
  </code></pre>
  
  <p>And: </p>
  
  <pre><code>&gt;&gt;&gt; wb2LL1(array([[5.,1.5],[5.1,1.6]]).T[...,newaxis], adata)
  array([-14.43743911, -14.68835298])
  </code></pre>
  
  <p>Admittedly, I always take it for granted that if something is already in <code>scipy</code>, it should be very well optimized and re-inventing the wheel is seldom going to make it faster. But here comes the surprise: if I <code>timeit</code>, 100000 calls of <code>wb2LL1(array([[5.,1.5],[5.1,1.6]])[...,newaxis], adata)</code> takes 2.156s while 100000 calls of <code>wb2LL(array([[5.,1.5],[5.1,1.6]])[...,newaxis], adata)</code> takes 5.219s, the build-in <code>stats.weibull_min.pdf()</code> is more than twice slower.</p>
  
  <p>Checking the source code <code>python_path/sitepackage/scipy/stat/distributions.py</code> didn't provides an easy answer, at least for me. If anything, from it I would expect <code>stats.weibull_min.pdf()</code> to be almost as fast as <code>wbp()</code>.</p>
  
  <p>Relevant source code: line 2999-3033:</p>
  
  <pre><code>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 &gt; 0``, ``c &gt; 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')
  </code></pre>
  
  <p>So, the question is: is <code>stats.weibull_min.pdf()</code> really slower? If so, how come?</p>