Un coeficiente de Gini ponderado más eficiente en Python

Según https://stackoverflow.com/a/48981834/1840471 , esta es una implementación del coeficiente de Gini ponderado en Python:

import numpy as np def gini(x, weights=None): if weights is None: weights = np.ones_like(x) # Calculate mean absolute deviation in two steps, for weights. count = np.multiply.outer(weights, weights) mad = np.abs(np.subtract.outer(x, x) * count).sum() / count.sum() rmad = mad / np.average(x, weights=weights) # Gini equals half the relative mean absolute deviation. return 0.5 * rmad 

Esto es limpio y funciona bien para matrices de tamaño mediano, pero como se advierte en su sugerencia inicial ( https://stackoverflow.com/a/39513799/1840471 ) es O (n 2 ). En mi computadora, eso significa que se rompe después de ~ 20k filas:

 n = 20000 # Works, 30000 fails. gini(np.random.rand(n), np.random.rand(n)) 

¿Se puede ajustar esto para que funcione con conjuntos de datos más grandes? El mío es ~ 150k filas.

Aquí hay una versión que es mucho más rápida que la que proporcionó anteriormente, y también utiliza una fórmula simplificada para el caso sin peso para obtener resultados aún más rápidos en ese caso.

 def gini(x, w=None): # The rest of the code requires numpy arrays. x = np.asarray(x) if w is not None: w = np.asarray(w) sorted_indices = np.argsort(x) sorted_x = x[sorted_indices] sorted_w = w[sorted_indices] # Force float dtype to avoid overflows cumw = np.cumsum(sorted_w, dtype=float) cumxw = np.cumsum(sorted_x * sorted_w, dtype=float) return (np.sum(cumxw[1:] * cumw[:-1] - cumxw[:-1] * cumw[1:]) / (cumxw[-1] * cumw[-1])) else: sorted_x = np.sort(x) n = len(x) cumx = np.cumsum(sorted_x, dtype=float) # The above formula, with all weights equal to 1 simplifies to: return (n + 1 - 2 * np.sum(cumx) / cumx[-1]) / n 

Aquí hay un código de prueba para verificar que obtengamos (en su mayoría) los mismos resultados:

 >>> x = np.random.rand(1000000) >>> w = np.random.rand(1000000) >>> gini_max_ghenis(x, w) 0.33376310938610521 >>> gini(x, w) 0.33376310938610382 

Pero la velocidad es muy diferente:

 %timeit gini(x, w) 203 ms ± 3.68 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) %timeit gini_max_ghenis(x, w) 55.6 s ± 3.35 s per loop (mean ± std. dev. of 7 runs, 1 loop each) 

Si eliminas las operaciones de pandas de la función, ya es mucho más rápido:

 %timeit gini_max_ghenis_no_pandas_ops(x, w) 1.62 s ± 75 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) 

Si desea obtener la última caída de rendimiento, puede usar numba o cython, pero eso solo ganaría un pequeño porcentaje porque la mayor parte del tiempo se dedica a la clasificación.

 %timeit ind = np.argsort(x); sx = x[ind]; sw = w[ind] 180 ms ± 4.82 ms per loop (mean ± std. dev. of 7 runs, 10 loops each) 

edit : gini_max_ghenis es el código utilizado en la respuesta de Max Ghenis

Adaptando la función StatsGini R desde aquí :

 import numpy as np import pandas as pd def gini(x, w=None): # Array indexing requires reset indexes. x = pd.Series(x).reset_index(drop=True) if w is None: w = np.ones_like(x) w = pd.Series(w).reset_index(drop=True) n = x.size wxsum = sum(w * x) wsum = sum(w) sxw = np.argsort(x) sx = x[sxw] * w[sxw] sw = w[sxw] pxi = np.cumsum(sx) / wxsum pci = np.cumsum(sw) / wsum g = 0.0 for i in np.arange(1, n): g = g + pxi.iloc[i] * pci.iloc[i - 1] - pci.iloc[i] * pxi.iloc[i - 1] return g 

Esto funciona para vectores grandes, al menos hasta 10M de filas:

 n = 1e7 gini(np.random.rand(n), np.random.rand(n)) # Takes ~15s. 

También produce el mismo resultado que la función proporcionada en la pregunta, por ejemplo, dando 0.2553 para este ejemplo:

 gini(np.array([3, 1, 6, 2, 1]), np.array([4, 2, 2, 10, 1]))