我有一个大小为n的样本。
我想要计算每个i:1 <= i <= n在numpy中sample[:i]
的中位数。例如,我可以为每个i计算平均值:
cummean = np.cumsum(sample) / np.arange(1, n + 1)
我能否做类似于中位数的事情,而不需要循环和理解?
我有一个大小为n的样本。
我想要计算每个i:1 <= i <= n在numpy中sample[:i]
的中位数。例如,我可以为每个i计算平均值:
cummean = np.cumsum(sample) / np.arange(1, n + 1)
我能否做类似于中位数的事情,而不需要循环和理解?
了解Python有一个heapq
模块,可以让您为可迭代对象保持运行的“最小值”,因此我搜索了heapq
和median
,并找到了各种各样的关于steaming medium
的内容。其中这个网址:
http://www.ardendertat.com/2011/11/03/programming-interview-questions-13-median-of-integer-stream/
有一个class streamMedian
,它维护两个heapq
,一个存储值的下半部分,另一个存储上半部分。中位数是其中一个的“顶部”或来自两个堆的值的平均值。该类具有insert
方法和getMedian
方法。大部分工作都在insert
方法中。
我将其复制到Ipython会话中,并定义:
def cummedian_stream(b):
S=streamMedian()
ret = []
for item in b:
S.insert(item)
ret.append(S.getMedian())
return np.array(ret)
测试:
In [155]: a = np.random.randint(0,100,(5000))
In [156]: amed = cummedian_stream(a)
In [157]: np.allclose(cummedian_sorted(a), amed)
Out[157]: True
In [158]: timeit cummedian_sorted(a)
1 loop, best of 3: 781 ms per loop
In [159]: timeit cummedian_stream(a)
10 loops, best of 3: 39.6 ms per loop
heapq
流式方法更快。
@Uriel
给出的列表推导式相对较慢。但如果我将 statistics.median
替换为 np.median
,它比 @Divakar
的排序解决方案更快:def fastloop(a):
return np.array([np.median(a[:i+1]) for i in range(len(a))])
In [161]: timeit fastloop(a)
1 loop, best of 3: 360 ms per loop
而且 @Paul Panzer
的分区方法也不错,但与流式处理相比仍然较慢。
In [165]: timeit cummedian_partition(a)
1 loop, best of 3: 391 ms per loop
如果需要,我可以复制StreamMedian
类到这个答案中。
2D
数组。然后,我们会将上三角区域填充为一个大数字,这样当我们对每行排序时,基本上会将所有元素排序到对角线元素,并模拟累计窗口。然后,按照选择中间一个或两个中间的平均值(对于偶数个元素)的median
的定义,我们将获得第一个位置的元素:(0,0)
,然后对于第二行: (1,0) & (1,1)
的均值,对于第三行:(2,1)
,对于第四行: (3,1) & (3,2)
的均值等等。因此,我们将从排序后的数组中提取出这些元素,从而获得中位数。
因此,实现如下 -
def cummedian_sorted(a):
n = a.size
maxn = a.max()+1
a_tiled_sorted = np.tile(a,n).reshape(-1,n)
mask = np.triu(np.ones((n,n),dtype=bool),1)
a_tiled_sorted[mask] = maxn
a_tiled_sorted.sort(1)
all_rows = a_tiled_sorted[np.arange(n), np.arange(n)//2].astype(float)
idx = np.arange(1,n,2)
even_rows = a_tiled_sorted[idx, np.arange(1,1+(n//2))]
all_rows[idx] += even_rows
all_rows[1::2] /= 2.0
return all_rows
运行时测试
方法 -
# Loopy solution from @Uriel's soln
def cummedian_loopy(arr):
return [median(a[:i]) for i in range(1,len(a)+1)]
# Nan-fill based solution from @Nickil Maveli's soln
def cummedian_nanfill(arr):
a = np.tril(arr).astype(float)
a[np.triu_indices(a.shape[0], k=1)] = np.nan
return np.nanmedian(a, axis=1)
时序 -
第一组:
In [43]: a = np.random.randint(0,100,(100))
In [44]: print np.allclose(cummedian_loopy(a), cummedian_sorted(a))
...: print np.allclose(cummedian_loopy(a), cummedian_nanfill(a))
...:
True
True
In [45]: %timeit cummedian_loopy(a)
...: %timeit cummedian_nanfill(a)
...: %timeit cummedian_sorted(a)
...:
1000 loops, best of 3: 856 µs per loop
1000 loops, best of 3: 778 µs per loop
10000 loops, best of 3: 200 µs per loop
第二组:
In [46]: a = np.random.randint(0,100,(1000))
In [47]: print np.allclose(cummedian_loopy(a), cummedian_sorted(a))
...: print np.allclose(cummedian_loopy(a), cummedian_nanfill(a))
...:
True
True
In [48]: %timeit cummedian_loopy(a)
...: %timeit cummedian_nanfill(a)
...: %timeit cummedian_sorted(a)
...:
10 loops, best of 3: 118 ms per loop
10 loops, best of 3: 47.6 ms per loop
100 loops, best of 3: 18.8 ms per loop
第三部分:
In [49]: a = np.random.randint(0,100,(5000))
In [50]: print np.allclose(cummedian_loopy(a), cummedian_sorted(a))
...: print np.allclose(cummedian_loopy(a), cummedian_nanfill(a))
True
True
In [54]: %timeit cummedian_loopy(a)
...: %timeit cummedian_nanfill(a)
...: %timeit cummedian_sorted(a)
...:
1 loops, best of 3: 3.36 s per loop
1 loops, best of 3: 583 ms per loop
1 loops, best of 3: 521 ms per loop
loopy
中获得了更快的时间,但我使用的是 np.median
,而不是统计学中的那个。 - hpaulj使用statistics.median
和累积列表推导(注意奇数索引包含偶数长度列表的中位数 - 其中中位数是两个中位数元素的平均值,因此通常会得到一个小数而不是整数):
>>> from statistics import median
>>> arr = [1, 3, 4, 2, 5, 3, 6]
>>> cum_median = [median(arr[:i+1]) for i in range(len(arr)-1)]
>>> cum_median
[1, 2.0, 3, 2.5, 3, 3.0]
np.median
在这个角色中比 statistics.median
快得多。 - hpaulj如果将值列表arr
视为概率质量函数,则可以使用np.cumsum(arr)
获取累积分布函数从而得到一个近似解。按定义,中位数仅适用于一半的概率。
arr[np.searchsorted(np.cumsum(arr), np.cumsum(arr)/2)]
arr
已经排序,因此您可以简单地查看每个i
的索引i//2
。 - hilberts_drinking_problem有没有晚点进入的空间?
def cummedian_partition(a):
n = len(a)
assert n%4 == 0 # for simplicity
mn = a.min() - 1
mx = a.max() + 1
h = n//2
N = n + h//2
work = np.empty((h, N), a.dtype)
work[:, :n] = a
work[:, n] = 2*mn - a[0]
i, j = np.tril_indices(h, -1)
work[i, n-1-j] = (2*mn - a[1:h+1])[j]
k, l = np.ogrid[:h, :h//2 - 1]
work[:, n+1:] = np.where(k > 2*l+1, mx, 2 * mn - mx)
out = np.partition(work, (N-n//2-1, N-n//2, h//2-1, h//2), axis=-1)
out = np.r_[2*mn-out[:, h//2: h//2-2:-1], out[::-1, N-n//2-1:N-n//2+1]]
out[::2, 0] = out[::2, 1]
return np.mean(out, axis=-1)
np.partition
不支持每行拆分点,因此需要进行一些技巧调整。所需的综合复杂度和内存是二次的。for j in (100, 1000, 5000):
a = np.random.randint(0, 100, (j,))
print('size', j)
print('correct', np.allclose(cummedian_partition(a), cummedian_sorted(a)))
print('Divakar', timeit(lambda: cummedian_sorted(a), number=10))
print('PP', timeit(lambda: cummedian_partition(a), number=10))
# size 100
# correct True
# Divakar 0.0022412699763663113
# PP 0.002393342030700296
# size 1000
# correct True
# Divakar 0.20881508802995086
# PP 0.10222102201078087
# size 5000
# correct True
# Divakar 6.158387024013791
# PP 3.437395485001616
numpy
也无法使一个O(n^2)
(或更糟糕的)算法比一个 -- 是什么 --O(n log n)
或者更好的算法更快。 - Paul Panzer