Tabla de contingencia de Python

Estoy generando muchas, muchas tablas de contingencia como parte de un proyecto que estoy escribiendo.

El flujo de trabajo es:

  • Tome una matriz de datos grande con filas continuas (flotantes) y conviértalas a valores enteros discretos al agruparlas (de modo que la fila resultante tenga valores de 0 a 9, por ejemplo)
  • Corte dos filas en vectores X e Y y genere una tabla de contingencia a partir de ellas, de modo que tenga la distribución de frecuencia bidimensional.
  • Por ejemplo, tendría una matriz de 10 x 10, contando el número de (xi, yi) que ocurren
  • Usa la tabla de contingencia para hacer un poco de teoría de la información matemática.

Inicialmente, escribí esto como:

def make_table(x, y, num_bins): ctable = np.zeros((num_bins, num_bins), dtype=np.dtype(int)) for xn, yn in zip(x, y): ctable[xn, yn] += 1 return ctable 

Esto funciona bien, pero es tan lento que consume casi el 90% del tiempo de ejecución de todo el proyecto.

La optimización de Python solo más rápida que he podido encontrar es la siguiente:

 def make_table(x, y, num_bins): ctable = np.zeros(num_bins ** 2, dtype=np.dtype(int)) reindex = np.dot(np.stack((x, y)).transpose(), np.array([num_bins, 1])) idx, count = np.unique(reindex, return_counts=True) for i, c in zip(idx, count): ctable[i] = c return ctable.reshape((num_bins, num_bins)) 

Eso es (de alguna manera) mucho más rápido, pero sigue siendo bastante caro para algo que no parece ser un cuello de botella. ¿Hay alguna forma eficiente de hacer esto que no esté viendo o debo rendirme y hacer esto en cython?

Además, aquí hay una función de evaluación comparativa.

 def timetable(func): size = 5000 bins = 10 repeat = 1000 start = time.time() for i in range(repeat): x = np.random.randint(0, bins, size=size) y = np.random.randint(0, bins, size=size) func(x, y, bins) end = time.time() print("Func {na}: {ti} Ms".format(na=func.__name__, ti=(end - start))) 

El truco inteligente para representar los elementos de np.stack((x, y)) como enteros se puede hacer más rápido:

 In [92]: %timeit np.dot(np.stack((x, y)).transpose(), np.array([bins, 1])) 109 µs ± 6.55 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each) In [94]: %timeit bins*x + y 12.1 µs ± 260 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each) 

Además, la última parte de su segunda solución se puede simplificar un poco, simplemente considerando

 np.unique(bins * x + y, return_counts=True)[1].reshape((bins, bins)) 

Además, dado que estamos tratando con enteros no negativos igualmente espaciados, np.bincount superará a np.unique ; Con eso, lo anterior se reduce a

 np.bincount(bins * x + y).reshape((bins, bins)) 

Con todo, esto proporciona bastante rendimiento sobre lo que está haciendo actualmente:

 In [78]: %timeit make_table(x, y, bins) # Your first solution 3.86 ms ± 159 µs per loop (mean ± std. dev. of 7 runs, 100 loops each) In [79]: %timeit make_table2(x, y, bins) # Your second solution 443 µs ± 23.7 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each) In [101]: %timeit np.unique(bins * x + y, return_counts=True)[1].reshape((bins, bins)) 307 µs ± 25 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each) In [118]: %timeit np.bincount(bins * x + y).reshape((10, 10)) 30.3 µs ± 3.44 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each) 

También es posible que desee conocer np.histogramdd que se encarga de redondear y agrupar a la vez, aunque es probable que sea más lento que redondear y usar np.bincount .