Comparing scipy.stats.t.sf vs GSL using Cython
<p>I would like to calculate the p values of a large 2D numpy t values array. However, this takes long time and I would like to improve its speed. I tried using GSL. 
Although a single gsl_cdf_tdist_P is much much faster than scipy.stats.t.sf, when iterating over the ndarray, the process is very slow. I would like help to improve this. 
See the code below.</p>

<p>GSL_Test.pyx</p>

<pre><code>import cython
cimport cython

import numpy
cimport numpy
from cython_gsl cimport *

DTYPE = numpy.float32
ctypedef numpy.float32_t DTYPE_t

cdef get_gsl_p(double t, double nu):
    return (1 - gsl_cdf_tdist_P(t, nu)) * 2

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
cdef get_gsl_p_for_2D_matrix(numpy.ndarray[DTYPE_t, ndim=2] t_matrix, int n):

    cdef unsigned int rows = t_matrix.shape[0]
    cdef numpy.ndarray[DTYPE_t, ndim=2] out = numpy.zeros((rows, rows), dtype='float32')
    cdef unsigned int row, col

    for row in range(rows):
        for col in range(rows):
            out[row, col] = get_gsl_p(t_matrix[row, col], n-2)
    return out

def get_gsl_p_for_2D_matrix_def(numpy.ndarray[DTYPE_t, ndim=2] t_matrix, int n):
    return get_gsl_p_for_2D_matrix(t_matrix, n)
</code></pre>

<p>ipython</p>

<pre><code>import GSL_Test
import numpy
import scipy.stats
a = numpy.random.rand(3544, 3544).astype('float32')
%timeit -n 1 GSL_Test.get_gsl_p_for_2D_matrix(a, 25)
1 loop, best of 3: 7.87 s per loop
%timeit -n 1 scipy.stats.t.sf(a, 25)*2
1 loop, best of 3: 4.66 s per loop
</code></pre>

<p><strong>UPDATE:</strong> Adding cdef declarations I was able to  reduce the computational time but not lower than scipy still. I modified the code to have the cdef declarations. </p>

<pre><code>%timeit -n 1 GSL_Test.get_gsl_p_for_2D_matrix_def(a, 25)
1 loop, best of 3: 6.73 s per loop
</code></pre>

<p>You can get some small gain in raw performance by using a raw special function instead of <code>stats.t.sf</code>. Looking at the source, you find  (<a href="https://github.com/scipy/scipy/blob/master/scipy/stats/_continuous_distns.py#L3849" rel="nofollow">https://github.com/scipy/scipy/blob/master/scipy/stats/_continuous_distns.py#L3849</a>)</p>

<pre><code>def _sf(self, x, df):
      return sc.stdtr(df, -x)
</code></pre>

<p>So that you can use <code>stdtr</code> directly:</p>

<pre><code>np.random.seed(1234)
x = np.random.random((3740, 374))
t1 = stats.t.sf(x, 25)
t2 = stdtr(25, -x)

1 loop, best of 3: 653 ms per loop
1 loop, best of 3: 562 ms per loop
</code></pre>

<p>If you do reach out for cython, the typed memoryview syntax often gives you faster code than the old ndarray syntax:</p>

<pre><code>from scipy.special.cython_special cimport stdtr
from numpy cimport npy_intp
import numpy as np

def tsf(double [:, ::1] x, int df=25):
    cdef double[:, ::1] out = np.empty_like(x)
    cdef npy_intp i, j
    cdef double tmp, xx

    for i in range(x.shape[0]):
        for j in range(x.shape[1]):
            xx = x[i, j]
            out[i, j] = stdtr(df, -xx)

    return np.asarray(out)
</code></pre>

<p>Here I'm also using the <code>cython_special</code> interface, which is only avaialble in the dev version of scipy (<a href="http://scipy.github.io/devdocs/special.cython_special.html#module-scipy.special.cython_special" rel="nofollow">http://scipy.github.io/devdocs/special.cython_special.html#module-scipy.special.cython_special</a>), but you can use GSL if you want.</p>

<p>Finally, if you suspect a bottleneck in iterations, don't forget to inspect the output of <code>cython -a</code> to see if there's some python overhead in the hot loops.</p>

python numpy scipy cython gsl