有效使用 numpy.random.choice 与重复数字和备选方案
Efficient use of numpy.random.choice with repeated numbers and alternatives
我需要生成一个包含重复元素的大数组，我的代码是：
np.repeat(xrange(x,y), data)
但是，数据是一个类型为 float64 的 numpy 数组（但它表示整数，那里没有 2.1），我收到错误
TypeError: Cannot cast array data from dtype('float64') to dtype('int64') according to the rule 'safe'
示例：
In [35]: x
Out[35]: 26
In [36]: y
Out[36]: 50
In [37]: data
Out[37]:
array([ 3269., 106., 5533., 317., 1512., 208., 502., 919.,
406., 421., 1690., 2236., 705., 505., 230., 213.,
307., 1628., 4389., 1491., 355., 103., 854., 424.])
In [38]: np.repeat(xrange(x,y), data)
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-38-105860821359> in <module>()
----> 1 np.repeat(xrange(x,y), data)
/home/pcadmin/anaconda2/lib/python2.7/site-packages/numpy /core/fromnumeric.pyc in repeat(a, repeats, axis)
394 repeat = a.repeat
395 except AttributeError:
--> 396 return _wrapit(a, 'repeat', repeats, axis)
397 return repeat(repeats, axis)
398
/home/pcadmin/anaconda2/lib/python2.7/site-packages/numpy /core/fromnumeric.pyc in _wrapit(obj, method, *args, **kwds)
46 except AttributeError:
47 wrap = None
---> 48 result = getattr(asarray(obj), method)(*args, **kwds)
49 if wrap:
50 if not isinstance(result, mu.ndarray):
TypeError: Cannot cast array data from dtype('float64') to dtype('int64') according to the rule 'safe'
我把代码改成
解决了
np.repeat(xrange(x,y), data.astype('int64'))
但是，现在这是我代码中最昂贵的行之一！！还有其他选择吗？
顺便说一句，我在里面用这个
np.random.choice(np.repeat(xrange(x,y), data.astype('int64')), z)
为了在不替换的情况下获得样本，其中 x 和 y 之间的整数的大小为 z，数据中给出了每个整数的数量。我想这也是最好的方法吧？
问题陈述
这个问题很有趣！只是为了让读者对这个问题有一个概念，而不是进入次要的数据转换问题，我们有一个值的范围，比方说 a = np.arange(5)，即
a = np.array([0,1,2,3,4])
现在，假设我们有另一个数组，其中列出了 a 中每个 5 数字的重复次数。所以，让这些成为：
reps = np.array([2,4,6,2,2])
接下来，我们将进行这些重复：
In [32]: rep_nums = np.repeat(a,reps)
In [33]: rep_nums
Out[33]: array([0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 4, 4])
最后，我们希望使用 np.random.choice() 从那些重复的数字中选择 z 个元素，而不用替换。
假设 z = 7 选择 7 个元素，那么对于 np.random.choice()，我们将有：
In [34]: np.random.choice(rep_nums,7,replace=False)
Out[34]: array([2, 4, 0, 2, 4, 1, 2])
现在，这个 without replacement 术语可能听起来令人困惑，因为我们已经在 rep_nums 中有重复的数字。但是，它的本质意思是，np.random.choice() 的输出不能包含例如多于两个 4's，因为 rep_nums 有两个 4's.
所以，问题是我们想要摆脱 np.repeat 部分，这可能是真正巨大数组的瓶颈。
提议的方法
查看 rep_nums 的输出，一个想法是生成 z = 7 个长度为 rep_nums 的唯一元素：
In [44]: np.random.choice(rep_nums.size,7,replace=False)
Out[44]: array([ 7, 2, 4, 10, 13, 8, 3])
这些数字表示该长度的索引。所以，我们只需要在 rep_nums 中寻找 bin（在 5 bin 之外），每个 7 数字都会进入其中。为此，我们可以使用 np.searchsorted。因此，我们将有一个实现来处理泛型 x、y，就像这样 -
# Get the intervals of those bins
intervals = data.astype(int).cumsum()
# Decide length of array if we had repeated with `np.repeat`
max_num = intervals[-1]
# Get unique numbers (indices in this case)
ids = np.random.choice(max_num,z,replace=False)
# Use searchsorted to get bin IDs and add in `x` offset
out = x+np.searchsorted(intervals,ids,'right')
运行时测试
函数：
def org_app(x,y,z,data):
rep_nums = np.repeat(range(x,y), data.astype('int64'))
out = np.random.choice(rep_nums, z,replace=False)
return out
def optimized_v1(x,y,z,data):
intervals = data.astype(int).cumsum()
max_num = intervals[-1]
ids = np.random.choice(max_num,z,replace=False)
out = x+np.searchsorted(intervals,ids,'right')
return out
全功能时间 -
In [79]: # Setup inputs
...: x = 100
...: y = 10010
...: z = 1000
...: data = np.random.randint(100,5000,(y-x)).astype(float)
...:
In [80]: %timeit org_app(x,y,z,data)
1 loop, best of 3: 7.17 s per loop
In [81]: %timeit optimized_v1(x,y,z,data)
1 loop, best of 3: 6.92 s per loop
看起来我们没有得到很好的加速。让我们深入挖掘一下，看看我们在替换 np.repeat!
上节省了多少
先关闭原来的方法-
In [82]: %timeit np.repeat(range(x,y), data.astype('int64'))
1 loop, best of 3: 227 ms per loop
让我们看看我们使用提议的方法对此有多大改进。因此，让我们在提议的方法中为 np.random.choice() 以外的所有内容计时 -
In [83]: intervals = data.astype(int).cumsum()
...: max_num = intervals[-1]
...: ids = np.random.choice(max_num,z,replace=False)
...: out = x+np.searchsorted(intervals,ids,'right')
...:
In [84]: %timeit data.astype(int).cumsum()
10000 loops, best of 3: 36.6 µs per loop
In [85]: %timeit intervals[-1]
10000000 loops, best of 3: 142 ns per loop
In [86]: %timeit x+np.searchsorted(intervals,ids,'right')
10000 loops, best of 3: 127 µs per loop
这比 np.repeat 的 227ms 好多了！！
因此，我们希望在非常大的数组中，删除 np.repeat 的好处会真正大放异彩，否则 np.random.choice() 本身看起来就是瓶颈。
为了完成，我还有一个替代实现。鉴于我们有 data，我们可以对每个 class:
使用超几何采样
- 反向计算
data.cumsum()
- 每 class 绘制
np.hypergeometric(data[pos], cumsum[pos]-data[pos], remain)
但是，当我们有很多 class 时，每个单元很少，这需要很长时间。
潜伏在问题中的是 multivariate hypergeometric distribution. In ，我实现了一个从该分布中抽取样本的函数。我怀疑它与答案中描述的解决方案@DiogoSantos 非常相似。 Diogo 说使用这种方法很慢，但我发现以下方法比 Divakar 的 optmized_v1.
更快
这是一个函数，它使用链接答案中的 sample(n, colors) 来实现与 Divakar 的函数具有相同签名的函数。
def hypergeom_version(x, y, z, data):
s = sample(z, data)
result = np.repeat(np.arange(x, y), s)
return result
(此 return 中的值按 排序 。如果您需要值随机排序，请在 [= 之前添加 np.random.shuffle(result) 38=] 语句。它不会显着改变执行时间。)
比较：
In [153]: x = 100
In [154]: y = 100100
In [155]: z = 10000
In [156]: data = np.random.randint(1, 125, (y-x)).astype(float)
Divakar 的 optimized_v1:
In [157]: %timeit optimized_v1(x, y, z, data)
1 loop, best of 3: 520 ms per loop
hypergeom_version:
In [158]: %timeit hypergeom_version(x, y, z, data)
1 loop, best of 3: 244 ms per loop
如果data中的值越大，相对性能就更好：
In [164]: data = np.random.randint(100, 500, (y-x)).astype(float)
In [165]: %timeit optimized_v1(x, y, z, data)
1 loop, best of 3: 2.91 s per loop
In [166]: %timeit hypergeom_version(x, y, z, data)
1 loop, best of 3: 246 ms per loop
我需要生成一个包含重复元素的大数组，我的代码是：
np.repeat(xrange(x,y), data)
但是，数据是一个类型为 float64 的 numpy 数组（但它表示整数，那里没有 2.1），我收到错误
TypeError: Cannot cast array data from dtype('float64') to dtype('int64') according to the rule 'safe'
示例：
In [35]: x
Out[35]: 26
In [36]: y
Out[36]: 50
In [37]: data
Out[37]:
array([ 3269., 106., 5533., 317., 1512., 208., 502., 919.,
406., 421., 1690., 2236., 705., 505., 230., 213.,
307., 1628., 4389., 1491., 355., 103., 854., 424.])
In [38]: np.repeat(xrange(x,y), data)
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-38-105860821359> in <module>()
----> 1 np.repeat(xrange(x,y), data)
/home/pcadmin/anaconda2/lib/python2.7/site-packages/numpy /core/fromnumeric.pyc in repeat(a, repeats, axis)
394 repeat = a.repeat
395 except AttributeError:
--> 396 return _wrapit(a, 'repeat', repeats, axis)
397 return repeat(repeats, axis)
398
/home/pcadmin/anaconda2/lib/python2.7/site-packages/numpy /core/fromnumeric.pyc in _wrapit(obj, method, *args, **kwds)
46 except AttributeError:
47 wrap = None
---> 48 result = getattr(asarray(obj), method)(*args, **kwds)
49 if wrap:
50 if not isinstance(result, mu.ndarray):
TypeError: Cannot cast array data from dtype('float64') to dtype('int64') according to the rule 'safe'
我把代码改成
解决了np.repeat(xrange(x,y), data.astype('int64'))
但是，现在这是我代码中最昂贵的行之一！！还有其他选择吗？
顺便说一句，我在里面用这个
np.random.choice(np.repeat(xrange(x,y), data.astype('int64')), z)
为了在不替换的情况下获得样本，其中 x 和 y 之间的整数的大小为 z，数据中给出了每个整数的数量。我想这也是最好的方法吧？
问题陈述
这个问题很有趣！只是为了让读者对这个问题有一个概念，而不是进入次要的数据转换问题，我们有一个值的范围，比方说 a = np.arange(5)，即
a = np.array([0,1,2,3,4])
现在，假设我们有另一个数组，其中列出了 a 中每个 5 数字的重复次数。所以，让这些成为：
reps = np.array([2,4,6,2,2])
接下来，我们将进行这些重复：
In [32]: rep_nums = np.repeat(a,reps)
In [33]: rep_nums
Out[33]: array([0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 4, 4])
最后，我们希望使用 np.random.choice() 从那些重复的数字中选择 z 个元素，而不用替换。
假设 z = 7 选择 7 个元素，那么对于 np.random.choice()，我们将有：
In [34]: np.random.choice(rep_nums,7,replace=False)
Out[34]: array([2, 4, 0, 2, 4, 1, 2])
现在，这个 without replacement 术语可能听起来令人困惑，因为我们已经在 rep_nums 中有重复的数字。但是，它的本质意思是，np.random.choice() 的输出不能包含例如多于两个 4's，因为 rep_nums 有两个 4's.
所以，问题是我们想要摆脱 np.repeat 部分，这可能是真正巨大数组的瓶颈。
提议的方法
查看 rep_nums 的输出，一个想法是生成 z = 7 个长度为 rep_nums 的唯一元素：
In [44]: np.random.choice(rep_nums.size,7,replace=False)
Out[44]: array([ 7, 2, 4, 10, 13, 8, 3])
这些数字表示该长度的索引。所以，我们只需要在 rep_nums 中寻找 bin（在 5 bin 之外），每个 7 数字都会进入其中。为此，我们可以使用 np.searchsorted。因此，我们将有一个实现来处理泛型 x、y，就像这样 -
# Get the intervals of those bins
intervals = data.astype(int).cumsum()
# Decide length of array if we had repeated with `np.repeat`
max_num = intervals[-1]
# Get unique numbers (indices in this case)
ids = np.random.choice(max_num,z,replace=False)
# Use searchsorted to get bin IDs and add in `x` offset
out = x+np.searchsorted(intervals,ids,'right')
运行时测试
函数：
def org_app(x,y,z,data):
rep_nums = np.repeat(range(x,y), data.astype('int64'))
out = np.random.choice(rep_nums, z,replace=False)
return out
def optimized_v1(x,y,z,data):
intervals = data.astype(int).cumsum()
max_num = intervals[-1]
ids = np.random.choice(max_num,z,replace=False)
out = x+np.searchsorted(intervals,ids,'right')
return out
全功能时间 -
In [79]: # Setup inputs
...: x = 100
...: y = 10010
...: z = 1000
...: data = np.random.randint(100,5000,(y-x)).astype(float)
...:
In [80]: %timeit org_app(x,y,z,data)
1 loop, best of 3: 7.17 s per loop
In [81]: %timeit optimized_v1(x,y,z,data)
1 loop, best of 3: 6.92 s per loop
看起来我们没有得到很好的加速。让我们深入挖掘一下，看看我们在替换 np.repeat!
先关闭原来的方法-
In [82]: %timeit np.repeat(range(x,y), data.astype('int64'))
1 loop, best of 3: 227 ms per loop
让我们看看我们使用提议的方法对此有多大改进。因此，让我们在提议的方法中为 np.random.choice() 以外的所有内容计时 -
In [83]: intervals = data.astype(int).cumsum()
...: max_num = intervals[-1]
...: ids = np.random.choice(max_num,z,replace=False)
...: out = x+np.searchsorted(intervals,ids,'right')
...:
In [84]: %timeit data.astype(int).cumsum()
10000 loops, best of 3: 36.6 µs per loop
In [85]: %timeit intervals[-1]
10000000 loops, best of 3: 142 ns per loop
In [86]: %timeit x+np.searchsorted(intervals,ids,'right')
10000 loops, best of 3: 127 µs per loop
这比 np.repeat 的 227ms 好多了！！
因此，我们希望在非常大的数组中，删除 np.repeat 的好处会真正大放异彩，否则 np.random.choice() 本身看起来就是瓶颈。
为了完成，我还有一个替代实现。鉴于我们有 data，我们可以对每个 class:
- 反向计算
data.cumsum() - 每 class 绘制
np.hypergeometric(data[pos], cumsum[pos]-data[pos], remain)
但是，当我们有很多 class 时，每个单元很少，这需要很长时间。
潜伏在问题中的是 multivariate hypergeometric distribution. In optmized_v1.
这是一个函数，它使用链接答案中的 sample(n, colors) 来实现与 Divakar 的函数具有相同签名的函数。
def hypergeom_version(x, y, z, data):
s = sample(z, data)
result = np.repeat(np.arange(x, y), s)
return result
(此 return 中的值按 排序 。如果您需要值随机排序，请在 [= 之前添加 np.random.shuffle(result) 38=] 语句。它不会显着改变执行时间。)
比较：
In [153]: x = 100
In [154]: y = 100100
In [155]: z = 10000
In [156]: data = np.random.randint(1, 125, (y-x)).astype(float)
Divakar 的 optimized_v1:
In [157]: %timeit optimized_v1(x, y, z, data)
1 loop, best of 3: 520 ms per loop
hypergeom_version:
In [158]: %timeit hypergeom_version(x, y, z, data)
1 loop, best of 3: 244 ms per loop
如果data中的值越大，相对性能就更好：
In [164]: data = np.random.randint(100, 500, (y-x)).astype(float)
In [165]: %timeit optimized_v1(x, y, z, data)
1 loop, best of 3: 2.91 s per loop
In [166]: %timeit hypergeom_version(x, y, z, data)
1 loop, best of 3: 246 ms per loop