Более эффективный весовой коэффициент Джини в Python

Question

Согласно /questions/19458335/vesovoj-koeffitsient-dzhini-v-python/19458343#19458343, это реализация взвешенного коэффициента Джини в 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

Это чисто и хорошо работает для массивов среднего размера, но, как и предупреждено в первоначальном предложении ( /questions/37585065/raschet-koeffitsienta-dzhini-v-python-numpy/37585070#37585070), это O (n²). На моем компьютере это означает, что он ломается после ~20k строк:

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

Можно ли это настроить для работы с большими наборами данных? Мой ~150 тыс рядов.

Answer 1

Вот версия, которая намного быстрее, чем та, которую вы предоставили выше, а также использует упрощенную формулу для случая без веса, чтобы получить еще более быстрые результаты в этом случае.

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

Вот некоторый тестовый код, чтобы проверить, что мы получаем (в основном) одинаковые результаты:

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

Но скорость совсем другая

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

%timeit gini_slow(x, w)
55.6 s ± 3.35 s per loop (mean ± std. dev. of 7 runs, 1 loop each)

Если вы удалите функции pandas из функции, это уже намного быстрее:

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

Если вы хотите получить последнее снижение производительности, вы можете использовать numba или cython, но это принесет лишь несколько процентов, потому что большая часть времени уходит на сортировку.

%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)

Answer 2

Адаптация StatsGini Функция R отсюда:

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

Это работает для больших векторов, по крайней мере, до 10 миллионов строк:

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

Он также дает тот же результат, что и функция, представленная в вопросе, например, дает 0,2553 для этого примера:

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