Python 中更有效的加权基尼系数
More efficient weighted Gini coefficient in Python
根据,这是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
这很干净,适用于中型阵列,但正如其最初建议 () 中所警告的那样,它的复杂度为 O(n2)。在我的电脑上,这意味着它在 ~20k 行后中断:
n = 20000 # Works, 30000 fails.
gini(np.random.rand(n), np.random.rand(n))
是否可以调整它以适用于更大的数据集?我的大约有 15 万行。
调整 StatsGini
来自 here 的 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
这适用于大向量,至少高达 1000 万行:
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]))
这是一个比您上面提供的版本快得多的版本,并且还对没有权重的情况使用了简化的公式,以便在这种情况下获得更快的结果。
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_max_ghenis(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_max_ghenis(x, w)
55.6 s ± 3.35 s per loop (mean ± std. dev. of 7 runs, 1 loop each)
如果从函数中删除 pandas 操作,它已经快得多:
%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)
如果你想获得最后一滴性能,你可以使用 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)
编辑:gini_max_ghenis 是 Max Ghenis 的回答中使用的代码
根据
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
这很干净,适用于中型阵列,但正如其最初建议 (
n = 20000 # Works, 30000 fails.
gini(np.random.rand(n), np.random.rand(n))
是否可以调整它以适用于更大的数据集?我的大约有 15 万行。
调整 StatsGini
来自 here 的 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
这适用于大向量,至少高达 1000 万行:
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]))
这是一个比您上面提供的版本快得多的版本,并且还对没有权重的情况使用了简化的公式,以便在这种情况下获得更快的结果。
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_max_ghenis(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_max_ghenis(x, w)
55.6 s ± 3.35 s per loop (mean ± std. dev. of 7 runs, 1 loop each)
如果从函数中删除 pandas 操作,它已经快得多:
%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)
如果你想获得最后一滴性能,你可以使用 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)
编辑:gini_max_ghenis 是 Max Ghenis 的回答中使用的代码